Luna Flow
简体中文

immut 教程

本教程演示如何使用不可变矩阵和向量：构造它们，在不丢失早期版本的情况下更新它们，计算幂和行列式等精确的整数结果，以及用惰性矩阵表示结构化数据。该包背后的道理见 immut 设计。

快速上手

moon add Luna-Flow/linear-algebra@0.5.0
///|
import {
  "Luna-Flow/linear-algebra/immut",
}
///|
test "first immutable matrix" {
  let a = @immut.Matrix::from_2d_array([[1, 2], [3, 4]])
  let b = a.set(0, 0, 10)
  inspect(a, content="|1, 2|\n|3, 4|")
  inspect(b, content="|10, 2|\n|3, 4|")
  inspect(a * b, content="|16, 10|\n|42, 22|")
}

set 返回新矩阵；a 保持原值。

日常任务

保留编辑历史

由于每次更新都返回新值，历史记录就是一个矩阵数组，这些矩阵共享大部分存储：

///|
test "undo by keeping old versions" {
  let history = [@immut.Matrix::new(2, 2, 0)]
  for step in 0..<3 {
    let current = history[history.length() - 1]
    history.push(current.set(step % 2, step / 2, step + 1))
  }
  inspect(history[3], content="|1, 3|\n|2, 0|")
  inspect(history[1], content="|1, 0|\n|0, 0|")
}

用矩阵幂计数路径

在邻接矩阵为 AA 的图中，AkA^k 的元素 (i,j)(i, j) 统计从 ii 到 jj 的长度为 kk 的通路数。整数幂是精确的：

///|
test "walks in a triangle graph" {
  let triangle = @immut.Matrix::from_2d_array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
  let walks = triangle.pow(4).unwrap()
  inspect(walks[0][0], content="6")
  inspect(walks[0][1], content="5")
  match triangle.pow(-1) {
    Err(e) => inspect(e.is_negative_exponent(), content="true")
    Ok(_) => fail("negative powers are rejected")
  }
}

精确行列式

determinant 使用无分数消元，因此整数矩阵得到精确结果，BigInt 矩阵在任意规模下都是精确的。把每个元素放大 101210^{12} 倍，5×55 \times 5 行列式就放大 106010^{60} 倍：

///|
test "exact determinant of an integer matrix" {
  let m = @immut.Matrix::from_2d_array([
    [2, 0, 1, 3, 1],
    [1, 3, 2, 0, 4],
    [0, 1, 4, 1, 2],
    [3, 2, 0, 5, 1],
    [1, 0, 2, 1, 3],
  ])
  inspect(m.determinant().unwrap(), content="62")
  let big = m.map(x => BigInt::from_int(x) * 1000000000000N)
  inspect(
    big.determinant().unwrap(),
    content="62000000000000000000000000000000000000000000000000000000000000",
  )
}

由向量构造矩阵

///|
test "outer products and diagonals" {
  let u = @immut.Vector::from_array([1, 2, 3])
  let ones = @immut.Vector::make(3, 1)
  let rank_one = u.tensor_product(ones)
  inspect(rank_one, content="|1, 1, 1|\n|2, 2, 2|\n|3, 3, 3|")
  let d = u.scaled_matrix()
  inspect(d * rank_one, content="|1, 1, 1|\n|4, 4, 4|\n|9, 9, 9|")
  inspect(u.to_row_matrix() * u.to_col_matrix(), content="|14|")
}

惰性地描述结构化矩阵

MatrixFn 存储的是规则而不是元素。它很适合由公式定义的矩阵（例如带状矩阵），尤其是只读取其中一部分时：

///|
test "a lazy tridiagonal matrix" {
  let n = 1000
  let band = @immut.MatrixFn::make(n, n, (i, j) => {
    if i == j {
      2
    } else if i - j == 1 || j - i == 1 {
      -1
    } else {
      0
    }
  })
  inspect(band[500][499], content="-1")
  inspect(band[500][700], content="0")
  let small = @immut.MatrixFn::make(3, 3, (i, j) => band[i][j])
  inspect(small, content="|2, -1, 0|\n|-1, 2, -1|\n|0, -1, 2|")
  inspect(small.determinant(), content="4")
}

从未分配过 1000×10001000 \times 1000 的数组。

进一步了解

泛型代码。 该具体类型没有实现 algebra trait。将其包装为 @default.ImmutableDenseMatrix，即可传给以 MatMulMatrix 为约束的函数（见 backends/default 教程）。

转换为 mutable。 对于求逆、分解或统计等数值工作，请用 container 适配器转换，再使用 mutable。

你自己的标量类型。 任何实现了 luna-generic trait 的类型都可以使用：带有 Num 和 Div 的有理数或模整数类型可通过 Bareiss 消元得到精确行列式，任何 Semiring 都可以使用 pow。

性能。 读取和单元素更新的开销为 O(log⁡32N)O(\log_{32} N)。对大矩阵做一长串编辑时，用 @mutable.Matrix 作为工作缓冲区更快；最后再把它冻结为 immut 值。

常见陷阱

  • 整数溢出。 Int 算术会回绕。此时 pow 在模 2322^{32} 意义下精确，但这很少是你想要的；对大数值请使用 Int64 或 BigInt。
  • 通过索引赋值。 这里不存在 m[r][c] = x；请使用 m.set(r, c, x) 并保留结果。
  • 向量相减。 Vector 没有 - 运算符；请写作 u + -v。
  • 读取惰性幂。 每次读取 MatrixFn 幂都会重新计算嵌套的乘积。如果要读取很多元素，请先用 Matrix::make(r, c, (i, j) => f[i][j]) 物化。
  • 运算符会中止。 +、-、* 在形状不匹配时中止；当形状来自输入时，请使用 matmul。

后续步骤