Linear algebra

Matrix products, reusable factorizations, eigenvalues and least-squares solutions.
import rustymath as rm

Solving and factoring

rm.linalg provides real and complex matrix solves and decompositions, including solve, det, lu, cholesky, qr, svd, eigh, eig, lstsq, pinv and matrix_rank. They take stacks of matrices, whose leading axes broadcast from the front as other functions’ arguments do. solve and pinv give least-norm least-squares solutions for singular or non-square matrices. permanent takes square matrices. echelon reduces one 2-D matrix. lu gives the packed factors and the row positions. Its .l and .u give the separate triangular matrices. Reuse the factors for solves, inverses and determinants:

a = rm.array([[4., 1], [1, 3]])
f = rm.linalg.lu(a)
f.solve([1., 2]), f.determinant(), rm.linalg.eigh(a).values, rm.linalg.solve([[1., 2], [2, 4]], [1., 2])
(Array([0.090909091, 0.63636364], dtype=float64),
 11.0,
 Array([2.381966, 4.618034], dtype=float64),
 Array([0.2, 0.4], dtype=float64))

Singular values and eigenvalues

svd returns Svd(u, s, v), where the reconstruction is U diag(s) V-adjoint. The last factor is V, not NumPy’s V-adjoint. cond gives the 2-norm condition number, the largest singular value divided by the smallest. Reuse these factors for solves or a pseudoinverse. eigh reads the lower triangle of a real symmetric or complex Hermitian matrix. eig takes a general square matrix. Vector signs and phases are not canonicalized:

factors = rm.linalg.svd(a)
factors.s, factors.solve([1., 2]), rm.linalg.cond(a), rm.linalg.eigh(a).values
(Array([4.618034, 2.381966], dtype=float64),
 Array([0.090909091, 0.63636364], dtype=float64),
 1.9387489019317514,
 Array([2.381966, 4.618034], dtype=float64))

Matrix products

rm.matmul and the @ operator multiply matrices. Integer matrices multiply in int64. A product that doesn’t fit raises OverflowError. Products and solves pair keyed axes by key. Their results keep the keys of the axes they index.

rm.linalg.inner(f, g) applies a native binary function g along the contracted axis, then folds with f from the right, as APL’s inner product f.g does. Recognized pairs use specialized kernels; other pairs use bounded blocks of terms. With add and equal or not_equal, it counts exact matches or mismatches; comparisons in other pairs use the default tolerance. Squaring a matrix of edge lengths with minimum and add gives the shortest paths of at most two edges:

d = rm.array([[0, 3, float('inf')], [float('inf'), 0, 1], [2, float('inf'), 0]])
minplus = rm.linalg.inner(rm.minimum, rm.add)
minplus(d, d)
Array([[0.0, 3.0, 4.0], [3.0, 0.0, 1.0], [2.0, 5.0, 0.0]], dtype=float64)

f need not be associative. With subtraction, the four products fold as 1 - (2 - (3 - 4)), not as a left-to-right subtraction. A contraction with one item returns that product unchanged; an empty contraction needs a known identity for f. The functions keep their usual dtype restrictions.

rm.linalg.inner(rm.subtract, rm.multiply)([1., 2, 3, 4], [1., 1, 1, 1])
-2.0