Numerical linear algebra

The following operations are supported in-place and allocation-free

• Matrix inversion in \(n^3\) flops. Recall that inversion by Cholesky costs \((5/3)n^3\) flops and needs to allocate a matrix of same size.

• Computation of determinant

• Check (strict) positive definiteness

import sweepystats as sw
import numpy as np

def random_symmetric_matrix(n, eigmin=float('-inf')):
    """
    Simulates a n*n symmetric matrix with minimum eigenvalue set to `eigmin`.
    The resulting matrix is stored in column-major format.
    """
    A = np.random.rand(n, n)
    A = 0.5 * (A + A.T)
    # force eigenvalues to be >= eigmin
    if eigmin > float('-inf'):
        eval, evec = np.linalg.eig(A)
        eval[np.where(eval < 0)[0]] = eigmin
        A = evec * np.diag(eval) * evec.T
    # convert to column major
    return np.array(A, order='F')

Matrix inverses

If A is full rank, A.sweep() converts A.A into -inv(A.A):

# simulated data
data = random_symmetric_matrix(100)

# true inverse
Ainv = np.linalg.inv(data)

# sweep
A = sw.SweepMatrix(data)
A.sweep()

# check answer is correct
np.allclose(A.A, -Ainv)

100%|███████████████████████████████████████| 100/100 [00:00<00:00, 13296.26it/s]

True

Determinants

Determinants are computed by A.det(). The original matrix is untouched by default.

# simulate data
data = random_symmetric_matrix(100)

# true determinant
Adet = np.linalg.det(data)

# sweep
A = sw.SweepMatrix(data)
det = A.det()

# check answer is correct
np.allclose(det, Adet)

100%|██████████████████████████████████████| 100/100 [00:00<00:00, 151255.10it/s]
100%|██████████████████████████████████████| 100/100 [00:00<00:00, 144931.03it/s]

True

Checking of positive definiteness

This is a PD matrix:

A = sw.SweepMatrix(random_symmetric_matrix(100, eigmin=0.00001)) # this is PD
A.isposdef(tol=1e-15, verbose=False)

True

Simulate a matrix with negative eigenvalues:

A = sw.SweepMatrix(random_symmetric_matrix(100, eigmin=-10)) # should return false
A.isposdef(verbose=False)

False

Simulate a matrix with 0 as eigenvalue:

A = sw.SweepMatrix(random_symmetric_matrix(100, eigmin=0.0)) # this is PSD, should return false
A.isposdef(verbose=False)

False

Matrix rank:

Less well known is that the sweep operator can also be used to compute the rank of a matrix, by counting how many diagonals are non-zero prior to sweeping.

A = sw.SweepMatrix(random_symmetric_matrix(123, eigmin=0.0))
A.rank()

100%|███████████████████████████████████████| 123/123 [00:00<00:00, 57836.25it/s]
100%|███████████████████████████████████████| 123/123 [00:00<00:00, 56351.65it/s]

62

Check answer:

np.linalg.matrix_rank(A.A)

np.int64(62)