Skip to main content

Linear algebra

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

Solving and factoring

rm.linalg solves and decomposes real or complex matrices. solve and pinv give least-norm least-squares solutions for singular or non-square matrices. For repeated square solves, reuse lu factors. Their methods also provide an inverse, determinant and the determinant’s sign and log magnitude:

a = rm.array([[1., 2], [2, 5]])
f = rm.linalg.lu(a)
f.solve([[1., 2], [2, 1]]), f.inverse(), f.determinant(), f.log_determinant()
(Array([[1.0, 8.0], [-0.0, -3.0]], dtype=float64),
 Array([[5.0, -2.0], [-2.0, 1.0]], dtype=float64),
 1.0,
 (1.0, 0.0))

The separate triangular factors satisfy l @ u == a[p], with p recording the row permutation.

assert rm.norm(f.l @ f.u - a[f.p]) < 1e-12
f.l, f.u, f.p
(Array([[1.0, 0.0], [0.5, 1.0]], dtype=float64),
 Array([[2.0, 5.0], [0.0, -0.5]], dtype=float64),
 Array([1, 0], dtype=int64))

For a Hermitian positive definite matrix, cholesky reads the lower triangle and returns L in A = L L-adjoint.

lower = rm.linalg.cholesky(a)
assert rm.norm(lower @ rm.conj(lower).T - a) < 1e-12
lower
Array([[1.0, 0.0], [2.0, 1.0]], dtype=float64)

Fitting

Fit a line to four observations with lstsq. Its two results are the least-norm solution and singular values, rather than NumPy’s four-item result. Compute residuals with @:

x = rm.array([0., 1, 2, 3])
y = rm.array([1., 2, 2, 4])
design = rm.stack([rm.ones(4), x], axis=1)
coefficients, singular_values = rm.linalg.lstsq(design, y)
coefficients, y - design @ coefficients, singular_values
(Array([0.9, 0.9], dtype=float64),
 Array([0.1, 0.2, -0.7, 0.4], dtype=float64),
 Array([4.1000304, 1.0907568], dtype=float64))

Thin QR gives a Q with as many columns as the smaller matrix dimension. complete=True makes Q square. Both reconstruct the design matrix. The diagonal of R is real, with no sign convention.

q, r = rm.linalg.qr(design)
q_full, r_full = rm.linalg.qr(design, complete=True)
assert rm.norm(q @ r - design) < 1e-12
assert rm.norm(q_full @ r_full - design) < 1e-12
q.shape, r.shape, q_full.shape, r_full.shape
((4, 2), (2, 2), (4, 4), (4, 2))

For dependent columns, matrix_rank counts the retained singular values. solve returns the least-norm solution. pinv provides a reusable matrix for the same calculation. See exact arithmetic for rational solves, row reduction and the permanent.

singular = rm.array([[1., 2], [2, 4]])
pseudoinverse = rm.linalg.pinv(singular)
rm.linalg.matrix_rank(singular), rm.linalg.solve(singular, [1., 2]), pseudoinverse @ [1., 2]
(1, Array([0.2, 0.4], dtype=float64), Array([0.2, 0.4], dtype=float64))

Singular values and eigenvalues

svd returns Svd(u, s, v) for A = U diag(s) V-adjoint. The last factor is V, not NumPy’s V-adjoint. Reuse the factors for solves or a pseudoinverse. cond gives the ratio of the largest singular value to the smallest, the condition number in the 2-norm. cond and norm also take the other matrix norms that NumPy names, 1, inf and "fro":

factors = rm.linalg.svd(a)
factors.s, factors.solve([1., 2]), rm.linalg.cond(a), rm.linalg.cond(a, 1)
(Array([5.8284271, 0.17157288], dtype=float64),
 Array([1.0, 7.7715612e-16], dtype=float64),
 33.97056274847704,
 49.0)

Use rank=k for an approximate low-rank representation. The randomized truncated algorithm computes the requested leading components. For these four measurements of three features, one component captures most of the variation. Scale the rows of V-adjoint by s[:, None] when reconstructing:

measurements = rm.array([[1., 2, 3], [2, 4, 5], [3, 6, 8], [4, 8, 11]])
reduced = rm.linalg.svd(measurements, rank=1)
approximation = reduced.u @ (reduced.s[:, None] * rm.conj(reduced.v).T)
relative_error = rm.norm(measurements - approximation) / rm.norm(measurements)
reduced.u.shape, reduced.s.shape, reduced.v.shape, relative_error
((4, 1), (1,), (3, 1), 0.018182420257431917)

eigh reads the lower triangle of a symmetric or Hermitian matrix and returns ascending real eigenvalues. Eigenvectors are columns. General eig returns complex eigenvalues and vectors, including for this real rotation matrix. Vector signs and phases are not canonicalized. Use eigvalsh, eigvals or svdvals when you need values without vectors.

rotation = rm.array([[0., -1], [1, 0]])
eigen = rm.linalg.eig(rotation)
residual = rotation @ eigen.vectors - eigen.vectors * eigen.values[None, :]
assert rm.norm(residual) < 1e-12
rm.linalg.eigh(a).values, eigen.values
(Array([0.17157288, 5.8284271], dtype=float64),
 Array([1j, -1j], dtype=complex128))

Matrix products

rm.matmul and the @ operator multiply matrices. Integer matrices multiply in int64 and wrap at its range, as BPL’s +.× does. Matrix functions take stacks in the last two axes. The preceding frame axes broadcast from the front.

Here three systems share one matrix of right-hand sides. Each column is a separate right-hand side:

matrices = rm.stack([a, 2*a, 3*a])
rhs = rm.array([[1., 2], [2, 1]])
solutions = rm.linalg.solve(matrices, rhs)
assert rm.norm(matrices @ solutions - rhs[None, :, :]) < 1e-12
solutions
Array([[[1.0, 8.0], [-0.0, -3.0]], [[0.5, 4.0], [-0.0, -1.5]], [[0.33333333, 2.6666667], [-0.0, -1.0]]], dtype=float64)

Only a one-dimensional right-hand side is interpreted as a vector. For a batch of vectors, add a singleton column axis, solve, then select that column:

rhs_vectors = rm.array([[1., 2], [2, 1], [1, 1]])
vector_solutions = rm.linalg.solve(matrices, rhs_vectors[..., None])
assert rm.norm(matrices @ vector_solutions - rhs_vectors[..., None]) < 1e-12
vector_solutions[..., 0]
Array([[1.0, -0.0], [4.0, -1.5], [1.0, -0.33333333]], dtype=float64)

Products and solves pair contracted axes by key when both have keys. Different key sets raise ValueError. The weights below are ordered y, x, while the data columns are x, y. The result retains its row keys.

features = rm.array({'Ada': {'x': 1., 'y': 2}, 'Bea': {'x': 3., 'y': 4}})
weights = rm.array({'y': 10., 'x': 100.})
weighted = features @ weights
assert weighted.tolist() == [120., 340.]
weighted, weighted.keys
(Array([120.0, 340.0], dtype=float64, keys=(('Ada', 'Bea'),)),
 (('Ada', 'Bea'),))

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 matches or mismatches. Comparisons are exact unless you pass a comparison from the Mathish namespace. 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