pylops.LinearOperator

class pylops.LinearOperator(Op=None, dtype=None, shape=None, dims=None, dimsd=None, clinear=None, explicit=None, name=None)[source]

Common interface for performing matrix-vector products.

This class is an overload of the scipy.sparse.linalg.LinearOperator class. It adds functionalities by overloading standard operators such as __truediv__ as well as creating convenience methods such as eigs, cond, and conj.

Note

End users of PyLops should not use this class directly but simply use operators that are already implemented. This class is meant for developers and it has to be used as the parent class of any new operator developed within PyLops. Find more details regarding implementation of new operators at Implementing new operators.

Parameters:
Op : scipy.sparse.linalg.LinearOperator or scipy.sparse.linalg._ProductLinearOperator or scipy.sparse.linalg._SumLinearOperator

Operator. If other arguments are provided, they will overwrite those obtained from Op.

shape : tuple(int, int), optional

Shape of operator. If not provided, obtained from dims and dimsd.

dims : tuple(int, ..., int), optional

New in version 2.0.0.

Dimensions of model. If not provided, (self.shape[1],) is used.

dimsd : tuple(int, ..., int), optional

New in version 2.0.0.

Dimensions of data. If not provided, (self.shape[0],) is used.

explicit : bool, optional

Operator contains a matrix that can be solved explicitly (True) or not (False). Defaults to False.

clinear : bool, optional

New in version 1.17.0.

Operator is complex-linear. Defaults to True.

name : str, optional

New in version 2.0.0.

Name of operator (to be used by pylops.utils.describe.describe)

Methods

__init__([Op, dtype, shape, dims, dimsd, …]) Initialize this LinearOperator.
adjoint() Hermitian adjoint.
apply_columns(cols) Apply subset of columns of operator
cond([uselobpcg]) Condition number of linear operator.
conj() Complex conjugate operator
div(y[, niter, densesolver]) Solve the linear problem \(\mathbf{y}=\mathbf{A}\mathbf{x}\).
dot(x) Matrix-matrix or matrix-vector multiplication.
eigs([neigs, symmetric, niter, uselobpcg]) Most significant eigenvalues of linear operator.
matmat(X) Matrix-matrix multiplication.
matvec(x) Matrix-vector multiplication.
reset_count() Reset counters
rmatmat(X) Matrix-matrix multiplication.
rmatvec(x) Adjoint matrix-vector multiplication.
todense([backend]) Return dense matrix.
toimag([forw, adj]) Imag operator
toreal([forw, adj]) Real operator
tosparse() Return sparse matrix.
trace([neval, method, backend]) Trace of linear operator.
transpose() Transpose this linear operator.
matvec(x)[source]

Matrix-vector multiplication.

Modified version of scipy matvec which does not consider the case where the input vector is np.matrix (the use np.matrix is now discouraged in numpy’s documentation).

Parameters:
x : numpy.ndarray

Input array of shape (N,) or (N,1)

Returns:
y : numpy.ndarray

Output array of shape (M,) or (M,1)

rmatvec(x)[source]

Adjoint matrix-vector multiplication.

Modified version of scipy rmatvec which does not consider the case where the input vector is np.matrix (the use np.matrix is now discouraged in numpy’s documentation).

Parameters:
y : numpy.ndarray

Input array of shape (M,) or (M,1)

Returns:
x : numpy.ndarray

Output array of shape (N,) or (N,1)

matmat(X)[source]

Matrix-matrix multiplication.

Modified version of scipy matmat which does not consider the case where the input vector is np.matrix (the use np.matrix is now discouraged in numpy’s documentation).

Parameters:
x : numpy.ndarray

Input array of shape (N,K)

Returns:
y : numpy.ndarray

Output array of shape (M,K)

rmatmat(X)[source]

Matrix-matrix multiplication.

Modified version of scipy rmatmat which does not consider the case where the input vector is np.matrix (the use np.matrix is now discouraged in numpy’s documentation).

Parameters:
x : numpy.ndarray

Input array of shape (M,K)

Returns:
y : numpy.ndarray

Output array of shape (N,K)

dot(x)[source]

Matrix-matrix or matrix-vector multiplication.

Parameters:
x : np.ndarray

Input array (or matrix)

Returns:
y : np.ndarray

Output array (or matrix) that represents the result of applying the linear operator on x.

div(y, niter=100, densesolver='scipy')[source]

Solve the linear problem \(\mathbf{y}=\mathbf{A}\mathbf{x}\).

Overloading of operator / to improve expressivity of Pylops when solving inverse problems.

Parameters:
y : np.ndarray

Data

niter : int, optional

Number of iterations (to be used only when explicit=False)

densesolver : str, optional

Use scipy (scipy) or numpy (numpy) dense solver

Returns:
xest : np.ndarray

Estimated model

todense(backend='numpy')[source]

Return dense matrix.

The operator is converted into its dense matrix equivalent. In order to do so, square or tall operators are applied to an identity matrix whose number of rows and columns is equivalent to the number of columns of the operator. Conversely, for skinny operators, the transpose operator is applied to an identity matrix whose number of rows and columns is equivalent to the number of rows of the operator and the resulting matrix is transposed (and complex conjugated).

Note that this operation may be costly for operators with large number of rows and columns and it should be used mostly as a way to inspect the structure of the matricial equivalent of the operator.

Parameters:
backend : str, optional

Backend used to densify matrix (numpy or cupy). Note that this must be consistent with how the operator has been created.

Returns:
matrix : numpy.ndarray or cupy.ndarray

Dense matrix.

tosparse()[source]

Return sparse matrix.

The operator in converted into its sparse (CSR) matrix equivalent. In order to do so, the operator is applied to series of unit vectors with length equal to the number of coloumns in the original operator.

Returns:
matrix : scipy.sparse.csr_matrix

Sparse matrix.

eigs(neigs=None, symmetric=False, niter=None, uselobpcg=False, **kwargs_eig)[source]

Most significant eigenvalues of linear operator.

Return an estimate of the most significant eigenvalues of the linear operator. If the operator has rectangular shape (shape[0]!=shape[1]), eigenvalues are first computed for the square operator \(\mathbf{A^H}\mathbf{A}\) and the square-root values are returned.

Parameters:
neigs : int

Number of eigenvalues to compute (if None, return all). Note that for explicit=False, only \(N-1\) eigenvalues can be computed where \(N\) is the size of the operator in the model space

symmetric : bool, optional

Operator is symmetric (True) or not (False). User should set this parameter to True only when it is guaranteed that the operator is real-symmetric or complex-hermitian matrices

niter : int, optional

Number of iterations for eigenvalue estimation

uselobpcg : bool, optional

Use scipy.sparse.linalg.lobpcg

**kwargs_eig

Arbitrary keyword arguments for scipy.sparse.linalg.eigs, scipy.sparse.linalg.eigsh, or scipy.sparse.linalg.lobpcg

Returns:
eigenvalues : numpy.ndarray

Operator eigenvalues.

Raises:
ValueError

If uselobpcg=True for a non-symmetric square matrix with complex type

Notes

Depending on the size of the operator, whether it is explicit or not and the number of eigenvalues requested, different algorithms are used by this routine.

More precisely, when only a limited number of eigenvalues is requested the scipy.sparse.linalg.eigsh method is used in case of symmetric=True and the scipy.sparse.linalg.eigs method is used symmetric=False. However, when the matrix is represented explicitly within the linear operator (explicit=True) and all the eigenvalues are requested the scipy.linalg.eigvals is used instead.

Finally, when only a limited number of eigenvalues is required, it is also possible to explicitly choose to use the scipy.sparse.linalg.lobpcg method via the uselobpcg input parameter flag.

Most of these algorithms are a port of ARPACK [1], a Fortran package which provides routines for quickly finding eigenvalues/eigenvectors of a matrix. As ARPACK requires only left-multiplication by the matrix in question, eigenvalues/eigenvectors can also be estimated for linear operators when the dense matrix is not available.

[1]http://www.caam.rice.edu/software/ARPACK/
cond(uselobpcg=False, **kwargs_eig)[source]

Condition number of linear operator.

Return an estimate of the condition number of the linear operator as the ratio of the largest and lowest estimated eigenvalues.

Parameters:
uselobpcg : bool, optional

Use scipy.sparse.linalg.lobpcg to compute eigenvalues

**kwargs_eig

Arbitrary keyword arguments for scipy.sparse.linalg.eigs, scipy.sparse.linalg.eigsh, or scipy.sparse.linalg.lobpcg

Returns:
eigenvalues : numpy.ndarray

Operator eigenvalues.

Notes

The condition number of a matrix (or linear operator) can be estimated as the ratio of the largest and lowest estimated eigenvalues:

\[k= \frac{\lambda_{max}}{\lambda_{min}}\]

The condition number provides an indication of the rate at which the solution of the inversion of the linear operator \(A\) will change with respect to a change in the data \(y\).

Thus, if the condition number is large, even a small error in \(y\) may cause a large error in \(x\). On the other hand, if the condition number is small then the error in \(x\) is not much bigger than the error in \(y\). A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned.

conj()[source]

Complex conjugate operator

Returns:
conjop : pylops.LinearOperator

Complex conjugate operator

apply_columns(cols)[source]

Apply subset of columns of operator

This method can be used to wrap a LinearOperator and mimic the action of a subset of columns of the operator on a reduced model in forward mode, and retrieve only the result of a subset of rows in adjoint mode.

Note that unless the operator has explicit=True, this is not optimal as the entire forward and adjoint passes of the original operator will have to be perfomed. It can however be useful for the implementation of solvers such as Orthogonal Matching Pursuit (OMP) that iteratively build a solution by evaluate only a subset of the columns of the operator.

Parameters:
cols : list

Columns to be selected

Returns:
colop : pylops.LinearOperator

Apply column operator

toreal(forw=True, adj=True)[source]

Real operator

Parameters:
forw : bool, optional

Apply real to output of forward pass

adj : bool, optional

Apply real to output of adjoint pass

Returns:
realop : pylops.LinearOperator

Real operator

toimag(forw=True, adj=True)[source]

Imag operator

Parameters:
forw : bool, optional

Apply imag to output of forward pass

adj : bool, optional

Apply imag to output of adjoint pass

Returns:
imagop : pylops.LinearOperator

Imag operator

trace(neval=None, method=None, backend='numpy', **kwargs_trace)[source]

Trace of linear operator.

Returns the trace (or its estimate) of the linear operator.

Parameters:
neval : int, optional

Maximum number of matrix-vector products compute. Default depends method.

method : str, optional

Should be one of the following:

Defaults to ‘explicit’ for explicit operators, and ‘Hutch++’ for the rest.

backend : str, optional

Backend used to densify matrix (numpy or cupy). Note that this must be consistent with how the operator has been created.

**kwargs_trace

Arbitrary keyword arguments passed to pylops.utils.trace_hutchinson, pylops.utils.trace_hutchpp, or pylops.utils.trace_nahutchpp

Returns:
trace : self.dtype

Operator trace.

Raises:
ValueError

If the operator has rectangular shape (shape[0] != shape[1])

NotImplementedError

If the method is not one of the available methods.

reset_count()[source]

Reset counters

When invoked all counters are set back to 0.