Linear Algebra

Linear algebra on matrices kept in their cheapest form.

A Operator holds its matrix in one of the forms here, and the solvers compute with them: a DenseOperator (a NumPy array), a SparseOperator (a SciPy sparse matrix), a DiagonalOperator (the diagonal only), a BlockDiagonalOperator (a form on each diagonal block, zeros between them) or a KroneckerOperator (Kronecker factors). None of them has labels: the labels are the Space on each axis of an Operator.

Each is a scipy.sparse.linalg.LinearOperator, so matvec, rmatvec, matmat and .T work as in SciPy. On top of that:

• A @ B of two forms keeps the cheapest form that is exact: diagonal times diagonal is diagonal; sparse or diagonal times sparse is sparse; a product with a dense side is dense; a Kronecker form times a sparse one (H @ S_0) is computed through the factors and is dense. A @ x with an array is an array, as in SciPy.

• A + B and A - B follow the same rule (an array counts as dense), and a scalar * or / keeps the form.

• A product of a sparse form and a dense array (H @ X, X @ H.T) runs on threads when it is large, one block of the sparse matrix’s rows to a thread, and is the same bit for bit as on one. The threads are the CPUs the process may run on (a Slurm job’s, not the node’s), at most OMP_NUM_THREADS when that is set.

• dense(), diagonal(), trace(), and for a symmetric positive-definite form solve(), logdet() and sqrt(), each done the way its form makes cheap.

• np.asarray(A) is dense(), so NumPy functions take a form as its matrix; np.array(A) is a copy, as for an array.

A form owns its data, as a Vector does: it copies the array or matrix it is given and makes the copy read-only, so changing the original afterwards leaves the form as it was, and writing into the form’s own array raises ValueError. The forms fips builds from one another (a product, a transpose, a selection) hold the arrays they computed or the read-only ones of their parents, never a second copy.

The forms are real (float64) and pickle with the standard protocol.

Forms

LinearOperator(shape)	A real matrix in one of fips's forms, for multiplying, adding, transposing and solving.
DenseOperator(array)	A matrix held as a dense 2D NumPy array.
SparseOperator(matrix)	A matrix held as a SciPy sparse matrix.
DiagonalOperator(values)	A square matrix that is zero off its main diagonal, held as the diagonal.
BlockDiagonalOperator(*args, **kwargs)	A square matrix made of forms on its diagonal blocks, zero between them.

Building

as_linear_operator(x)

Return x as the form that fits how it is stored.

fips.linalg.LinearOperatorLike

The type of what as_linear_operator() takes: a form, a SciPy sparse matrix or array, or a 2D array-like.

Linear algebra

LinearOperator.dot(x)	Return self @ x, or the form scaled by a scalar x.
LinearOperator.solve(B)	Solve self @ X = B for a symmetric positive-definite form.
LinearOperator.logdet()	Return the logarithm of the determinant, for a symmetric positive-definite form.
LinearOperator.sqrt()	Return a square root \(L\) with \(L L^T\) equal to this symmetric positive-definite form.
LinearOperator.dense()	Return the form as a dense NumPy array.
LinearOperator.diagonal()	Return the main diagonal; a sparse or diagonal form reads it from its entries.
LinearOperator.trace()	Return the sum of the main diagonal.

Each form’s page lists its own way of doing these, where it has one (a diagonal solves by division, a sparse form with a sparse LU).

Stored data

DenseOperator.array	The matrix, a read-only array.
SparseOperator.matrix	The matrix, as a SciPy sparse array whose arrays are read-only.
DiagonalOperator.values	The diagonal, a read-only array.
BlockDiagonalOperator.blocks	The diagonal blocks, each in its own form.

Testing

dottest(op[, rtol, seed])

Check that op.rmatvec is the adjoint of op.matvec.

Kronecker factors

Kronecker covariances: a covariance kept as the product of its marginal correlations.

A covariance over a grid of times and places is often the product of a correlation in time and one in space. KroneckerOperator keeps the two (or more) factors and never forms the n x n matrix: products, solves, the diagonal and the log-determinant go through the factors. It is one of the forms of fips.linalg, so a Covariance holds one as it holds a dense or a diagonal matrix.

KroneckerOperator(factors[, std, positions])

A covariance kept as the Kronecker product of its marginal correlations.

Linear algebra

KroneckerOperator.solve(B)	Solve self @ X = B through the factors.
KroneckerOperator.logdet()	Return the log-determinant through the factors.
KroneckerOperator.dense()	Return the operator as a dense NumPy array.
KroneckerOperator.diagonal()	Return the main diagonal, std**2 * prod_k diag(C_k)[pos_k].

Stored data

KroneckerOperator.factors	The marginal correlation matrices, one per marginal, read-only.
KroneckerOperator.positions	For each factor, the position of each row in it (pos_k), read-only.
KroneckerOperator.std	The standard deviation of each row, read-only.