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|")
}

行列のべき乗で経路を数える

AkA^k の要素 (i,j)(i, j) は、隣接行列 AA を持つグラフにおける 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 を実装していません。MatMulMatrix で境界付けられた関数に渡すには @default.ImmutableDenseMatrix でラップしてください(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 を使ってください。

次のステップ