Source code for pylops.basicoperators.Diagonal

import numpy as np

from pylops import LinearOperator
from pylops.utils.backend import get_array_module, to_cupy_conditional

[docs]class Diagonal(LinearOperator): r"""Diagonal operator. Applies element-wise multiplication of the input vector with the vector ``diag`` in forward and with its complex conjugate in adjoint mode. This operator can also broadcast; in this case the input vector is reshaped into its dimensions ``dims`` and the element-wise multiplication with ``diag`` is perfomed on the direction ``dir``. Note that the vector ``diag`` will need to have size equal to ``dims[dir]``. Parameters ---------- diag : :obj:`numpy.ndarray` Vector to be used for element-wise multiplication. dims : :obj:`list`, optional Number of samples for each dimension (``None`` if only one dimension is available) dir : :obj:`int`, optional Direction along which multiplication is applied. dtype : :obj:`str`, optional Type of elements in input array. Attributes ---------- shape : :obj:`tuple` Operator shape explicit : :obj:`bool` Operator contains a matrix that can be solved explicitly (``True``) or not (``False``) Notes ----- Element-wise multiplication between the model :math:`\mathbf{x}` and/or data :math:`\mathbf{y}` vectors and the array :math:`\mathbf{d}` can be expressed as .. math:: y_i = d_i x_i \quad \forall i=1,2,\ldots,N This is equivalent to a matrix-vector multiplication with a matrix containing the vector :math:`\mathbf{d}` along its main diagonal. For real-valued ``diag``, the Diagonal operator is self-adjoint as the adjoint of a diagonal matrix is the diagonal matrix itself. For complex-valued ``diag``, the adjoint is equivalent to the element-wise multiplication with the complex conjugate elements of ``diag``. """ def __init__(self, diag, dims=None, dir=0, dtype="float64"): ncp = get_array_module(diag) self.diag = diag.ravel() self.complex = True if ncp.iscomplexobj(self.diag) else False if dims is None: self.shape = (len(self.diag), len(self.diag)) self.dims = None self.reshape = False else: diagdims = [1] * len(dims) diagdims[dir] = dims[dir] self.diag = self.diag.reshape(diagdims) self.shape = (, self.dims = dims self.reshape = True self.dtype = np.dtype(dtype) self.explicit = False def _matvec(self, x): if type(self.diag) != type(x): self.diag = to_cupy_conditional(x, self.diag) if not self.reshape: y = self.diag * x.ravel() else: x = x.reshape(self.dims) y = self.diag * x return y.ravel() def _rmatvec(self, x): if type(self.diag) != type(x): self.diag = to_cupy_conditional(x, self.diag) if self.complex: diagadj = self.diag.conj() else: diagadj = self.diag if not self.reshape: y = diagadj * x.ravel() else: x = x.reshape(self.dims) y = diagadj * x return y.ravel()
[docs] def matrix(self): """Return diagonal matrix as dense :obj:`numpy.ndarray` Returns ------- densemat : :obj:`numpy.ndarray` Dense matrix. """ ncp = get_array_module(self.diag) densemat = ncp.diag(self.diag.squeeze()) return densemat
[docs] def todense(self): """Fast implementation of todense based on known structure of the operator Returns ------- densemat : :obj:`numpy.ndarray` Dense matrix. """ return self.matrix()