sporco.linalg

Linear algebra functions

Functions

complex_dtype(dtype) Construct the corresponding complex dtype for a given real dtype, e.g.
pyfftw_byte_aligned(array[, dtype, n]) Construct a byte-aligned array for efficient use by pyfftw.
pyfftw_empty_aligned(shape, dtype[, order, n]) Construct an empty byte-aligned array for efficient use by pyfftw.
pyfftw_rfftn_empty_aligned(shape, axes, dtype) Construct an empty byte-aligned array for efficient use by pyfftw functions pyfftw.interfaces.numpy_fft.rfftn and pyfftw.interfaces.numpy_fft.irfftn.
fftn(a[, s, axes, norm]) Compute the multi-dimensional discrete Fourier transform.
ifftn(a[, s, axes, norm]) Compute the multi-dimensional inverse discrete Fourier transform.
rfftn(a[, s, axes, norm]) Compute the multi-dimensional discrete Fourier transform for real input.
irfftn(a[, s, axes, norm]) Compute the inverse of the multi-dimensional discrete Fourier transform for real input.
dctii(x[, axes]) Compute a multi-dimensional DCT-II over specified array axes.
idctii(x[, axes]) Compute a multi-dimensional inverse DCT-II over specified array axes.
fftconv(a, b[, axes]) Compute a multi-dimensional convolution via the Discrete Fourier Transform.
inner(x, y[, axis]) Compute inner product of x and y on specified axis, equivalent to np.sum(x * y, axis=axis, keepdims=True).
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 Conjugate Gradient (CG) via scipy.sparse.linalg.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.
zpad(x, pd, ax) Zero-pad array x with pd = (leading, trailing) zeros on axis ax.
Gax(x, ax) Compute gradient of x along axis ax.
GTax(x, ax) Compute transpose of gradient of x along axis ax.
GradientFilters(ndim, axes, axshp[, dtype]) Construct a set of filters for computing gradients in the frequency domain.
zdivide(x, y) Return x/y, with 0 instead of NaN where y is 0.
shrink1(x, alpha) Scalar shrinkage/soft thresholding function
shrink2(x, alpha[, axis]) Vector shrinkage/soft thresholding function
shrink12(x, alpha, beta[, axis]) Compound shrinkage/soft thresholding function [38] [10]
proj_l2ball(b, s, r[, axes]) Project \(\mathbf{b}\) into the \(\ell_2\) ball of radius \(r\) about \(\mathbf{s}\), i.e.
promote16(u[, fn]) Utility function for use with functions that do not support arrays of dtype np.float16.
atleast_nd(n, u) If the input array has fewer than n dimensions, append singleton dimensions so that it is n dimensional.
split(u[, axis]) Split an array into a list of arrays on the specified axis.
blockcirculant(A) Construct a block circulant matrix from a tuple of arrays.
fl2norm2(xf[, axis]) Compute the squared \(\ell_2\) norm in the DFT domain, taking into account the unnormalised DFT scaling, i.e.
rfl2norm2(xf, xs[, axis]) Compute the squared \(\ell_2\) norm in the DFT domain, taking into account the unnormalised DFT scaling, i.e.
rrs(ax, b) Compute relative residual \(\|\mathbf{b} - A \mathbf{x}\|_2 / \|\mathbf{b}\|_2\) of the solution to a linear equation \(A \mathbf{x} = \mathbf{b}\).

Function Descriptions

sporco.linalg.complex_dtype(dtype)[source]

Construct the corresponding complex dtype for a given real dtype, e.g. the complex dtype corresponding to np.float32 is np.complex64.

Parameters:
dtype : dtype

A real dtype, e.g. np.float32, np.float64

Returns:
cdtype : dtype

The complex dtype corresponding to the input dtype

sporco.linalg.pyfftw_byte_aligned(array, dtype=None, n=None)[source]

Construct a byte-aligned array for efficient use by pyfftw. This function is a wrapper for pyfftw.byte_align

Parameters:
array : ndarray

Input array

dtype : dtype, optional (default None)

Output array dtype

n : int, optional (default None)

Output array should be aligned to n-byte boundary

Returns:
a : ndarray

Array with required byte-alignment

sporco.linalg.pyfftw_empty_aligned(shape, dtype, order='C', n=None)

Construct an empty byte-aligned array for efficient use by pyfftw. This function is a wrapper for pyfftw.empty_aligned

Parameters:
shape : sequence of ints

Output array shape

dtype : dtype

Output array dtype

order : {‘C’, ‘F’}, optional (default ‘C’)

Specify whether arrays should be stored in row-major (C-style) or column-major (Fortran-style) order

n : int, optional (default None)

Output array should be aligned to n-byte boundary

Returns:
a : ndarray

Empty array with required byte-alignment

sporco.linalg.pyfftw_rfftn_empty_aligned(shape, axes, dtype, order='C', n=None)

Construct an empty byte-aligned array for efficient use by pyfftw functions pyfftw.interfaces.numpy_fft.rfftn and pyfftw.interfaces.numpy_fft.irfftn. The shape of the empty array is appropriate for the output of pyfftw.interfaces.numpy_fft.rfftn applied to an array of the shape specified by parameter shape, and for the input of the corresponding pyfftw.interfaces.numpy_fft.irfftn call that reverses this operation.

Parameters:
shape : sequence of ints

Output array shape

axes : sequence of ints

Axes on which the FFT will be computed

dtype : dtype

Real dtype from which the complex dtype of the output array is derived

order : {‘C’, ‘F’}, optional (default ‘C’)

Specify whether arrays should be stored in row-major (C-style) or column-major (Fortran-style) order

n : int, optional (default None)

Output array should be aligned to n-byte boundary

Returns:
a : ndarray

Empty array with required byte-alignment

sporco.linalg.fftn(a, s=None, axes=None, norm=None)[source]

Compute the multi-dimensional discrete Fourier transform. This function is a wrapper for pyfftw.interfaces.numpy_fft.fftn, with an interface similar to that of numpy.fft.fftn.

Parameters:
a : array_like

Input array (can be complex)

s : sequence of ints, optional (default None)

Shape of the output along each transformed axis (input is cropped or zero-padded to match).

axes : sequence of ints, optional (default None)

Axes over which to compute the DFT.

Returns:
af : complex ndarray

DFT of input array

sporco.linalg.ifftn(a, s=None, axes=None, norm=None)[source]

Compute the multi-dimensional inverse discrete Fourier transform. This function is a wrapper for pyfftw.interfaces.numpy_fft.ifftn, with an interface similar to that of numpy.fft.ifftn.

Parameters:
a : array_like

Input array (can be complex)

s : sequence of ints, optional (default None)

Shape of the output along each transformed axis (input is cropped or zero-padded to match).

axes : sequence of ints, optional (default None)

Axes over which to compute the inverse DFT.

Returns:
af : complex ndarray

Inverse DFT of input array

sporco.linalg.rfftn(a, s=None, axes=None, norm=None)[source]

Compute the multi-dimensional discrete Fourier transform for real input. This function is a wrapper for pyfftw.interfaces.numpy_fft.rfftn, with an interface similar to that of numpy.fft.rfftn.

Parameters:
a : array_like

Input array (taken to be real)

s : sequence of ints, optional (default None)

Shape of the output along each transformed axis (input is cropped or zero-padded to match).

axes : sequence of ints, optional (default None)

Axes over which to compute the DFT.

Returns:
af : complex ndarray

DFT of input array

sporco.linalg.irfftn(a, s=None, axes=None, norm=None)[source]

Compute the inverse of the multi-dimensional discrete Fourier transform for real input. This function is a wrapper for pyfftw.interfaces.numpy_fft.irfftn, with an interface similar to that of numpy.fft.irfftn.

Parameters:
a : array_like

Input array

s : sequence of ints

Shape of the output along each transformed axis (input is cropped or zero-padded to match). This parameter is not optional because, unlike ifftn, the output shape cannot be uniquely determined from the input shape.

axes : sequence of ints, optional (default None)

Axes over which to compute the inverse DFT.

Returns:
af : ndarray

Inverse DFT of input array

sporco.linalg.dctii(x, axes=None)[source]

Compute a multi-dimensional DCT-II over specified array axes. This function is implemented by calling the one-dimensional DCT-II scipy.fftpack.dct with normalization mode ‘ortho’ for each of the specified axes.

Parameters:
a : array_like

Input array

axes : sequence of ints, optional (default None)

Axes over which to compute the DCT-II.

Returns:
y : ndarray

DCT-II of input array

sporco.linalg.idctii(x, axes=None)[source]

Compute a multi-dimensional inverse DCT-II over specified array axes. This function is implemented by calling the one-dimensional inverse DCT-II scipy.fftpack.idct with normalization mode ‘ortho’ for each of the specified axes.

Parameters:
a : array_like

Input array

axes : sequence of ints, optional (default None)

Axes over which to compute the inverse DCT-II.

Returns:
y : ndarray

Inverse DCT-II of input array

sporco.linalg.fftconv(a, b, axes=(0, 1))[source]

Compute a multi-dimensional convolution via the Discrete Fourier Transform.

Parameters:
a : array_like

Input array

b : array_like

Input array

axes : sequence of ints, optional (default (0, 1))

Axes on which to perform convolution

Returns:
ab : ndarray

Convolution of input arrays, a and b, along specified axes

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

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

Parameters:
x : array_like

Input array x

y : array_like

Input array y

axis : int, optional (default -1)

Axis over which to compute the sum

Returns:
y : ndarray

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

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

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 [32])

\[(\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:
ah : array_like

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

rho : float

Linear system parameter \(\rho\)

b : array_like

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

c : array_like, optional (default None)

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

axis : int, optional (default 4)

Axis along which to solve the linear system

Returns:
x : ndarray

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

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

Compute cached component used by solvedbi_sm.

Parameters:
ah : array_like

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

a : array_like

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

rho : float

Linear system parameter \(\rho\)

axis : int, optional (default 4)

Axis along which to solve the linear system

Returns:
c : ndarray

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

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

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 [32])

\[(\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:
ah : array_like

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

d : array_like

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

b : array_like

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

c : array_like, optional (default None)

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

axis : int, optional (default 4)

Axis along which to solve the linear system

Returns:
x : ndarray

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

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

Compute cached component used by solvedbd_sm.

Parameters:
ah : array_like

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

a : array_like

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

d : array_like

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

axis : int, optional (default 4)

Axis along which to solve the linear system

Returns:
c : ndarray

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

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

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 [32])

\[(\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:
ah : array_like

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

rho : float

Linear system parameter \(\rho\)

b : array_like

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

axisM : int

Axis in input corresponding to index m in linear system

axisK : int

Axis in input corresponding to index k in linear system

Returns:
x : ndarray

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

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

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 inner products are taken.

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

\[(\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:
ah : array_like

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

rho : float

Linear system parameter \(\rho\)

b : array_like

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

axisK : int

Axis in input corresponding to index k in linear system

dimN : int, optional (default 2)

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

Returns:
x : ndarray

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

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

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 [32])

\[(\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.

Parameters:
ah : array_like

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

rho : float

Parameter rho

b : array_like

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

axisM : int

Axis in input corresponding to index m in linear system

axisK : int

Axis in input corresponding to index k in linear system

tol : float

CG tolerance

mit : int

CG maximum iterations

isn : array_like

CG initial solution

Returns:
x : ndarray

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

cgit : int

Number of CG iterations

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

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

Parameters:
A : array_like

Array \(A\)

rho : float

Scalar \(\rho\)

check_finite : bool, optional (default False)

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

Returns:
lu : ndarray

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

piv : ndarray

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)[source]

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:
A : array_like

Matrix \(A\)

rho : float

Scalar \(\rho\)

b : array_like

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

lu : array_like

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

piv : array_like

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

check_finite : bool, optional (default False)

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

Returns:
x : ndarray

Solution to the linear system

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

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:
A : array_like

Matrix \(A\)

rho : float

Scalar \(\rho\)

b : array_like

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

lu : array_like

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

piv : array_like

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

check_finite : bool, optional (default False)

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

Returns:
x : ndarray

Solution to the linear system

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

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

Parameters:
A : array_like

Array \(A\)

rho : float

Scalar \(\rho\)

lower : bool, optional (default False)

Flag indicating whether lower or upper triangular factors are computed

check_finite : bool, optional (default False)

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

Returns:
c : ndarray

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

lwr : bool

Flag indicating whether the factor is lower or upper triangular

sporco.linalg.cho_solve_ATAI(A, rho, b, c, lwr, check_finite=True)[source]

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:
A : array_like

Matrix \(A\)

rho : float

Scalar \(\rho\)

b : array_like

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

c : array_like

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

lwr : bool

Flag indicating whether the factor is lower or upper triangular

Returns:
x : ndarray

Solution to the linear system

sporco.linalg.cho_solve_AATI(A, rho, b, c, lwr, check_finite=True)[source]

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:
A : array_like

Matrix \(A\)

rho : float

Scalar \(\rho\)

b : array_like

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

c : array_like

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

lwr : bool

Flag indicating whether the factor is lower or upper triangular

Returns:
x : ndarray

Solution to the linear system

sporco.linalg.zpad(x, pd, ax)[source]

Zero-pad array x with pd = (leading, trailing) zeros on axis ax.

Parameters:
x : array_like

Array to be padded

pd : tuple

Sequence of two ints (leading,trailing) specifying number of zeros for padding

ax : int

Axis to be padded

Returns:
xp : array_like

Padded array

sporco.linalg.Gax(x, ax)[source]

Compute gradient of x along axis ax.

Parameters:
x : array_like

Input array

ax : int

Axis on which gradient is to be computed

Returns:
xg : ndarray

Output array

sporco.linalg.GTax(x, ax)[source]

Compute transpose of gradient of x along axis ax.

Parameters:
x : array_like

Input array

ax : int

Axis on which gradient transpose is to be computed

Returns:
xg : ndarray

Output array

sporco.linalg.GradientFilters(ndim, axes, axshp, dtype=None)[source]

Construct a set of filters for computing gradients in the frequency domain.

Parameters:
ndim : integer

Total number of dimensions in array in which gradients are to be computed

axes : tuple of integers

Axes on which gradients are to be computed

axshp : tuple of integers

Shape of axes on which gradients are to be computed

dtype : dtype

Data type of output arrays

Returns:
Gf : ndarray

Frequency domain gradient operators \(\hat{G}_i\)

GHGf : ndarray

Sum of products \(\sum_i \hat{G}_i^H \hat{G}_i\)

sporco.linalg.zdivide(x, y)[source]

Return x/y, with 0 instead of NaN where y is 0.

Parameters:
x : array_like

Numerator

y : array_like

Denominator

Returns:
z : ndarray

Quotient x/y

sporco.linalg.shrink1(x, alpha)[source]

Scalar shrinkage/soft thresholding function

\[\mathcal{S}_{1,\alpha}(\mathbf{x}) = \mathrm{sign}(\mathbf{x}) \odot \max(0, |\mathbf{x}| - \alpha) = \mathrm{prox}_f(\mathbf{x}) \;\; \text{where} \;\; f(\mathbf{u}) = \alpha \|\mathbf{u}\|_1 \;\;.\]
Parameters:
x : array_like

Input array \(\mathbf{x}\)

alpha : float or array_like

Shrinkage parameter \(\alpha\)

Returns:
x : ndarray

Output array

sporco.linalg.shrink2(x, alpha, axis=-1)[source]

Vector shrinkage/soft thresholding function

\[\mathcal{S}_{2,\alpha}(\mathbf{x}) = \frac{\mathbf{x}}{\|\mathbf{x}\|_2} \max(0, \|\mathbf{x}\|_2 - \alpha) = \mathrm{prox}_f(\mathbf{x}) \;\; \text{where} \;\; f(\mathbf{u}) = \alpha \|\mathbf{u}\|_2 \;\;.\]

The \(\ell_2\) norm is applied over the specified axis of a multi-dimensional input (the last axis by default).

Parameters:
x : array_like

Input array \(\mathbf{x}\)

alpha : float or array_like

Shrinkage parameter \(\alpha\)

axis : int, optional (default -1)

Axis of x over which the \(\ell_2\) norm

Returns:
x : ndarray

Output array

sporco.linalg.shrink12(x, alpha, beta, axis=-1)[source]

Compound shrinkage/soft thresholding function [38] [10]

\[\mathcal{S}_{1,2,\alpha,\beta}(\mathbf{x}) = \mathcal{S}_{2,\beta}(\mathcal{S}_{1,\alpha}(\mathbf{x})) = \mathrm{prox}_f(\mathbf{x}) \;\; \text{where} \;\; f(\mathbf{u}) = \alpha \|\mathbf{u}\|_1 + \beta \|\mathbf{u}\|_2 \;\;.\]

The \(\ell_2\) norm is applied over the specified axis of a multi-dimensional input (the last axis by default).

Parameters:
x : array_like

Input array \(\mathbf{x}\)

alpha : float or array_like

Shrinkage parameter \(\alpha\)

beta : float or array_like

Shrinkage parameter \(\beta\)

axis : int, optional (default -1)

Axis of x over which the \(\ell_2\) norm

Returns:
x : ndarray

Output array

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

Project \(\mathbf{b}\) into the \(\ell_2\) ball of radius \(r\) about \(\mathbf{s}\), i.e. \(\{ \mathbf{x} : \|\mathbf{x} - \mathbf{s} \|_2 \leq r \}\).

Parameters:
b : array_like

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

s : array_like

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

r : float

Radius of ball

axes : sequence of ints, optional (default all axes)

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

Returns:
x : ndarray

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

sporco.linalg.promote16(u, fn=None, *args, **kwargs)[source]

Utility function for use with functions that do not support arrays of dtype np.float16. This function has two distinct modes of operation. If called with only the u parameter specified, the returned value is either u itself if u is not of dtype np.float16, or u promoted to np.float32 dtype if it is. If the function parameter fn is specified then u is conditionally promoted as described above, passed as the first argument to function fn, and the returned values are converted back to dtype np.float16 if u is of that dtype. Note that if parameter fn is specified, it may not be be specified as a keyword argument if it is followed by any non-keyword arguments.

Parameters:
u : array_like

Array to be promoted to np.float32 if it is of dtype np.float16

fn : function or None, optional (default None)

Function to be called with promoted u as first parameter and *args and **kwargs as additional parameters

*args

Variable length list of arguments for function fn

**kwargs

Keyword arguments for function fn

Returns:
up : ndarray

Conditionally dtype-promoted version of u if fn is None, or value(s) returned by fn, converted to the same dtype as u, if fn is a function

sporco.linalg.atleast_nd(n, u)[source]

If the input array has fewer than n dimensions, append singleton dimensions so that it is n dimensional. Note that the interface differs substantially from that of numpy.atleast_3d etc.

Parameters:
n : int

Minimum number of required dimensions

u : array_like

Input array

Returns:
v : ndarray

Output array with at least n dimensions

sporco.linalg.split(u, axis=0)[source]

Split an array into a list of arrays on the specified axis. The length of the list is the shape of the array on the specified axis, and the corresponding axis is removed from each entry in the list. This function does not have the same behaviour as numpy.split.

Parameters:
u : array_like

Input array

axis : int, optional (default 0)

Axis on which to split the input array

Returns:
v : list of ndarray

List of arrays

sporco.linalg.blockcirculant(A)[source]

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

Parameters:
A : tuple of array_like

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

Returns:
B : ndarray

Output array

sporco.linalg.fl2norm2(xf, axis=(0, 1))[source]

Compute the squared \(\ell_2\) norm in the DFT domain, taking into account the unnormalised DFT scaling, i.e. given the DFT of a multi-dimensional array computed via fftn, return the squared \(\ell_2\) norm of the original array.

Parameters:
xf : array_like

Input array

axis : sequence of ints, optional (default (0,1))

Axes on which the input is in the frequency domain

Returns:
x : float

\(\|\mathbf{x}\|_2^2\) where the input array is the result of applying fftn to the specified axes of multi-dimensional array \(\mathbf{x}\)

sporco.linalg.rfl2norm2(xf, xs, axis=(0, 1))[source]

Compute the squared \(\ell_2\) norm in the DFT domain, taking into account the unnormalised DFT scaling, i.e. given the DFT of a multi-dimensional array computed via rfftn, return the squared \(\ell_2\) norm of the original array.

Parameters:
xf : array_like

Input array

xs : sequence of ints

Shape of original array to which rfftn was applied to obtain the input array

axis : sequence of ints, optional (default (0,1))

Axes on which the input is in the frequency domain

Returns:
x : float

\(\|\mathbf{x}\|_2^2\) where the input array is the result of applying rfftn to the specified axes of multi-dimensional array \(\mathbf{x}\)

sporco.linalg.rrs(ax, b)[source]

Compute relative residual \(\|\mathbf{b} - A \mathbf{x}\|_2 / \|\mathbf{b}\|_2\) of the solution to a linear equation \(A \mathbf{x} = \mathbf{b}\). Returns 1.0 if \(\mathbf{b} = 0\).

Parameters:
ax : array_like

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

b : array_like

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

Returns:
x : float

Relative residual