Luna Flow
简体中文

mutable 教程

本教程演示如何把 @mutable.Matrix 和 @mutable.Vector 用作工作缓冲区：原地构造和编辑它们，通过行、列和转置视图进行操作，并在妥善处理错误的前提下运行数值例程（逆矩阵、行列式、Cholesky、特征值、统计）。算法及其精度在 mutable 设计中解释。

快速上手

moon add Luna-Flow/linear-algebra@0.5.0
///|
import {
  "Luna-Flow/linear-algebra/mutable",
}
///|
test "first mutable matrix" {
  let m = @mutable.Matrix::from_2d_array([[1.0, 2.0], [3.0, 4.0]])
  m.set(0, 1, 9.0)
  inspect(m, content="|1, 9|\n|3, 4|")
  inspect(m.determinant().unwrap(), content="-23")
}

set 原地修改 m；determinant 返回 Result，因为它只对方阵有定义。

日常任务

通过视图编辑行和列

///|
test "normalize each row to sum one" {
  let m = @mutable.Matrix::from_2d_array([[1.0, 3.0], [2.0, 2.0], [0.0, 5.0]])
  for i in 0..<m.row() {
    let row = m.row_view(i)
    let mut sum = 0.0
    row.each(x => sum = sum + x)
    row.map_inplace(x => x / sum)
  }
  inspect(m, content="|0.25, 0.75|\n|0.5, 0.5|\n|0, 1|")
  let first_col = m.col_view(0)
  first_col[2] = 1.0
  inspect(m.get(2, 0), content="1")
}

求解小型线性方程组

没有针对右端项的公开求解器；对于小型方程组，请乘以逆矩阵，并处理奇异的情形：

///|
fn mut_tut_solve(
  a : @mutable.Matrix[Double],
  b : @mutable.Vector[Double],
) -> Result[@mutable.Vector[Double], @la_error.LinearAlgebraError] {
  let inv = match a.inverse() {
    Ok(inv) => inv
    Err(e) => return Err(e)
  }
  inv.mul_vec(b)
}

///|
test "solve 2x + y = 5, x + 3y = 10" {
  let a = @mutable.Matrix::from_2d_array([[2.0, 1.0], [1.0, 3.0]])
  let x = mut_tut_solve(a, @mutable.Vector::from_array([5.0, 10.0])).unwrap()
  inspect((x[0] - 1.0).abs() < 1.0e-12, content="true")
  inspect((x[1] - 3.0).abs() < 1.0e-12, content="true")
  let singular = @mutable.Matrix::from_2d_array([[1.0, 2.0], [2.0, 4.0]])
  match mut_tut_solve(singular, @mutable.Vector::from_array([1.0, 2.0])) {
    Err(e) => inspect(e.is_singular_matrix(), content="true")
    Ok(_) => fail("singular systems have no unique solution")
  }
}

如上所示，请带容差比较浮点结果，而不是使用 ==。

检查正定性并进行分解

///|
test "covariance-like matrix" {
  let c = @mutable.Matrix::from_2d_array([
    [4.0, 2.0, 0.0],
    [2.0, 5.0, 1.0],
    [0.0, 1.0, 3.0],
  ])
  inspect(c.is_symmetric(), content="true")
  match c.cholesky_decomposition() {
    Some(l) => {
      let back = l * l.transpose()
      let mut err = 0.0
      back.each_row_col((i, j, x) => err = err + (x - c.get(i, j)).abs())
      inspect(err < 1.0e-12, content="true")
    }
    None => fail("c is positive definite")
  }
  let indefinite = @mutable.Matrix::from_2d_array([[1.0, 2.0], [2.0, 1.0]])
  inspect(indefinite.cholesky_decomposition() is None, content="true")
}

对称矩阵的特征值

///|
test "vibration modes of a three-mass chain" {
  let k = @mutable.Matrix::from_2d_array([
    [2.0, -1.0, 0.0],
    [-1.0, 2.0, -1.0],
    [0.0, -1.0, 2.0],
  ])
  let (values, vectors) = k.eigen()
  let sorted = values.iter().to_array()
  sorted.sort()
  let expected = [2.0 - 1.4142135623730951, 2.0, 2.0 + 1.4142135623730951]
  for i in 0..<3 {
    inspect((sorted[i] - expected[i]).abs() < 1.0e-12, content="true")
  }
  let check = vectors.transpose() * vectors
  inspect((check.get(0, 0) - 1.0).abs() < 1.0e-12, content="true")
  inspect(check.get(0, 1).abs() < 1.0e-12, content="true")
}

返回的特征值未排序；当顺序重要时，请自行排序。

汇总统计

///|
test "statistics of a measurement grid" {
  let grid = @mutable.Matrix::from_2d_array([
    [2.0, 4.0],
    [4.0, 4.0],
    [5.0, 5.0],
    [7.0, 9.0],
  ])
  inspect(grid.mean().unwrap(), content="5")
  inspect(grid.std_dev().unwrap(), content="2")
  inspect(grid.min_element().unwrap(), content="2")
  let empty : @mutable.Matrix[Double] = @mutable.Matrix::new(0, 2, 0.0)
  inspect(empty.mean() is Err(_), content="true")
}

进一步了解

避免复制。 Matrix::from_array 和 Vector::from_array 直接接管你传入的数组。这很快，但意味着数组和矩阵会一起变化；需要独立时请调用 copy()。

不复制的转置。 to_transpose() 是 O(1)O(1) 的视图；两个视图的乘积会复用矩阵内核。需要独立的矩阵时，请使用 transpose() 或 materialize()。

在热循环中使用非受检形式。 在确立前置条件之后（例如由你自己的代码构造的方阵），unchecked_* 方法会跳过校验。切勿在未检查过的形状上使用 unchecked_matmul：它不做校验，可能返回错误的结果。

缩放。 这些例程以 10−1110^{-11} 的绝对阈值判断“零”。请先把数据缩放到量级为一；否则很小但正则的矩阵会被报告为奇异。参见设计页面。

泛型代码。 把矩阵包装为 @default.DenseMatrix，即可传给以 algebra trait 为约束的函数；用 inner() 取回。

常见陷阱

  • 忘记 reduce_row_elimination 会修改矩阵。 它原地变换接收者并将其返回。请先调用 copy()。
  • 期望 2×22 \times 2 特征向量已归一化。 对 2×22 \times 2 输入，特征向量列未归一化；更大的输入返回标准正交的列。
  • 在整数上使用数值例程。 determinant、inverse、rank、eigen 等需要 Tolerance，而它只为 Float 和 Double 提供。要得到精确的整数行列式，请使用 @immut.Matrix::determinant。
  • Float 精度。 容差远低于 Float 的精度；接近奇异的 Float 矩阵不会被检测出来。请优先使用 Double。
  • 在对称谱上使用幂法。 模相等的特征值 ±λ\pm\lambda 会阻碍收敛；此时 power_method 返回 None。

后续步骤