Luna Flow
English

immut tutorial

This tutorial shows how to work with immutable matrices and vectors: build them, update them without losing earlier versions, compute exact integer results such as powers and determinants, and use lazy matrices for structured data. The reasoning behind the package is in the immut design.

Quick start

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 returns a new matrix; a keeps its value.

Everyday tasks

Keep a history of edits

Because every update returns a new value, a history is just an array of matrices that share most of their storage:

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

Count paths with a matrix power

Entry (i,j)(i, j) of AkA^k counts the walks of length kk from ii to jj in the graph with adjacency matrix AA. Integer powers are exact:

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

Exact determinants

determinant uses fraction-free elimination, so integer matrices get exact results, and BigInt matrices are exact at any size. Scaling every entry by 101210^{12} scales the 5×55 \times 5 determinant by 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",
  )
}

Build matrices from vectors

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

Describe a structured matrix lazily

A MatrixFn stores a rule instead of entries. It is a good fit for matrices defined by a formula, such as a band matrix, when you read only part of them:

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

No 1000×10001000 \times 1000 array is ever allocated.

Going further

Generic code. The concrete type does not implement the algebra traits. Wrap it in @default.ImmutableDenseMatrix to pass it to functions bounded by MatMulMatrix (see the backends/default tutorial).

Converting to mutable. For numerical work such as inverses, decompositions or statistics, convert with the container adapters and use mutable.

Your own scalar type. Any type implementing the luna-generic traits works: a rational or modular-integer type with Num and Div gets exact determinants through Bareiss elimination, and any Semiring gets pow.

Performance. Reads and single-entry updates cost O(log⁡32N)O(\log_{32} N). For a long series of edits on a large matrix, a @mutable.Matrix working buffer is faster; freeze it into an immut value at the end.

Common pitfalls

  • Integer overflow. Int arithmetic wraps. pow is then exact modulo 2322^{32}, which is rarely what you want; use Int64 or BigInt for large values.
  • Assigning through indexing. m[r][c] = x does not exist here; use m.set(r, c, x) and keep the result.
  • Subtracting vectors. Vector has no - operator; write u + -v.
  • Reading lazy powers. Every read of a MatrixFn power recomputes the nested products. Materialize first with Matrix::make(r, c, (i, j) => f[i][j]) if you read many entries.
  • Operators abort. +, -, * abort on shape mismatch; use matmul when shapes come from input.

Next steps