algebra 教程

本教程演示如何针对 algebra trait 只编写一次线性代数辅助函数,并在任何实现了这些 trait 的矩阵或向量类型上运行:可以是仓库的稠密包装类型,也可以是你自己的类型。你需要了解 MoonBit 泛型;数学部分尽量从简,详见 algebra 设计。

快速上手

添加模块以及你使用的包:

moon add Luna-Flow/linear-algebra@0.5.0
///|
import {
  "Luna-Flow/linear-algebra/algebra",
  "Luna-Flow/linear-algebra/backends/default",
}

最小的有用程序通过 trait 查询矩阵的形状,因此它适用于每个 MatrixShape 类型:

///|
fn[M : @algebra.MatrixShape] alg_tut_describe(matrix : M) -> String {
  let (rows, cols) = @algebra.MatrixShape::shape(matrix)
  "\{rows}x\{cols}"
}

///|
test "describe a dense matrix" {
  let m = @default.DenseMatrix::from_2d_array([[1, 2, 3], [4, 5, 6]])
  inspect(alg_tut_describe(m), content="2x3")
}

输出为 2x3。

日常任务

为任意加法向量类型计算残差

残差 r=b−Axr = b - A x 只需要向量减法。要求 AdditiveVector,该辅助函数就能同时接受可变和不可变的稠密向量:

///|
fn[V : @algebra.AdditiveVector] alg_tut_residual(
  observed : V,
  predicted : V,
) -> V {
  observed - predicted
}

///|
test "residual over two vector representations" {
  let b = @default.DenseVector::from_array([3, 5, 7])
  let p = @default.DenseVector::from_array([1, 5, 9])
  let r = alg_tut_residual(b, p)
  inspect(r.inner(), content="|2, 0, -2|")
  let bi = @default.ImmutableDenseVector::from_array([3, 5, 7])
  let pi = @default.ImmutableDenseVector::from_array([1, 5, 9])
  inspect(alg_tut_residual(bi, pi).inner(), content="|2, 0, -2|")
}

构造 Gram 矩阵 ATAA^{\mathsf T} A

Gram 矩阵需要转置和矩阵乘积,这恰好就是 MatMulMatrix:

///|
fn[M : @algebra.MatMulMatrix] alg_tut_gram(a : M) -> M {
  @algebra.TransposeMatrix::transpose(a) * a
}

///|
test "gram matrix of a 3x2 design matrix" {
  let a = @default.ImmutableDenseMatrix::from_2d_array([[1, 0], [1, 1], [1, 2]])
  let g = alg_tut_gram(a)
  debug_inspect(@algebra.MatrixShape::shape(g), content="(2, 2)")
  inspect(g.inner(), content="|3, 3|\n|3, 5|")
}

无论 AA 是什么形状,ATAA^{\mathsf T}A 总有定义:m×nm \times n 矩阵得到 n×nn \times n 的结果。

相乘之前检查可复合性

在运行时确定形状的矩阵上,* 是部分运算,稠密包装类型在不匹配时会中止。当形状来自数据时,请先用 MatrixShape 检查:

///|
fn[M : @algebra.MatMulMatrix] alg_tut_try_product(a : M, b : M) -> M? {
  let (_, inner_a) = @algebra.MatrixShape::shape(a)
  let (inner_b, _) = @algebra.MatrixShape::shape(b)
  if inner_a == inner_b {
    Some(a * b)
  } else {
    None
  }
}

///|
test "product only for composable shapes" {
  let a = @default.DenseMatrix::from_2d_array([[1, 2]])
  let b = @default.DenseMatrix::from_2d_array([[3], [4]])
  debug_inspect(
    alg_tut_try_product(a, b).map(m => m.inner().get(0, 0)),
    content="Some(11)",
  )
  inspect(alg_tut_try_product(a, a) is None, content="true")
}

使用自己的定长矩阵

2×22 \times 2 矩阵类型的乘积是全函数,因此它可以实现每一层矩阵 trait,而不带任何运行时前置条件。先实现运算符 trait,再声明各层:

///|
struct AlgTutMat2 {
  a : Int
  b : Int
  c : Int
  d : Int
}

///|
impl @algebra.MatrixShape for AlgTutMat2 with fn shape(_) {
  (2, 2)
}

///|
impl @algebra.TransposeMatrix for AlgTutMat2 with fn transpose(m) {
  { a: m.a, b: m.c, c: m.b, d: m.d, }
}

///|
impl Add for AlgTutMat2 with fn add(x, y) {
  { a: x.a + y.a, b: x.b + y.b, c: x.c + y.c, d: x.d + y.d, }
}

///|
impl Neg for AlgTutMat2 with fn neg(x) {
  { a: -x.a, b: -x.b, c: -x.c, d: -x.d, }
}

///|
impl Sub for AlgTutMat2 with fn sub(x, y) {
  x + -y
}

///|
impl Mul for AlgTutMat2 with fn mul(x, y) {
  {
    a: x.a * y.a + x.b * y.c,
    b: x.a * y.b + x.b * y.d,
    c: x.c * y.a + x.d * y.c,
    d: x.c * y.b + x.d * y.d,
  }
}

///|
impl @algebra.AdditiveMatrix for AlgTutMat2

///|
impl @algebra.MatMulMatrix for AlgTutMat2

///|
test "the generic gram helper runs on a custom type" {
  let m : AlgTutMat2 = { a: 1, b: 2, c: 3, d: 4, }
  let g = alg_tut_gram(m)
  debug_inspect((g.a, g.b, g.c, g.d), content="(10, 14, 14, 20)")
}

上一个任务中的辅助函数 alg_tut_gram 无需任何修改。

进一步了解

使用 backends/default 的泛型辅助函数。 @default.shape_of、@default.transpose 和 @default.matmul 是这三种运算带 trait 约束的版本;当你需要函数值而不是带 trait 限定的调用时,它们很方便。

与 container 结合使用。 algebra trait 描述整体对象运算,从不暴露元素。当算法还需要读取或构造单个元素时,请在 trait 约束之外配合使用 container 的操作字典;这两层相互独立。

标量要求留在具体类型上。 这些 trait 不涉及标量类型。如果算法需要标量乘法或点积,请接收具体类型(T : AddMonoid + Mul 的 @default.DenseVector[T]),或把标量运算作为函数传入。

受检乘积。 当失败必须作为值返回时,请转换为具体类型并使用受检方法,例如返回 Result[_, LinearAlgebraError] 的 @immut.Matrix::matmul。参见错误教程。

常见陷阱

  • 在泛型代码中用点语法调用 trait 方法。 在 fn[M : @algebra.MatMulMatrix] 内部,应写作 @algebra.TransposeMatrix::transpose(m)。在 MoonBit 0.10 中,对类型参数用点语法调用父 trait 方法已被弃用。
  • 用 Hadamard * 实现 MatMulMatrix。 MatMulMatrix 上的 * 必须是矩阵乘积;* 为逐元素乘法的类型应止步于 AdditiveMatrix。
  • 期望对所有标量都有 (AB)^T = B^T A^T。 它需要可交换的标量;四元数矩阵就违反了这一点。
  • 精确比较浮点结果。 Double 矩阵的乘积只在舍入误差范围内一致。请带容差比较(见 arithmetic 教程)。

后续步骤