Constructors¶

Module: pycsou.linop.base

Classes for constructing linear operators.

Interfaces with common scientific computing Python objects

PyLopLinearOperator(PyLop[, is_symmetric, …])

Construct a linear operator from a pylops.LinearOperator instance.

ExplicitLinearOperator(array[, is_symmetric])

Construct an explicit linear operator.

DenseLinearOperator(ndarray[, is_symmetric])

Construct a linear operator from a Numpy array.

SparseLinearOperator(spmatrix[, is_symmetric])

Construct a linear operator from a sparse Scipy matrix (scipy.sparse.spmatrix).

DaskLinearOperator(dask_array[, is_symmetric])

Construct a linear operator from a Dask array (dask.array.core.Array).

Basic operators

DiagonalOperator(diag)

Construct a diagonal operator.

PolynomialLinearOperator(LinOp, coeffs)

Polynomial linear operator \(P(L)\).

IdentityOperator(size[, dtype])

Square identity operator.

NullOperator(shape[, dtype])

Null operator.

Special Structures

LinOpStack(*linops, axis[, n_jobs, …])

Stack linear operators together.

LinOpVStack(*linops[, n_jobs, joblib_backend])

Alias for vertical stacking, equivalent to LinOpStack(*linops, axis=0).

LinOpHStack(*linops[, n_jobs, joblib_backend])

Alias for horizontal stacking, equivalent to LinOpStack(*linops, axis=1).

BlockOperator(linops[, n_jobs])

Construct a block operator from N lists of M linear operators each.

BlockDiagonalOperator(*linops[, n_jobs])

Construct a block diagonal operator from N linear operators.

KroneckerProduct(linop1, linop2)

Kronecker product \(\otimes\) of two operators.

KroneckerSum(linop1, linop2)

Kronecker sum \(\oplus\) of two operators.

KhatriRaoProduct(linop1, linop2)

Khatri-Rao product \(\circ\) of two operators.

class PyLopLinearOperator(PyLop: pylops.LinearOperator.LinearOperator, is_symmetric: bool = False, is_dense: bool = False, is_sparse: bool = False, lipschitz_cst: float = inf)[source]¶

Bases: pycsou.core.linop.LinearOperator

Construct a linear operator from a pylops.LinearOperator instance.

__init__(PyLop: pylops.LinearOperator.LinearOperator, is_symmetric: bool = False, is_dense: bool = False, is_sparse: bool = False, lipschitz_cst: float = inf)[source]¶

Parameters

PyLop (pylops.LinearOperator) – Pylops linear operator.
is_symmetric (bool) – Whether the linear operator is symmetric or not.
is_dense (bool) – If True, the linear operator is specified explicitly in terms of a Numpy array.
is_sparse (bool) – If True, the linear operator is specified explicitly in terms of a Scipy sparse matrix.
lipschitz_cst (float) – Lipschitz constant of the operator.

adjoint(y: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class ExplicitLinearOperator(array: Union[numpy.ndarray, scipy.sparse.base.spmatrix, dask.array.core.Array], is_symmetric: bool = False)[source]¶

Bases: pycsou.core.linop.LinearOperator

Construct an explicit linear operator.

Explicit operators can be built from a Numpy array/Scipy sparse matrix/Dask array. The array is stored in the attribute self.mat.

__init__(array: Union[numpy.ndarray, scipy.sparse.base.spmatrix, dask.array.core.Array], is_symmetric: bool = False)[source]¶

Parameters

array (Union[np.ndarray, sparse.spmatrix, da.core.Array]) – Numpy array, Scipy sparse matrix or Dask array from which to construct the linear operator.
is_symmetric (bool) – Whether the linear operator is symmetric or not.

adjoint(y: Union[numbers.Number, numpy.ndarray, dask.array.core.Array]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class DenseLinearOperator(ndarray: numpy.ndarray, is_symmetric: bool = False)[source]¶

Bases: pycsou.linop.base.ExplicitLinearOperator

Construct a linear operator from a Numpy array.

The array is stored in the attribute self.mat.

__init__(ndarray: numpy.ndarray, is_symmetric: bool = False)[source]¶

Parameters

ndarray (numpy.ndarray) – Numpy array from which to construct the linear operator.
is_symmetric (bool) – Whether the linear operator is symmetric or not.

class SparseLinearOperator(spmatrix: scipy.sparse.base.spmatrix, is_symmetric: bool = False)[source]¶

Bases: pycsou.linop.base.ExplicitLinearOperator

Construct a linear operator from a sparse Scipy matrix (scipy.sparse.spmatrix).

The array is stored in the attribute self.mat.

__init__(spmatrix: scipy.sparse.base.spmatrix, is_symmetric: bool = False)[source]¶

Parameters

spmatrix (scipy.sparse.spmatrix) – Scipy sparse matrix from which to construct the linear operator.
is_symmetric (bool) – Whether the linear operator is symmetric or not.

class DaskLinearOperator(dask_array: dask.array.core.Array, is_symmetric: bool = False)[source]¶

Bases: pycsou.linop.base.ExplicitLinearOperator

Construct a linear operator from a Dask array (dask.array.core.Array).

The array is stored in the attribute self.mat.

__init__(dask_array: dask.array.core.Array, is_symmetric: bool = False)[source]¶

Parameters

dask_array (dask.array.core.Array) – Dask array from which to construct the linear operator.
is_symmetric (bool) – Whether the linear operator is symmetric or not.

class LinOpStack(*linops, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.core.linop.LinearOperator, pycsou.core.map.DiffMapStack

Stack linear operators together.

This class constructs a linear operator by stacking multiple linear operators together, either vertically (axis=0) or horizontally (axis=1):

Vertical stacking: Consider a collection \(\{L_i:\mathbb{R}^{N}\to \mathbb{R}^{M_i}, i=1,\ldots, k\}\) of linear operators. Their vertical stacking is defined as the operator

\[\begin{split}V:\begin{cases}\mathbb{R}^{N}\to \mathbb{R}^{M_1}\times \cdots \times\mathbb{R}^{M_k}\\ \mathbf{x}\mapsto (L_1\mathbf{x},\ldots, L_k\mathbf{x}). \end{cases}\end{split}\]

The adjoint of \(V\) is moreover given by:

\[V^\ast(\mathbf{y}_1, \ldots, \mathbf{y}_k)=\sum_{i=1}^k L_i^\ast \mathbf{y}_i, \quad \forall (\mathbf{y}_1, \ldots, \mathbf{y}_k)\in \mathbb{R}^{M_1}\times \cdots \times\mathbb{R}^{M_k}.\]

The Lipschitz constant of the vertically stacked operator can be bounded by \(\sqrt{\sum_{i=1}^k \|L_i\|_2^2}\).
Horizontal stacking: Consider a collection \(\{L_i:\mathbb{R}^{N_i}\to \mathbb{R}^{M}, i=1,\ldots, k\}\) of linear operators. Their horizontal stacking is defined as the operator

\[\begin{split}H:\begin{cases}\mathbb{R}^{N_1}\times \cdots \times\mathbb{R}^{N_k}\to \mathbb{R}^{M}\\ (\mathbf{x}_1,\ldots, \mathbf{x}_k)\mapsto \sum_{i=1}^k L_i \mathbf{x}_i. \end{cases}\end{split}\]

The adjoint of \(H\) is moreover given by:

\[H^\ast(\mathbf{y})=(L_1^\ast \mathbf{y},\ldots, L_k^\ast \mathbf{y}) \quad \forall \mathbf{y}\in \mathbb{R}^{M}.\]

The Lipschitz constant of the horizontally stacked operator can be bounded by \({\max_{i=1}^k \|L_i\|_2}\).

Examples

We can form the 2D gradient operator by stacking two 1D derivative operators:

>>> from pycsou.linop.diff import FirstDerivative, Gradient
>>> x = np.linspace(-2.5, 2.5, 100)
>>> X,Y = np.meshgrid(x,x)
>>> Z = peaks(X, Y)
>>> D1 = FirstDerivative(size=Z.size, shape=Z.shape, axis=0, kind='centered')
>>> D2 = FirstDerivative(size=Z.size, shape=Z.shape, axis=1, kind='centered')
>>> G1 = LinOpStack(D1, D2, axis=0)
>>> G2 = Gradient(shape=Z.shape, kind='centered')
>>> Z_d = D2*Z.flatten()
>>> np.allclose(G1*Z.flatten(), G2 * Z.flatten())
True
>>> np.allclose(G1.adjoint(G1*Z.flatten()), G2.adjoint(G2 * Z.flatten()))
True
>>> G3 = LinOpStack(D1.H, D2.H, axis=1)
>>> np.allclose(G1.adjoint(G1*Z.flatten()), (G3 * G1) * Z.flatten())
True
>>> parG1 = LinOpStack(D1, D2, axis=0, n_jobs=-1)
>>> parG3 = LinOpStack(D1.H, D2.H, axis=1, n_jobs=-1)
>>> np.allclose(G1.adjoint(G1*Z.flatten()), parG1.adjoint(parG1*Z.flatten()))
True
>>> np.allclose((G3 * G1) * Z.flatten(), (parG3 * parG1) * Z.flatten())
True

See also

LinOpVStack, LinOpHStack

__init__(*linops, axis: int, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

linops (LinearOperator) – List of linear operators to stack.
axis – Stacking direction: 0 for vertical and 1 for horizontal stacking.
n_jobs (int) – Number of cores to be used for parallel evaluation of the linear operator stack and its adjoint. If n_jobs==1, the operator stack and its adjoint are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) – Joblib backend (more details here).

adjoint(y: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class LinOpVStack(*linops, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.linop.base.LinOpStack

Alias for vertical stacking, equivalent to LinOpStack(*linops, axis=0).

__init__(*linops, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

linops (LinearOperator) – List of linear operators to stack.
n_jobs (int) – Number of cores to be used for parallel evaluation of the linear operator stack and its adjoint. If n_jobs==1, the operator stack and its adjoint are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

class LinOpHStack(*linops, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Bases: pycsou.linop.base.LinOpStack

Alias for horizontal stacking, equivalent to LinOpStack(*linops, axis=1).

__init__(*linops, n_jobs: int = 1, joblib_backend: str = 'loky')[source]¶

Parameters

linops (LinearOperator) – List of linear operators to stack.
n_jobs (int) – Number of cores to be used for parallel evaluation of the linear operator stack and its adjoint. If n_jobs==1, the operator stack and its adjoint are evaluated sequentially, otherwise they are evaluated in parallel. Setting n_jobs=-1 uses all available cores.
joblib_backend (str) –
Joblib backend (more details here).

BlockOperator(linops: List[List[pycsou.core.linop.LinearOperator]], n_jobs: int = 1) → pycsou.linop.base.PyLopLinearOperator[source]¶

Construct a block operator from N lists of M linear operators each.

Parameters

linops (List[List[LinearOperator]]) – List of lists of linear operators to be combined in block fashion. Alternatively, numpy.ndarray or scipy.sparse.spmatrix can be passed in place of one or more operators.
n_jobs (int) – Number of processes used to evaluate the N operators in parallel using multiprocessing. If n_jobs=1 (default), work in serial mode.

Returns

Block linear operator.

Return type

PyLopLinearOperator

Examples

>>> from pycsou.linop.base import BlockOperator
>>> from pycsou.linop.diff import SecondDerivative
>>> Nv, Nh = 11, 21
>>> D2hop = SecondDerivative(size=Nv * Nh, shape=(Nv,Nh), axis=1)
>>> D2vop = SecondDerivative(size=Nv * Nh, shape=(Nv,Nh), axis=0)
>>> Dblock = BlockOperator([[D2vop, 0.5 * D2vop, - D2hop], [D2hop, 2 * D2hop, D2vop]])
>>> x = np.zeros((Nv, Nh)); x[int(Nv//2), int(Nh//2)] = 1; z = np.tile(x, (3,1)).flatten()
>>> np.allclose(Dblock(z), np.concatenate(((D2vop + 0.5 * D2vop - D2hop)(x.flatten()), (D2hop + 2 * D2hop + D2vop)(x.flatten()))))
True

Notes

In mathematics, a block or a partitioned matrix is a matrix that is interpreted as being broken into sections called blocks or submatrices. Similarly a block operator is composed of N sets of M linear operators each such that its application in forward mode leads to

\[\begin{split}\begin{bmatrix} \mathbf{L_{1,1}} & \mathbf{L_{1,2}} & \cdots & \mathbf{L_{1,M}} \\ \mathbf{L_{2,1}} & \mathbf{L_{2,2}} & \cdots & \mathbf{L_{2,M}} \\ \vdots & \vdots & \cdots & \vdots \\ \mathbf{L_{N,1}} & \mathbf{L_{N,2}} & \cdots & \mathbf{L_{N,M}} \\ \end{bmatrix} \begin{bmatrix} \mathbf{x}_{1} \\ \mathbf{x}_{2} \\ \vdots \\ \mathbf{x}_{M} \end{bmatrix} = \begin{bmatrix} \mathbf{L_{1,1}} \mathbf{x}_{1} + \mathbf{L_{1,2}} \mathbf{x}_{2} + \mathbf{L_{1,M}} \mathbf{x}_{M} \\ \mathbf{L_{2,1}} \mathbf{x}_{1} + \mathbf{L_{2,2}} \mathbf{x}_{2} + \mathbf{L_{2,M}} \mathbf{x}_{M} \\ \vdots \\ \mathbf{L_{N,1}} \mathbf{x}_{1} + \mathbf{L_{N,2}} \mathbf{x}_{2} + \mathbf{L_{N,M}} \mathbf{x}_{M} \\ \end{bmatrix}\end{split}\]

while its application in adjoint mode leads to

\[\begin{split}\begin{bmatrix} \mathbf{L_{1,1}}^\ast & \mathbf{L_{2,1}}^\ast & \cdots & \mathbf{L_{N,1}}^\ast \\ \mathbf{L_{1,2}}^\ast & \mathbf{L_{2,2}}^\ast & \cdots & \mathbf{L_{N,2}}^\ast \\ \vdots & \vdots & \cdots & \vdots \\ \mathbf{L_{1,M}}^\ast & \mathbf{L_{2,M}}^\ast & \cdots & \mathbf{L_{N,M}}^\ast \\ \end{bmatrix} \begin{bmatrix} \mathbf{y}_{1} \\ \mathbf{y}_{2} \\ \vdots \\ \mathbf{y}_{N} \end{bmatrix} = \begin{bmatrix} \mathbf{L_{1,1}}^\ast \mathbf{y}_{1} + \mathbf{L_{2,1}}^\ast \mathbf{y}_{2} + \mathbf{L_{N,1}}^\ast \mathbf{y}_{N} \\ \mathbf{L_{1,2}}^\ast \mathbf{y}_{1} + \mathbf{L_{2,2}}^\ast \mathbf{y}_{2} + \mathbf{L_{N,2}}^\ast \mathbf{y}_{N} \\ \vdots \\ \mathbf{L_{1,M}}^\ast \mathbf{y}_{1} + \mathbf{L_{2,M}}^\ast \mathbf{y}_{2} + \mathbf{L_{N,M}}^\ast \mathbf{y}_{N} \\ \end{bmatrix}\end{split}\]

The Lipschitz constant of the block operator can be bounded by \(\max_{j=1}^M\sqrt{\sum_{i=1}^N \|\mathbf{L}_{i,j}\|_2^2}\).

Warning

The parameter n_jobs is currently unused and is there for compatibility with the future API of PyLops. The code should be updated when the next version on PyLops is released.

See also

BlockDiagonalOperator, LinOpStack

BlockDiagonalOperator(*linops, n_jobs: int = 1) → pycsou.linop.base.PyLopLinearOperator[source]¶

Construct a block diagonal operator from N linear operators.

Parameters

linops (LinearOperator) – Linear operators forming the diagonal blocks. Alternatively, numpy.ndarray or scipy.sparse.spmatrix can be passed in place of one or more operators.
n_jobs (int) – Number of processes used to evaluate the N operators in parallel using multiprocessing. If n_jobs=1 (default), work in serial mode.

Returns

Block diagonal linear operator.

Return type

PyLopLinearOperator

Examples

>>> from pycsou.linop.base import BlockDiagonalOperator
>>> from pycsou.linop.diff import SecondDerivative
>>> Nv, Nh = 11, 21
>>> D2hop = SecondDerivative(size=Nv * Nh, shape=(Nv,Nh), axis=1)
>>> D2vop = SecondDerivative(size=Nv * Nh, shape=(Nv,Nh), axis=0)
>>> Dblockdiag = BlockDiagonalOperator(D2vop, 0.5 * D2vop, -1 * D2hop)
>>> x = np.zeros((Nv, Nh)); x[int(Nv//2), int(Nh//2)] = 1; z = np.tile(x, (3,1)).flatten()
>>> np.allclose(Dblockdiag(z), np.concatenate((D2vop(x.flatten()), 0.5 * D2vop(x.flatten()), - D2hop(x.flatten()))))
True

Notes

A block-diagonal operator composed of N linear operators is created such as its application in forward mode leads to

\[\begin{split}\begin{bmatrix} \mathbf{L_1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{L_2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{L_N} \end{bmatrix} \begin{bmatrix} \mathbf{x}_{1} \\ \mathbf{x}_{2} \\ \vdots \\ \mathbf{x}_{N} \end{bmatrix} = \begin{bmatrix} \mathbf{L_1} \mathbf{x}_{1} \\ \mathbf{L_2} \mathbf{x}_{2} \\ \vdots \\ \mathbf{L_N} \mathbf{x}_{N} \end{bmatrix}\end{split}\]

while its application in adjoint mode leads to

\[\begin{split}\begin{bmatrix} \mathbf{L_1}^\ast & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{L_2}^\ast & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{L_N}^\ast \end{bmatrix} \begin{bmatrix} \mathbf{y}_{1} \\ \mathbf{y}_{2} \\ \vdots \\ \mathbf{y}_{N} \end{bmatrix} = \begin{bmatrix} \mathbf{L_1}^\ast \mathbf{y}_{1} \\ \mathbf{L_2}^\ast \mathbf{y}_{2} \\ \vdots \\ \mathbf{L_N}^\ast \mathbf{y}_{N} \end{bmatrix}\end{split}\]

The Lipschitz constant of the block-diagonal operator can be bounded by \({\max_{i=1}^N \|\mathbf{L}_{i}\|_2}\).

Warning

The parameter n_jobs is currently unused and is there for compatibility with the future API of PyLops. The code should be updated when the next version on PyLops is released.

See also

BlockOperator, LinOpStack

class DiagonalOperator(diag: Union[numbers.Number, numpy.ndarray])[source]¶

Bases: pycsou.core.linop.LinearOperator

Construct a diagonal operator.

__init__(diag: Union[numbers.Number, numpy.ndarray])[source]¶

Parameters: diag (Union[Number, np.ndarray]) – Diagonal of the operator.

adjoint(y: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class IdentityOperator(size: int, dtype: Optional[type] = None)[source]¶

Bases: pycsou.linop.base.DiagonalOperator

Square identity operator.

__init__(size: int, dtype: Optional[type] = None)[source]¶

Parameters

size (int) – Dimension of the domain.
dtype (Optional[type]) – Data type of the operator.

class NullOperator(shape: Tuple[int, int], dtype: Optional[type] = <class 'numpy.float64'>)[source]¶

Bases: pycsou.core.linop.LinearOperator

Null operator.

__init__(shape: Tuple[int, int], dtype: Optional[type] = <class 'numpy.float64'>)[source]¶

Parameters

shape (Tuple[int, int]) – Shape of the linear operator.
dtype (Optional[type]) – Data type of the linear operator.
is_explicit (bool) – If True, the linear operator is specified explicitly in terms of a Numpy/Scipy/Dask array.
is_dense (bool) – If True, the linear operator is specified explicitly in terms of a Numpy array.
is_sparse (bool) – If True, the linear operator is specified explicitly in terms of a Scipy sparse matrix.
is_dask (bool) – If True, the linear operator is specified explicitly in terms of a Dask array.
is_symmetric (bool) – Whether the linear operator is symmetric or not.
lipschitz_cst (float) – Lipschitz constant (maximal singular value) of the linear operator if known. Default to \(+\infty\).

adjoint(y: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

eigenvals(k: int, which='LM', **kwargs) → numpy.ndarray[source]¶

Find k eigenvalues of a square operator.

Parameters

k (int) – The number of eigenvalues and eigenvectors desired. k must be strictly smaller than the dimension of the operator. It is not possible to compute all eigenvectors of a matrix.
which (str, [‘LM’ | ‘SM’ | ‘LR’ | ‘SR’ | ‘LI’ | ‘SI’]) –
Which k eigenvalues to find:
- ‘LM’ : largest magnitude
- ‘SM’ : smallest magnitude
- ‘LR’ : largest real part
- ‘SR’ : smallest real part
- ‘LI’ : largest imaginary part
- ‘SI’ : smallest imaginary part
kwargs (dict) – A dict of additional keyword arguments values accepted by Scipy’s functions scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh().

Returns

Array containing the k requested eigenvalues.

Return type

np.ndarray

Raises

NotImplementedError – If the linear operator is not square.

Notes

This function calls one of the two Scipy’s functions: scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh(). See the documentation of these two functions for more information on their behaviour and the underlying ARPACK functions they rely on.

See also

svds()

singularvals(k: int, which='LM', **kwargs) → numpy.ndarray[source]¶

Compute the largest or smallest k singular values of an operator. The order of the singular values is not guaranteed.

Parameters

k (int) – Number of singular values to compute. Must be 1 <= k < min(self.shape).
which (str, [‘LM’ | ‘SM’]) –
Which k singular values to find:
- ‘LM’ : largest magnitude
- ‘SM’ : smallest magnitude
kwargs (dict) – A dict of additional keyword arguments values accepted by Scipy’s function scipy.sparse.linalg.svds().

Returns

Array containing the k requested singular values.

Return type

np.ndarray

Examples

>>> sig = np.repeat([0., 1., 0.], 10)
>>> filter = signal.hann(5); filter[filter.size//2:] = 0
>>> ConvOp = Convolve1D(size=sig.size, filter=filter)
>>> np.round((ConvOp.singularvals(k=3, which='LM', tol=1e-3)), 2)
array([0.5, 0.5, 0.5])

Notes

This function calls the Scipy’s function: scipy.sparse.linalg.svds(). See the documentation of this function for more information on its behaviour and the underlying ARPACK/LOBPCG functions it relies on.

See also

eigenvals()

class HomothetyMap(size: int, constant: numbers.Number)[source]¶

Bases: pycsou.linop.base.DiagonalOperator

__init__(size: int, constant: numbers.Number)[source]¶

Parameters: diag (Union[Number, np.ndarray]) – Diagonal of the operator.

jacobianT(arg: Optional[numbers.Number] = None) → numbers.Number[source]¶

Transpose of the Jacobian matrix of the differentiable map evaluated at arg.

Parameters: arg (Union[Number, np.ndarray]) – Point at which the transposed Jacobian matrix is evaluated.
Returns: Linear operator associated to the transposed Jacobian matrix.
Return type: LinearOperator

class PolynomialLinearOperator(LinOp: pycsou.core.linop.LinearOperator, coeffs: Union[numpy.ndarray, list, tuple])[source]¶

Bases: pycsou.core.linop.LinearOperator

Polynomial linear operator \(P(L)\).

Base class for polynomial operators. Useful for implementing generalised differential operators.

Given a polynomial \(P(x)=\sum_{k=0}^N a_k x^k\) and a square linear operator \(\mathbf{L}:\mathbb{R}^N\to \mathbb{R}^N,\) we define the polynomial linear operator \(P(\mathbf{L}):\mathbb{R}^N\to \mathbb{R}^N\) as:

\[P(\mathbf{L})=\sum_{k=0}^N a_k \mathbf{L}^k,\]

where \(\mathbf{L}^0\) is the identity matrix. The adjoint of \(P(\mathbf{L})\) is given by:

\[P(\mathbf{L})^\ast=\sum_{k=0}^N a_k (\mathbf{L}^\ast)^k.\]

Examples

>>> from pycsou.linop import DenseLinearOperator, PolynomialLinearOperator
>>> L = DenseLinearOperator(np.arange(64).reshape(8,8))
>>> PL = PolynomialLinearOperator(LinOp=L, coeffs=[1/2 ,2, 1])
>>> x = np.arange(8)
>>> np.allclose(PL(x), x/2 + 2 * L(x) + (L**2)(x))
True

__init__(LinOp: pycsou.core.linop.LinearOperator, coeffs: Union[numpy.ndarray, list, tuple])[source]¶

Parameters

LinOp (pycsou.core.LinearOperator) – Square linear operator \(\mathbf{L}\).
coeffs (Union[np.ndarray, list, tuple]) – Coefficients \(\{a_0,\ldots, a_N\}\) of the polynomial \(P\).

adjoint(x: Union[numbers.Number, numpy.ndarray]) → Union[numbers.Number, numpy.ndarray][source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class KroneckerProduct(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Bases: pycsou.core.linop.LinearOperator

Kronecker product \(\otimes\) of two operators.

Examples

>>> from pycsou.linop.base import KroneckerProduct
>>> from pycsou.linop.diff import SecondDerivative
>>> Nv, Nh = 11, 21
>>> D2hop = SecondDerivative(size=Nh)
>>> D2vop = SecondDerivative(size=Nv)
>>> Dkron = KroneckerProduct(D2hop, D2vop)
>>> x = np.zeros((Nv, Nh)); x[int(Nv//2), int(Nh//2)] = 1
>>> np.allclose(Dkron(x.flatten()), D2vop.apply_along_axis(D2hop.apply_along_axis(x.transpose(), axis=0).transpose(), axis=0).flatten())
True

Notes

The Kronecker product between two operators \(\mathbf{A}\in \mathbb{R}^{k\times l}\) and \(\mathbf{B}\in \mathbb{R}^{n\times m}\) is defined as:

\[\begin{split}\mathbf{A} \otimes \mathbf{B}=\left[ \begin{array}{ccc} A_{11}\mathbf{B} & \cdots & A_{1l}\mathbf{B} \\ \vdots & \ddots & \vdots \\ A_{k1}\mathbf{B} & \cdots & A_{kl}\mathbf{B} \\ \end{array} \right] \in \mathbb{R}^{kn\times lm}\end{split}\]

Let \(\mathbf{X}\in \mathbb{R}^{m\times l}\) and \(\mathbf{Y}\in \mathbb{R}^{n\times k}\). Then we have:

\[(\mathbf{A} \otimes \mathbf{B})\mbox{vec}(\mathbf{X})= \mbox{vec}\left(\mathbf{B}\mathbf{X}\mathbf{A}^T\right)\]

and

\[(\mathbf{A} \otimes \mathbf{B})^\ast\mbox{vec}(\mathbf{Y})= \mbox{vec}\left(\mathbf{B}^\ast\mathbf{Y}\overline{\mathbf{A}}\right)\]

where \(\mbox{vec}\) denotes the vectorisation operator. Such operations are leveraged to implement the linear operator in matrix-free form (i.e. the matrix \(\mathbf{A} \otimes \mathbf{B}\) is not explicitely constructed) both in forward and adjoint mode.

We have also \(\|\mathbf{A} \otimes \mathbf{B}\|_2=\|\mathbf{A}\|_2\|\mathbf{B}\|_2\) and \((\mathbf{A} \otimes \mathbf{B})^\dagger= \mathbf{A}^\dagger \otimes \mathbf{B}^\dagger\) which we use to compute efficiently self.lipschitz_cst and self.PinvOp.

See also

KroneckerSum, KhatriRaoProduct

__init__(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Parameters

linop1 (LinearOperator) – Linear operator on the left hand-side of the Kronecker product (multiplicand).
linop2 (LinearOperator) – Linear operator on the right hand-side of the Kronecker product (multiplier).

adjoint(y: numpy.ndarray) → numpy.ndarray[source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

property PinvOp¶: Return the pseudo-inverse of the operator as a LinearOperator.

class KroneckerSum(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Bases: pycsou.core.linop.LinearOperator

Kronecker sum \(\oplus\) of two operators.

Examples

>>> from pycsou.linop.base import KroneckerProduct, KroneckerSum, DiagonalOperator
>>> m1=np.linspace(0,3,5); m2=np.linspace(-3,2,7)
>>> D1=DiagonalOperator(diag=m1); ExpD1=DiagonalOperator(diag=np.exp(m1))
>>> D2=DiagonalOperator(diag=m2); ExpD2=DiagonalOperator(diag=np.exp(m2))
>>> Expkronprod=KroneckerProduct(ExpD1,ExpD2)
>>> Kronsum=KroneckerSum(D1,D2)
>>> np.allclose(np.diag(Expkronprod.todense().mat), np.exp(np.diag(Kronsum.todense().mat)))
True

Notes

The Kronecker sum between two operators \(\mathbf{A}\in \mathbb{R}^{k\times l}\) and \(\mathbf{B}\in \mathbb{R}^{n\times m}\) is defined as:

\[\mathbf{A} \oplus \mathbf{B}=\mathbf{A} \otimes \mathbf{I}_{n\times m} + \mathbf{I}_{k\times l} \otimes \mathbf{B} \in \mathbb{R}^{kn\times lm}.\]

Let \(\mathbf{X}\in \mathbb{R}^{m\times l}\) and \(\mathbf{Y}\in \mathbb{R}^{n\times k}\). Then we have:

\[(\mathbf{A} \oplus \mathbf{B})\mbox{vec}(\mathbf{X})= \mbox{vec}\left(\mathbf{X}\mathbf{A}^T + \mathbf{B}\mathbf{X}\right)\]

and

\[(\mathbf{A} \oplus \mathbf{B})^\ast\mbox{vec}(\mathbf{Y})= \mbox{vec}\left(\mathbf{Y}\overline{\mathbf{A}} + \mathbf{B}^\ast\mathbf{Y}\right)\]

where \(\mbox{vec}\) denotes the vectorisation operator. Such operations are leveraged to implement the linear operator in matrix-free form (i.e. the matrix \(\mathbf{A} \oplus \mathbf{B}\) is not explicitely constructed) both in forward and adjoint mode.

The Lipschitz constant of the Kronecker sum can be bounded by \(\|\mathbf{A}\|_2+ \|\mathbf{B}\|_2\).

See also

KroneckerSum, KhatriRaoProduct

__init__(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Parameters

linop1 (LinearOperator) – Linear operator on the left hand-side of the Kronecker sum.
linop2 (LinearOperator) – Linear operator on the right hand-side of the Kronecker sum.

adjoint(y: numpy.ndarray) → numpy.ndarray[source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]

class KhatriRaoProduct(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Bases: pycsou.core.linop.LinearOperator

Khatri-Rao product \(\circ\) of two operators.

Examples

>>> from pycsou.linop.base import KhatriRaoProduct
>>> from pycsou.linop.diff import SecondDerivative
>>> D1 = SecondDerivative(size=11)
>>> D2 = SecondDerivative(size=11)
>>> Dkrao = KhatriRaoProduct(D1, D2)
>>> x = np.arange(11)
>>> Dkrao(x).shape
(121,)
>>> np.allclose(Dkrao(x), ((D1.todense().mat * x[None, :]) @ D2.todense().mat.transpose()).flatten())
True

Notes

The Khatri-Rao product between two operators \(\mathbf{A}\in \mathbb{R}^{k\times l}\) and \(\mathbf{B}\in \mathbb{R}^{n\times l}\) is defined as the column-wise Kronecker product:

\[\mathbf{A} \circ \mathbf{B}=\left[ \begin{array}{ccc} \mathbf{A}_1\otimes \mathbf{B}_1 & \cdots & \mathbf{A}_l\otimes \mathbf{B}_l \end{array} \right] \in \mathbb{R}^{kn\times l}\]

Let \(\mathbf{x}\in \mathbb{R}^{l}\) and \(\mathbf{Y}\in \mathbb{R}^{n\times k}\). Then we have:

\[(\mathbf{A} \circ \mathbf{B})\mathbf{x}= \mbox{vec}\left(\mathbf{B}\mbox{diag}(\mathbf{x})\mathbf{A}^T\right)\]

and

\[(\mathbf{A} \circ \mathbf{B})^\ast\mbox{vec}(\mathbf{Y})= \mbox{diag}\left(\mathbf{B}^\ast\mathbf{Y}\overline{\mathbf{A}}\right)\]

where \(\mbox{diag}\), \(\mbox{vec}\) denote the diagonal and vectorisation operators respectively. Such operations are leveraged to implement the linear operator in matrix-free form (i.e. the matrix \(\mathbf{A} \circ \mathbf{B}\) is not explicitely constructed) both in forward and adjoint mode.

The Lipschitz constant of the Khatri-Rao product can be bounded by \(\|\mathbf{A}\|_2\|\mathbf{B}\|_2\).

See also

KroneckerProduct, KroneckerSum

__init__(linop1: pycsou.core.linop.LinearOperator, linop2: pycsou.core.linop.LinearOperator)[source]¶

Parameters

linop1 (LinearOperator) – Linear operator on the left hand-side of the Khatri-Rao product (multiplicand).
linop2 (LinearOperator) – Linear operator on the right hand-side of the Khatri-Rao product (multiplier).

Raises

ValueError – If linop1.shape[1] != self.linop2.shape[1].

adjoint(y: numpy.ndarray) → numpy.ndarray[source]¶

Evaluates the adjoint of the operator at a point.

Parameters: y (Union[Number, np.ndarray]) – Point at which the adjoint should be evaluated.
Returns: Evaluation of the adjoint at y.
Return type: Union[Number, np.ndarray]