L11 — Matrices

Goal: drive the matrix domain from Scala without being surprised by it. The mathematics is the ordinary kind; what needs a lesson is the shape of the API — two carriers where you expect one, a parser that decides on structure before anything is evaluated, and results that come back as Left although nothing went wrong.

import it.grypho.scala.leonardo.core.*
import it.grypho.scala.leonardo.matrix.*
import it.grypho.scala.leonardo.parser.Parser

val env = Environment()
def parse(s: String): _Expression = Parser.parse(s).get

Two carriers

A matrix of numbers is a _MatrixValue: a dense, row-major array of Doubles, the thing the numeric kernels (multiply, LU, QR, eigenvalues) work on.

val dense = _MatrixValue(2, 2, Array(1.0, 2.0, 3.0, 4.0))
// dense: _MatrixValue = [[1.0, 2.0], [3.0, 4.0]]
dense(1, 0)
// res0: Double = 3.0

A matrix of expressions is a _Matrix: its cells are any _Expression — a free variable, a function, another matrix. It is what the parser builds from a literal:

val symbolic = parse("[[x, 1], [0, sin(y)]]")
// symbolic: _Expression = _Matrix(
//   rows = 2,
//   cols = 2,
//   elems = Vector(
//     _Variable("x"),
//     _Number(1.0),
//     _Number(0.0),
//     Sin(_Variable("y"))
//   )
// )

The two meet in eval. A _Matrix whose cells all reduce to numbers collapses to a _MatrixValue; one that keeps a free name stays a _Matrix, with every cell reduced as far as it goes:

parse("[[1, 2], [3, 1 + 3]]").eval(env)
// res1: Either[_Expression, _Value] = Right([[1.0, 2.0], [3.0, 4.0]])
symbolic.eval(env)
// res2: Either[_Expression, _Value] = Left(
//   _Matrix(
//     rows = 2,
//     cols = 2,
//     elems = Vector(
//       _Variable("x"),
//       _Number(1.0),
//       _Number(0.0),
//       Sin(_Variable("y"))
//     )
//   )
// )
symbolic.eval(env.withBinding("x", _Number(2.0)).withBinding("y", _Number(0.0)))
// res3: Either[_Expression, _Value] = Right([[2.0, 1.0], [0.0, 0.0]])

So your code reads a matrix result from both sides of the Either — Right(_MatrixValue) when it is numeric, Left(_Matrix) when it is not — and the Left here is a complete answer, not a failure. L2’s rule, met again: Left means still symbolic, and a matrix with a free name in it is exactly that.

There is one more case worth knowing before it surprises you. An exact matrix never collapses, even when every cell is a number: the dense carrier holds Doubles, and collapsing would round the exact cells away. So in exact mode a numeric matrix is a Left(_Matrix) of rationals:

Parser.parse("[[1, 2], [3, 4]] * [[1/3, 0], [0, 1]]", Some(30)).get.eval(env)
// res4: Either[_Expression, _Value] = Left(
//   _Matrix(rows = 2, cols = 2, elems = Vector(1/3, 2, 1, 4))
// )

The parser decides on structure

* means two different things — a matrix product and a scalar product — and the parser has to choose a node before anything is evaluated, when all it can see is the shape of each operand. A matrix literal, or a matrix operation, looks like a matrix:

parse("[[1, 2], [3, 4]] * [[0, 1], [1, 0]]")
// res5: _Expression = MatProduct(
//   a = _Matrix(
//     rows = 2,
//     cols = 2,
//     elems = Vector(_Number(1.0), _Number(2.0), _Number(3.0), _Number(4.0))
//   ),
//   b = _Matrix(
//     rows = 2,
//     cols = 2,
//     elems = Vector(_Number(0.0), _Number(1.0), _Number(1.0), _Number(0.0))
//   )
// )

Two plain names do not. A * B is an ordinary scalar Product — and it still multiplies two matrices correctly, because the scalar node’s evaluation handles dense matrix values when it meets them:

val ab = parse("A * B")
// ab: _Expression = Product(a = _Variable("A"), b = _Variable("B"))
ab.eval(env.withBinding("A", dense).withBinding("B", _MatrixValue(2, 2, Array(0.0, 1.0, 1.0, 0.0))))
// res6: Either[_Expression, _Value] = Right([[2.0, 1.0], [4.0, 3.0]])

The case that is easy to get wrong is a literal matrix times a name:

parse("[[1, 2], [3, 4]] * y")
// res7: _Expression = MatProduct(
//   a = _Matrix(
//     rows = 2,
//     cols = 2,
//     elems = Vector(_Number(1.0), _Number(2.0), _Number(3.0), _Number(4.0))
//   ),
//   b = _Variable("y")
// )

It is a MatProduct, not a scaling of the matrix by y, and the reason is the whole point of this section: y might be bound to a matrix by the time the expression is evaluated, and a matrix product does not commute — turning M·y into y·M would be a different answer. Only a literal number on the right is known to commute, so only that becomes a MatScale:

parse("[[1, 2], [3, 4]] * 2")
// res8: _Expression = MatScale(
//   k = _Number(2.0),
//   m = _Matrix(
//     rows = 2,
//     cols = 2,
//     elems = Vector(_Number(1.0), _Number(2.0), _Number(3.0), _Number(4.0))
//   )
// )

When you build trees yourself (L4), make the same choice: MatProduct whenever an operand could be a matrix, keeping the order you were given.

Results that are rows of matrices

A decomposition answers several matrices at once, and the library returns them as a 1×n _Matrix — a row whose cells are matrices. LU, for instance, is [[L, U, P]] with P·A = L·U:

val lu = _LUDecomposition(dense).eval(env)
// lu: Either[_Expression, _Value] = Left(
//   _Matrix(
//     rows = 1,
//     cols = 3,
//     elems = Vector(
//       [[1.0, 0.0], [0.33333, 1.0]],
//       [[3.0, 4.0], [0.0, 0.66667]],
//       [[0.0, 1.0], [1.0, 0.0]]
//     )
//   )
// )

The row is a Left, because a matrix whose cells are matrices is not a dense matrix of numbers — the same reason as before. Read it by matching the carrier and indexing the cells:

def factors(result: Either[_Expression, _Value]): Option[Vector[_Expression]] = result match
  case Left(row: _Matrix) if row.rows == 1 => Some(row.elems)
  case _                                   => None

factors(lu).map(_.size)
// res9: Option[Int] = Some(3)
factors(lu).map(fs => fs(1))
// res10: Option[_Expression] = Some([[3.0, 4.0], [0.0, 0.66667]])

At the prompt the same row is indexed with at(r, 1, k), and tuple assignment L, U, P := lu(A) binds the three names at once — both are spellings of reading this row.

What stays as written

A matrix operation that has no answer stays an unevaluated node rather than throwing, so the tree you get back tells you what was asked:

Inverse(_MatrixValue(2, 2, Array(1.0, 2.0, 2.0, 4.0))).eval(env)
// res11: Either[_Expression, _Value] = Left(Inverse([[1.0, 2.0], [2.0, 4.0]]))
parse("[[1, 2]] * [[1, 2]]").eval(env)
// res12: Either[_Expression, _Value] = Left(
//   MatProduct(
//     a = _Matrix(rows = 1, cols = 2, elems = Vector(_Number(1.0), _Number(2.0))),
//     b = _Matrix(rows = 1, cols = 2, elems = Vector(_Number(1.0), _Number(2.0)))
//   )
// )

The first is singular, the second does not conform (a 1×2 times a 1×2). Both come back as Left of the node, which in L9’s catalogue is a refusal. Telling it apart from a symbolic answer is your code’s job, and the test is the node’s type: a Left(_Matrix) is a matrix with symbols in it, a Left of an operation node is an operation that did not happen.

Try it yourself

Write solve2(a, b) that solves A·x = b for a dense A and column b through Inverse and MatProduct, and returns Option[_MatrixValue] — None when A is singular.

Answers

```scala def solve2(a: _MatrixValue, b: _MatrixValue): Option[_MatrixValue] = MatProduct(Inverse(a), b).eval(env) match case Right(x: _MatrixValue) => Some(x) case _ => None ``` ```scala solve2(_MatrixValue(2, 2, Array(2.0, 1.0, 1.0, 3.0)), _MatrixValue(2, 1, Array(3.0, 5.0))) // res13: Option[_MatrixValue] = Some([[0.8], [1.4]]) solve2(_MatrixValue(2, 2, Array(1.0, 2.0, 2.0, 4.0)), _MatrixValue(2, 1, Array(3.0, 5.0))) // res14: Option[_MatrixValue] = None ``` The singular system comes back `None` because the inner `Inverse` stayed unevaluated, so the product had nothing to multiply and stayed unevaluated too: a refusal propagates outward through the tree without any code of yours checking for it, and the one `match` at the end is the only place that has to know.

What you learned

• Two carriers: _MatrixValue for numbers, _Matrix for expressions. eval collapses the second into the first when every cell is a number and the cells are not exact.
• A matrix result can be a Left — a matrix with free names, an exact matrix, or a row of matrices from a decomposition — and that Left is an answer.
• The parser chooses MatProduct from shape alone, and keeps the order, because a name may be a matrix by the time it is evaluated. Build your own trees the same way.
• A refused operation stays as its node; the node’s type is what distinguishes it from a symbolic result.

Where to go next

• The Matrix reference lists every operation, decomposition and constructor in the language these examples use.
• L9 — Refusals is the catalogue the last section added to.

---

Previous: L10 — A small application, end to end · Next: — the most recent lesson · All lessons