sporco.linalg

Linear algebra functions.

Functions

inner(x, y[, axis])	Inner product of x and y on specified axis.
dot(a, b[, axis])	Matrix product of a and the specified axes of b.
valid_adjoint(A, AT, Ashape, ATshape[, eps])	Validate a transform and adjoint transform pair.
solvedbi_sm(ah, rho, b[, c, axis])	Solve a diagonal block linear system with a scaled identity term using the Sherman-Morrison equation.
solvedbi_sm_c(ah, a, rho[, axis])	Compute cached component used by solvedbi_sm.
solvedbd_sm(ah, d, b[, c, axis])	Solve a diagonal block linear system with a diagonal term using the Sherman-Morrison equation.
solvedbd_sm_c(ah, a, d[, axis])	Compute cached component used by solvedbd_sm.
solvemdbi_ism(ah, rho, b, axisM, axisK)	Solve a multiple diagonal block linear system with a scaled identity term by iterated application of the Sherman-Morrison equation.
solvemdbi_rsm(ah, rho, b, axisK[, dimN])	Solve a multiple diagonal block linear system with a scaled identity term by repeated application of the Sherman-Morrison equation.
solvemdbi_cg(ah, rho, b, axisM, axisK[, ...])	Solve a multiple diagonal block linear system with a scaled identity term using CG.
lu_factor(A, rho[, check_finite])	Compute LU factorisation of either \(A^T A + \rho I\) or \(A A^T + \rho I\), depending on which matrix is smaller.
lu_solve_ATAI(A, rho, b, lu, piv[, check_finite])	Solve the linear system \((A^T A + \rho I)\mathbf{x} = \mathbf{b}\) or \((A^T A + \rho I)X = B\) using scipy.linalg.lu_solve.
lu_solve_AATI(A, rho, b, lu, piv[, check_finite])	Solve the linear system \((A A^T + \rho I)\mathbf{x} = \mathbf{b}\) or \((A A^T + \rho I)X = B\) using scipy.linalg.lu_solve.
cho_factor(A, rho[, lower, check_finite])	Compute Cholesky factorisation of either \(A^T A + \rho I\) or \(A A^T + \rho I\), depending on which matrix is smaller.
cho_solve_ATAI(A, rho, b, c[, lwr, check_finite])	Solve the linear system \((A^T A + \rho I)\mathbf{x} = \mathbf{b}\) or \((A^T A + \rho I)X = B\) using scipy.linalg.cho_solve.
cho_solve_AATI(A, rho, b, c[, lwr, check_finite])	Solve the linear system \((A A^T + \rho I)\mathbf{x} = \mathbf{b}\) or \((A A^T + \rho I)X = B\) using scipy.linalg.cho_solve.
solve_symmetric_sylvester(A, B, C, alpha)	Solve a variant of the Sylvester equation with symmetric matrices.
block_circulant(A)	Construct a block circulant matrix from a tuple of arrays.
rrs(ax, b)	Relative residual of the solution to a linear equation.
pca(U[, centre])	Compute the PCA basis for columns of input array U.
nkp(A, bshape, cshape)	Solve the Nearest Kronecker Product problem.
kpsvd(A, bshape, cshape)	Compute the Kronecker Product SVD.
proj_l2ball(b, s, r[, axes])	Projection onto the \(\ell_2\) ball.

Function Descriptions

sporco.linalg.inner(x, y, axis=- 1)

Inner product of x and y on specified axis.

Compute inner product of x and y on specified axis, equivalent to np.sum(x * y, axis=axis, keepdims=True).

Parameters

xarray_like

Input array x

yarray_like

Input array y

axisint, optional (default -1)

Axis over which to compute the sum

Returns

yndarray

Inner product array equivalent to summing x*y over the specified axis

sporco.linalg.dot(a, b, axis=- 2)

Matrix product of a and the specified axes of b.

Compute the matrix product of a and the specified axes of b, with broadcasting over the remaining axes of b. This function is a generalisation of numpy.dot, supporting sum product over an arbitrary axis instead of just over the last axis.

If a and b are both 2D arrays, dot gives the same result as numpy.dot. If b has more than 2 axes, the result is obtained as follows (where a has shape (M0, M1) and b has shape (N0, N1, ..., M1, Nn, ...)):

1. Reshape a to shape ( 1,  1, ..., M0, M1,  1, ...)

2. Reshape b to shape (N0, N1, ...,  1, M1, Nn, ...)

3. Take the broadcast product and sum over the specified axis (the axis with dimension M1 in this example) to give an array of shape (N0, N1, ...,  M0,  1, Nn, ...)

4. Remove the singleton axis created by the summation to give an array of shape (N0, N1, ...,  M0, Nn, ...)

Parameters

aarray_like, 2D

First component of product

barray_like, 2D or greater

Second component of product

axisinteger, optional (default -2)

Axis of b over which sum is to be taken

Returns

prodndarray

Matrix product of a and specified axes of b, with broadcasting over the remaining axes of b

sporco.linalg.valid_adjoint(A, AT, Ashape, ATshape, eps=1e-07)

Validate a transform and adjoint transform pair.

Check whether transform AT is the adjoint of A. The test exploits the identity

\[\mathbf{y}^T (A \mathbf{x}) = (\mathbf{y}^T A) \mathbf{x} = (A^T \mathbf{y})^T \mathbf{x}\]

by computing \(\mathbf{u} = A \mathbf{x}\) and \(\mathbf{v} = A^T \mathbf{y}\) for random \(\mathbf{x}\) and \(\mathbf{y}\) and confirming that \(\| \mathbf{y}^T \mathbf{u} - \mathbf{v}^T \mathbf{x} \|_2 < \epsilon\) since

\[\mathbf{y}^T \mathbf{u} = \mathbf{y}^T (A \mathbf{x}) = (A^T \mathbf{y})^T \mathbf{x} = \mathbf{v}^T \mathbf{x}\]

when \(A^T\) is a valid adjoint of \(A\).

Parameters

Afunction

Primary function

ATfunction

Adjoint function

Ashapetuple

Shape of input array expected by function A

ATshapetuple

Shape of input array expected by function AT

epsfloat or None, optional (default 1e-7)

Error threshold for validation of AT as adjoint of A. If None, the relative error is returned instead of a boolean value.

Returns

errboolean or float

Boolean value indicating that validation passed, or relative error of test, depending on type of parameter eps

sporco.linalg.solvedbi_sm(ah, rho, b, c=None, axis=4)

Solve a diagonal block linear system with a scaled identity term using the Sherman-Morrison equation.

The solution is obtained by independently solving a set of linear systems of the form (see [53])

\[(\rho I + \mathbf{a} \mathbf{a}^H ) \; \mathbf{x} = \mathbf{b} \;\;.\]

In this equation inner products and matrix products are taken along the specified axis of the corresponding multi-dimensional arrays; the solutions are independent over the other axes.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

rhofloat

Linear system parameter \(\rho\)

barray_like

Linear system component \(\mathbf{b}\)

carray_like, optional (default None)

Solution component \(\mathbf{c}\) that may be pre-computed using solvedbi_sm_c and cached for re-use.

axisint, optional (default 4)

Axis along which to solve the linear system

Returns

xndarray

Linear system solution \(\mathbf{x}\)

sporco.linalg.solvedbi_sm_c(ah, a, rho, axis=4)

Compute cached component used by solvedbi_sm.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

aarray_like

Linear system component \(\mathbf{a}\)

rhofloat

Linear system parameter \(\rho\)

axisint, optional (default 4)

Axis along which to solve the linear system

Returns

cndarray

Argument \(\mathbf{c}\) used by solvedbi_sm

sporco.linalg.solvedbd_sm(ah, d, b, c=None, axis=4)

Solve a diagonal block linear system with a diagonal term using the Sherman-Morrison equation.

The solution is obtained by independently solving a set of linear systems of the form (see [53])

\[(\mathbf{d} + \mathbf{a} \mathbf{a}^H ) \; \mathbf{x} = \mathbf{b} \;\;.\]

In this equation inner products and matrix products are taken along the specified axis of the corresponding multi-dimensional arrays; the solutions are independent over the other axes.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

darray_like

Linear system parameter \(\mathbf{d}\)

barray_like

Linear system component \(\mathbf{b}\)

carray_like, optional (default None)

Solution component \(\mathbf{c}\) that may be pre-computed using solvedbd_sm_c and cached for re-use.

axisint, optional (default 4)

Axis along which to solve the linear system

Returns

xndarray

Linear system solution \(\mathbf{x}\)

sporco.linalg.solvedbd_sm_c(ah, a, d, axis=4)

Compute cached component used by solvedbd_sm.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

aarray_like

Linear system component \(\mathbf{a}\)

darray_like

Linear system parameter \(\mathbf{d}\)

axisint, optional (default 4)

Axis along which to solve the linear system

Returns

cndarray

Argument \(\mathbf{c}\) used by solvedbd_sm

sporco.linalg.solvemdbi_ism(ah, rho, b, axisM, axisK)

Solve a multiple diagonal block linear system with a scaled identity term by iterated application of the Sherman-Morrison equation.

The computation is performed in a way that avoids explictly constructing the inverse operator, leading to an \(O(K^2)\) time cost.

The solution is obtained by independently solving a set of linear systems of the form (see [53])

\[(\rho I + \mathbf{a}_0 \mathbf{a}_0^H + \mathbf{a}_1 \mathbf{a}_1^H + \; \ldots \; + \mathbf{a}_{K-1} \mathbf{a}_{K-1}^H) \; \mathbf{x} = \mathbf{b}\]

where each \(\mathbf{a}_k\) is an \(M\)-vector. The sums, inner products, and matrix products in this equation are taken along the \(M\) and \(K\) axes of the corresponding multi-dimensional arrays; the solutions are independent over the other axes.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

rhofloat

Linear system parameter \(\rho\)

barray_like

Linear system component \(\mathbf{b}\)

axisMint

Axis in input corresponding to index m in linear system

axisKint

Axis in input corresponding to index k in linear system

Returns

xndarray

Linear system solution \(\mathbf{x}\)

sporco.linalg.solvemdbi_rsm(ah, rho, b, axisK, dimN=2)

Solve a multiple diagonal block linear system with a scaled identity term by repeated application of the Sherman-Morrison equation.

The computation is performed by explictly constructing the inverse operator, leading to an \(O(K)\) time cost and \(O(M^2)\) memory cost, where \(M\) is the dimension of the axis over which \(M\) inner products are taken.

The solution is obtained by independently solving a set of linear systems of the form (see [53])

\[(\rho I + \mathbf{a}_0 \mathbf{a}_0^H + \mathbf{a}_1 \mathbf{a}_1^H + \; \ldots \; + \mathbf{a}_{K-1} \mathbf{a}_{K-1}^H) \; \mathbf{x} = \mathbf{b}\]

where each \(\mathbf{a}_k\) is an \(M\)-vector. The sums, inner products, and matrix products in this equation are taken along the \(M\) and \(K\) axes of the corresponding multi-dimensional arrays; the solutions are independent over the other axes.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

rhofloat

Linear system parameter \(\rho\)

barray_like

Linear system component \(\mathbf{b}\)

axisKint

Axis in input corresponding to index k in linear system

dimNint, optional (default 2)

Number of spatial dimensions arranged as leading axes in input array. Axis M is taken to be at dimN+2.

Returns

xndarray

Linear system solution \(\mathbf{x}\)

sporco.linalg.solvemdbi_cg(ah, rho, b, axisM, axisK, tol=1e-05, mit=1000, isn=None)

Solve a multiple diagonal block linear system with a scaled identity term using CG.

Solve a multiple diagonal block linear system with a scaled identity term using Conjugate Gradient (CG) via scipy.sparse.linalg.cg.

The solution is obtained by independently solving a set of linear systems of the form (see [53])

\[(\rho I + \mathbf{a}_0 \mathbf{a}_0^H + \mathbf{a}_1 \mathbf{a}_1^H + \; \ldots \; + \mathbf{a}_{K-1} \mathbf{a}_{K-1}^H) \; \mathbf{x} = \mathbf{b}\]

where each \(\mathbf{a}_k\) is an \(M\)-vector. The inner products and matrix products in this equation are taken along the \(M\) and \(K\) axes of the corresponding multi-dimensional arrays; the solutions are independent over the other axes.

Warning: This function is not supported under Windows due to an access violation error of unknown origin.

Parameters

aharray_like

Linear system component \(\mathbf{a}^H\)

rhofloat

Parameter rho

barray_like

Linear system component \(\mathbf{b}\)

axisMint

Axis in input corresponding to index m in linear system

axisKint

Axis in input corresponding to index k in linear system

tolfloat

CG tolerance

mitint

CG maximum iterations

isnarray_like

CG initial solution

Returns

xndarray

Linear system solution \(\mathbf{x}\)

cgitint

Number of CG iterations

sporco.linalg.lu_factor(A, rho, check_finite=True)

Compute LU factorisation of either \(A^T A + \rho I\) or \(A A^T + \rho I\), depending on which matrix is smaller.

Parameters

Aarray_like

Array \(A\)

rhofloat

Scalar \(\rho\)

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

lundarray

Matrix containing U in its upper triangle, and L in its lower triangle, as returned by scipy.linalg.lu_factor

pivndarray

Pivot indices representing the permutation matrix P, as returned by scipy.linalg.lu_factor

sporco.linalg.lu_solve_ATAI(A, rho, b, lu, piv, check_finite=True)

Solve the linear system \((A^T A + \rho I)\mathbf{x} = \mathbf{b}\) or \((A^T A + \rho I)X = B\) using scipy.linalg.lu_solve.

Parameters

Aarray_like

Matrix \(A\)

rhofloat

Scalar \(\rho\)

barray_like

Vector \(\mathbf{b}\) or matrix \(B\)

luarray_like

Matrix containing U in its upper triangle, and L in its lower triangle, as returned by scipy.linalg.lu_factor

pivarray_like

Pivot indices representing the permutation matrix P, as returned by scipy.linalg.lu_factor

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

xndarray

Solution to the linear system

sporco.linalg.lu_solve_AATI(A, rho, b, lu, piv, check_finite=True)

Solve the linear system \((A A^T + \rho I)\mathbf{x} = \mathbf{b}\) or \((A A^T + \rho I)X = B\) using scipy.linalg.lu_solve.

Parameters

Aarray_like

Matrix \(A\)

rhofloat

Scalar \(\rho\)

barray_like

Vector \(\mathbf{b}\) or matrix \(B\)

luarray_like

Matrix containing U in its upper triangle, and L in its lower triangle, as returned by scipy.linalg.lu_factor

pivarray_like

Pivot indices representing the permutation matrix P, as returned by scipy.linalg.lu_factor

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

xndarray

Solution to the linear system

sporco.linalg.cho_factor(A, rho, lower=False, check_finite=True)

Compute Cholesky factorisation of either \(A^T A + \rho I\) or \(A A^T + \rho I\), depending on which matrix is smaller.

Parameters

Aarray_like

Array \(A\)

rhofloat

Scalar \(\rho\)

lowerbool, optional (default False)

Flag indicating whether lower or upper triangular factors are computed

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

cndarray

Matrix containing lower or upper triangular Cholesky factor, as returned by scipy.linalg.cho_factor

lwrbool

Flag indicating whether the factor is lower or upper triangular

sporco.linalg.cho_solve_ATAI(A, rho, b, c, lwr=False, check_finite=True)

Solve the linear system \((A^T A + \rho I)\mathbf{x} = \mathbf{b}\) or \((A^T A + \rho I)X = B\) using scipy.linalg.cho_solve.

Parameters

Aarray_like

Matrix \(A\)

rhofloat

Scalar \(\rho\)

barray_like

Vector \(\mathbf{b}\) or matrix \(B\)

carray_like

Matrix containing lower or upper triangular Cholesky factor, as returned by scipy.linalg.cho_factor

lwrbool, optional (default False)

Flag indicating whether lower or upper triangular factors are computed

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

xndarray

Solution to the linear system

sporco.linalg.cho_solve_AATI(A, rho, b, c, lwr=False, check_finite=True)

Solve the linear system \((A A^T + \rho I)\mathbf{x} = \mathbf{b}\) or \((A A^T + \rho I)X = B\) using scipy.linalg.cho_solve.

Parameters

Aarray_like

Matrix \(A\)

rhofloat

Scalar \(\rho\)

barray_like

Vector \(\mathbf{b}\) or matrix \(B\)

carray_like

Matrix containing lower or upper triangular Cholesky factor, as returned by scipy.linalg.cho_factor

lwrbool, optional (default False)

Flag indicating whether lower or upper triangular factors are computed

check_finitebool, optional (default False)

Flag indicating whether the input array should be checked for Inf and NaN values

Returns

xndarray

Solution to the linear system

sporco.linalg.solve_symmetric_sylvester(A, B, C, alpha)

Solve a variant of the Sylvester equation with symmetric matrices.

Solve for \(X\) in the equation \(A X B + \alpha X = C\), where \(A\) and \(B\) are diagonal or real symmetric matrices.

Parameters

Aarray_like or tuple

Matrix \(A\). If \(A\) is diagonal, this should be a 1d array or a \(N \times 1\) 2d array representing the diagonal. If \(A\) is symmetric, this should be a 2d array or a tuple (LambdaA, QA) representing the eigenvalue decomposition of \(A\) such that \(Q_A \Lambda_A Q_A^T = A\)

Barray_like or tuple

Matrix \(B\). If \(B\) is diagonal, this should be a 1d array or a \(1 \times M\) 2d array representing the diagonal. If \(B\) is symmetric, this should be a 2d array or a tuple (LambdaB, QB) representing the eigenvalue decomposition of \(B\) such that \(Q_B \Lambda_B Q_B^T = B\)

Carray_like

Matrix \(C\) as a 2d array

alphafloat

Scalar \(\alpha\)

Returns

Xndarray

Solution to the linear system

sporco.linalg.block_circulant(A)

Construct a block circulant matrix from a tuple of arrays.

Construct a block circulant matrix from a tuple of arrays. This is a block-matrix variant of scipy.linalg.circulant.

Parameters

Atuple of array_like

Tuple of arrays corresponding to the first block column of the output block matrix

Returns

Bndarray

Output array

sporco.linalg.rrs(ax, b)

Relative residual of the solution to a linear equation.

The standard relative residual for the linear system \(A \mathbf{x} = \mathbf{b}\) is \(\|\mathbf{b} - A \mathbf{x}\|_2 / \|\mathbf{b}\|_2\). This function computes a variant \(\|\mathbf{b} - A \mathbf{x}\|_2 / \max(\|A\mathbf{x}\|_2, \|\mathbf{b}\|_2)\) that is robust to the case \(\mathbf{b} = 0\).

Parameters

axarray_like

Linear component \(A \mathbf{x}\) of equation

barray_like

Constant component \(\mathbf{b}\) of equation

Returns

xfloat

Relative residual

sporco.linalg.pca(U, centre=False)

Compute the PCA basis for columns of input array U.

Parameters

Uarray_like

2D data array with rows corresponding to different variables and columns corresponding to different observations

centerbool, optional (default False)

Flag indicating whether to centre data

Returns

Bndarray

A 2D array representing the PCA basis; each column is a PCA component. B.T is the analysis transform into the PCA representation, and B is the corresponding synthesis transform

Sndarray

The eigenvalues of the PCA components

Cndarray or None

None if centering is disabled, otherwise the mean of the data matrix subtracted in performing the centering

sporco.linalg.nkp(A, bshape, cshape)

Solve the Nearest Kronecker Product problem.

Given matrix \(A\), find matrices \(B\) and \(C\), of the specified sizes, such that \(B\) and \(C\) solve the problem [47]

\[\mathrm{argmin}_{B, C} \| A - B \otimes C \| \;.\]

Parameters

Aarray_like

2D input array

bshapetuple (Mb, Nb)

The desired shape of returned array \(B\)

cshapetuple (Mc, Nc)

The desired shape of returned array \(C\)

Returns

Bndarray

2D output array \(B\)

Cndarray

2D output array \(C\)

sporco.linalg.kpsvd(A, bshape, cshape)

Compute the Kronecker Product SVD.

Given matrix \(A\), find matrices \(B_i\) and \(C_i\), of the specified sizes, such that

\[A = \sum_i \sigma_i B_i \otimes C_i\]

and \(\sum_i^n \sigma_i B_i \otimes C_i\) is the best \(n\) term approximation to \(A\) [47].

Parameters

Aarray_like

2D input array

bshapetuple (Mb, Nb)

The desired shape of arrays \(B_i\)

cshapetuple (Mc, Nc)

The desired shape of arrays \(C_i\)

Returns

Sndarray

1D array of \(\sigma_i\) values

Bndarray

3D array of \(B_i\) matrices with index \(i\) on the last axis

Cndarray

3D array of \(C_i\) matrices with index \(i\) on the last axis

sporco.linalg.proj_l2ball(b, s, r, axes=None)

Projection onto the \(\ell_2\) ball.

Project \(\mathbf{b}\) onto the \(\ell_2\) ball of radius \(r\) about \(\mathbf{s}\), i.e. \(\{ \mathbf{x} : \|\mathbf{x} - \mathbf{s} \|_2 \leq r \}\). Note that proj_l2ball(b, s, r) is equivalent to prox.proj_l2 (b - s, r) + s.

NB: This function is to be deprecated; please use prox.proj_l2 instead (see note above about interface differences).

Parameters

barray_like

Vector \(\mathbf{b}\) to be projected

sarray_like

Centre of \(\ell_2\) ball \(\mathbf{s}\)

rfloat

Radius of ball

axessequence of ints, optional (default all axes)

Axes over which to compute \(\ell_2\) norms

Returns

xndarray

Projection of \(\mathbf{b}\) into ball