Matrices

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

val env = new Environment()

Constructing a matrix

_Matrix.ofRows accepts row vectors as Vector[_Expression] varargs:

// 2×2 identity matrix
val I2 = _Matrix.ofRows(
  Vector(_Number(1.0), _Number(0.0)),
  Vector(_Number(0.0), _Number(1.0))
)

I2.toString
// res0: String = "[[1.0, 0.0], [0.0, 1.0]]"

The parser accepts the same structure as a [[…], […]] literal:

Parser.parse("[[1, 0], [0, 1]]").get.toString
// res1: String = "[[1.0, 0.0], [0.0, 1.0]]"

A row vector uses double brackets: [[1, 2, 3]].

Element types

Matrix elements can be arbitrary expressions — variables, functions, even functionals. The matrix is symbolic until every element reduces to a number:

val x = _Variable("x")
val symM = _Matrix.ofRows(
  Vector(x,                        Sin(x)),
  Vector(Power(x, _Number(2.0)),   _Number(1.0))
)

symM.toString
// res2: String = "[[x, sin(x)], [(x ^ 2.0), 1.0]]"

Evaluate with x = 1:

val envX = new Environment(5, Map("x" -> _Number(1.0)))
// envX: Environment = it.grypho.scala.leonardo.core.Environment@2e698805
symM.eval(envX)
// res3: Either[_Expression, _Value] = Right(
//   value = [[1.0, 0.84147], [1.0, 1.0]]
// )

Arithmetic

Standard +, *, scalar multiplication, transpose, determinant and inverse work on both symbolic and concrete matrices:

val A = Parser.parse("[[1, 2], [3, 4]]").get
val B = Parser.parse("[[5, 6], [7, 8]]").get

// A + B
MatSum(A, B).eval(env)
// res4: Either[_Expression, _Value] = Right(
//   value = [[6.0, 8.0], [10.0, 12.0]]
// )

// A * B
MatProduct(A, B).eval(env)
// res5: Either[_Expression, _Value] = Right(
//   value = [[19.0, 22.0], [43.0, 50.0]]
// )

// 2 * A
MatScale(_Number(2.0), A).eval(env)
// res6: Either[_Expression, _Value] = Right(value = [[2.0, 4.0], [6.0, 8.0]])

// Aᵀ
Transpose(A).eval(env)
// res7: Either[_Expression, _Value] = Right(value = [[1.0, 3.0], [2.0, 4.0]])

// det(A) — a scalar result (LU decomposition with partial pivoting on dense values)
Determinant(A).eval(env)
// res8: Either[_Expression, _Value] = Right(value = _Number(d = -2.0))

// inv(A) — the inverse matrix (Gauss–Jordan elimination); 1/A and M/N parse to this too
Inverse(A).eval(env)
// res9: Either[_Expression, _Value] = Right(
//   value = [[-2.0, 1.0], [1.5, -0.5]]
// )

det(A) reduces to a scalar _Number; a non-square or (for the inverse) singular matrix has no result and stays symbolic, the same x/0 contract used elsewhere. The reciprocal spellings 1 / A and M / N (= M · N⁻¹) route through the same Inverse node. Small symbolic matrices expand by cofactors (determinant) and adjugate/det (inverse). Algorithm references: LU decomposition, Gauss–Jordan elimination.

Dense evaluation

When all elements reduce to numbers the result is a _MatrixValue — a dense row-major Array[Double]. The multiply kernel is block-tiled for cache efficiency and runs in parallel above a 2¹⁶-element work threshold:

MatProduct(A, B).eval(env) match {
  case Right(m: _MatrixValue) => m.toVector
  case other                  => other
}
// res10: Object & Equals = Vector(19.0, 22.0, 43.0, 50.0)

Calculus on matrices

derive, simplify, and expand distribute element-wise over _Matrix nodes (the _ElementWise marker enables this without domain-specific cases):

val exprM = _Matrix.ofRows(
  Vector(Power(x, _Number(2.0)),  Sin(x)),
  Vector(Exp(x),                  _Number(1.0))
)

derive(exprM, x).toString
// res11: String = "[[(2.0 * x), cos(x)], [exp(x), 0.0]]"

Matrix products do not distribute automatically under differentiation — they need the product rule applied explicitly.

Functions on matrices

Scalar functions (sin, exp, ln, …) applied to a matrix distribute element-wise. This works for a dense value and for a symbolic matrix alike — numeric cells fold, free-variable cells stay as f(cell) until bound:

// exp over a symbolic matrix: exp(x) stays symbolic, exp(0) folds to 1.0
Exp(_Matrix.ofRows(Vector(x, _Number(0.0)))).eval(env) match {
  case Left(m)  => m.toString
  case Right(v) => v.toString
}
// res12: String = "[[exp(x), 1.0]]"

Binding the free variables lets the whole matrix collapse to a dense value:

Exp(_Matrix.ofRows(Vector(x, _Number(0.0)))).eval(envX) match {
  case Right(v) => v.toString
  case Left(m)  => m.toString
}
// res13: String = "[[2.71828, 1.0]]"

Decompositions

Each decomposition returns its factors bundled as a 1×n row of matrices, so one indexing mechanism (at(result, 1, k), or tuple assignment in the REPL) serves them all:

Call	Result	Algorithm
lu(A)	[[L, U, P]], P·A = L·U	LU decomposition with partial pivoting
qr(A)	[[Q, R]], A = Q·R	QR decomposition by modified Gram–Schmidt
eigen(A)	[[λ₁, …, λₙ]]	QR algorithm with Wilkinson shifts; complex pairs come back as _Complex
eig(A)	[[V, D]], A·V = V·D	Eigendecomposition: eigenvector columns and the diagonal eigenvalue matrix
jordan(A)	[[P, J]], A = P·J·P⁻¹	Jordan normal form for diagonalizable input; a defective matrix stays symbolic

Parser.parse("eigen([[2, 1], [1, 2]])").get.eval(env).toExpression.toString
// res14: String = "[[3.0, 1.0]]"

Parser.parse("lu([[4, 3], [6, 3]])").get.eval(env).toExpression.toString
// res15: String = "[[[[1.0, 0.0], [0.66667, 1.0]], [[6.0, 3.0], [0.0, 1.0]], [[0.0, 1.0], [1.0, 0.0]]]]"

Matrix exponential

expm(A) computes e^A — the matrix exponential, a different operation from the integer power A^n. It uses scaling and squaring with a degree-13 Padé approximant, which is why it also works for defective matrices — input the eigendecomposition route has no basis for:

// nilpotent: the series terminates, so the result is exact
Parser.parse("expm([[0, 1], [0, 0]])").get.eval(env).toExpression.toString
// res16: String = "[[1.0, 1.0], [0.0, 1.0]]"

// a rotation generator exponentiates to the rotation matrix
Parser.parse("expm([[0, -1], [1, 0]])").get.eval(env).toExpression.toString
// res17: String = "[[0.5403, -0.84147], [0.84147, 0.5403]]"