import rustymath as rmLinear algebra
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