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) のビューです。2 つのビューの積は行列カーネルを再利用します。独立した行列が必要な場合は transpose() か materialize() を使ってください。

ホットループでの検査なし形式。 前提条件を確立した後（たとえば自分のコードで構築した正方行列）であれば、unchecked_* メソッドは検証を省略します。検査していない形状に unchecked_matmul を使ってはいけません。検証を行わず、誤った結果を返すことがあります。

スケーリング。 ルーチンは 10−1110^{-11} という絶対しきい値で「ゼロ」を判定します。先にデータを 1 程度の大きさにスケーリングしてください。そうしないと、小さいだけの正則な行列が特異と報告されます。設計ページ を参照してください。

汎用コード。 algebra の trait で境界付けられた関数に渡すには、行列を @default.DenseMatrix でラップしてください。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 を返します。

次のステップ