brahmap.base.LinearOperator

Bases: BaseLinearOperator

A generic linear operator class.

A linear operator constructed from a matrix-vector multiplication matvec, \(x \mapsto A(x)=Ax\) and possibly with a transposed-matrix-vector operation rmatvec, \(x \mapsto A(x)=A^T x\). If symmetric is True, rmatvec is ignored. All other keyword arguments are passed directly to the superclass.

Parameters:

Name	Type	Description	Default
nargin	int	Size of the input vector \(x\)	required
nargout	int	Size of the output vector \(A(x)\)	required
matvec	Callable	A function that defines the matrix-vector product \(x \mapsto A(x)=Ax\)	required
rmatvec	Optional[Callable]	A function that defines the transposed-matrix-vector product \(x \mapsto A(x)=A^T x\), by default None	None
**kwargs	Any	Extra keywords arguments	{}

Source code in brahmap/base/linop.py

class LinearOperator(BaseLinearOperator):
    """
    A generic linear operator class.

    A linear operator constructed from a matrix-vector multiplication `matvec`,
    $x \\mapsto A(x)=Ax$ and possibly with a transposed-matrix-vector
    operation `rmatvec`, $x \\mapsto A(x)=A^T x$. If `symmetric` is `True`,
    `rmatvec` is ignored. All other keyword arguments are passed directly to
    the superclass.

    Parameters
    ----------
    nargin : int
        Size of the input vector $x$
    nargout : int
        Size of the output vector $A(x)$
    matvec : Callable
        A function that defines the matrix-vector product $x \\mapsto A(x)=Ax$
    rmatvec : Optional[Callable], optional
        A function that defines the transposed-matrix-vector product
        $x \\mapsto A(x)=A^T x$, by default `None`
    **kwargs : Any
        Extra keywords arguments
    """

    def __init__(
        self,
        nargin: int,
        nargout: int,
        matvec: Callable,
        rmatvec: Optional[Callable] = None,
        **kwargs,
    ):
        super(LinearOperator, self).__init__(
            nargin,
            nargout,
            **kwargs,
        )
        adjoint_of = kwargs.get("adjoint_of", None) or kwargs.get("transpose_of", None)
        rmatvec = rmatvec or kwargs.get("matvec_transp", None)

        self.__matvec = matvec

        if self.symmetric:
            self.__H = self
        else:
            if adjoint_of is None:
                if rmatvec is not None:
                    # Create 'pointer' to transpose operator.
                    self.__H = LinearOperator(
                        nargout,
                        nargin,
                        matvec=rmatvec,
                        rmatvec=matvec,
                        adjoint_of=self,
                        **kwargs,
                    )
                else:
                    self.__H = None
            else:
                # Use operator supplied as transpose operator.
                if isinstance(adjoint_of, BaseLinearOperator):
                    self.__H = adjoint_of
                else:
                    msg = (
                        "kwarg adjoint_of / transpose_of must be of type"
                        " LinearOperator."
                    )
                    msg += " Got " + str(adjoint_of.__class__)
                    raise ValueError(msg)

    @property
    def T(self):
        """The transpose operator"""
        return self.__H

    @property
    def H(self):
        """The adjoint operator"""
        return self.__H

    def matvec(self, x):
        """
        Matrix-vector multiplication.

        The matvec property encapsulates the `matvec` routine specified at
        construct time, to ensure the consistency of the input and output
        arrays with the operator's shape.
        """
        x = np.asanyarray(x)
        M, N = self.shape

        # check input data consistency
        N = int(N)
        try:
            x = x.reshape(N)
        except ValueError:
            msg = (
                "The size of the input array is incompatible with the "
                "operator dimensions\n"
                f"size of the input array: {len(x)}\n"
                f"shape of the operator: {self.shape}"
            )
            raise ValueError(msg)

        y = self.__matvec(x)

        # check output data consistency
        M = int(M)
        try:
            y = y.reshape(M)
        except ValueError:
            msg = (
                "The size of the output array is incompatible with the "
                "operator dimensions\n"
                f"size of the output array: {len(y)}\n"
                f"shape of the operator: {self.shape}"
            )
            raise ValueError(msg)

        return y

    def to_array(self) -> np.ndarray:
        """Returns the dense form of the linear operator as a 2D NumPy array

        !!! Warning

            This method first allocates a NumPy array of shape `self.shape`
            and data-type `self.dtype`, and then fills them with numbers. As
            such it can occupy an enormous amount of memory. Don't use it
            unless you understand the risk!
        """
        n, m = self.shape
        H = np.empty((n, m), dtype=self.dtype)
        ej = np.zeros(m, dtype=self.dtype)
        for j in range(m):
            ej[j] = 1.0
            H[:, j] = self * ej
            ej[j] = 0.0
        return H

    def __mul_scalar(self, x):
        # Product between a linear operator and a scalar
        result_type = np.result_type(self.dtype, type(x))

        if x != 0:

            def matvec(y):
                return x * (self(y))

            def rmatvec(y):
                return x * (self.H(y))

            return LinearOperator(
                self.nargin,
                self.nargout,
                symmetric=self.symmetric,
                matvec=matvec,
                rmatvec=rmatvec,
                dtype=result_type,
            )
        else:
            return ZeroOperator(self.nargin, self.nargout, dtype=result_type)

    def __mul_linop(self, op):
        # Product between two linear operators
        if self.nargin != op.nargout:
            msg = (
                "Cannot multiply the two operators together\n"
                f"shape of the first operator: {self.shape}\n"
                f"shape of the second operator: {op.shape}"
            )
            raise ShapeError(msg)

        def matvec(x):
            return self(op(x))

        def rmatvec(x):
            return op.T(self.H(x))

        result_type = np.result_type(self.dtype, op.dtype)

        return LinearOperator(
            op.nargin,
            self.nargout,
            symmetric=False,  # Generally.
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __mul_vector(self, x):
        # Product between a linear operator and a vector
        self._nMatvec += 1
        result_type = np.result_type(self.dtype, x.dtype)
        return self.matvec(x).astype(result_type, copy=False)

    def __mul__(self, x):
        # Returns a linear operator if x is a scalar or a linear operator
        # Returns a vector if x is an array
        if isinstance(x, numbers.Number):
            return self.__mul_scalar(x)
        elif isinstance(x, BaseLinearOperator):
            return self.__mul_linop(x)
        elif isinstance(x, np.ndarray):
            return self.__mul_vector(x)
        else:
            raise ValueError("Invalid multiplier! Cannot multiply")

    def __rmul__(self, x):
        if np.isscalar(x):
            return self.__mul__(x)
        raise ValueError("Invalid operation! Cannot multiply")

    def __add__(self, other):
        if not isinstance(other, BaseLinearOperator):
            raise ValueError("Invalid operation! Cannot add")
        if self.shape != other.shape:
            msg = (
                "Cannot add the two operators together\n"
                f"shape of the first operator: {self.shape}\n"
                f"shape of the second operator: {other.shape}"
            )
            raise ShapeError(msg)

        def matvec(x):
            return self(x) + other(x)

        def rmatvec(x):
            return self.H(x) + other.T(x)

        result_type = np.result_type(self.dtype, other.dtype)

        return LinearOperator(
            self.nargin,
            self.nargout,
            symmetric=self.symmetric and other.symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __neg__(self):
        return self * (-1)

    def __sub__(self, other):
        if not isinstance(other, BaseLinearOperator):
            raise ValueError("Invalid operation! Cannot subtract")
        if self.shape != other.shape:
            msg = (
                "Cannot subtract one operator from the other\n"
                f"shape of the first operator: {self.shape}\n"
                f"shape of the second operator: {other.shape}"
            )
            raise ShapeError(msg)

        def matvec(x):
            return self(x) - other(x)

        def rmatvec(x):
            return self.H(x) - other.T(x)

        result_type = np.result_type(self.dtype, other.dtype)

        return LinearOperator(
            self.nargin,
            self.nargout,
            symmetric=self.symmetric and other.symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=result_type,
        )

    def __truediv__(self, other):
        if np.isscalar(other):
            return self * (1 / other)
        else:
            raise ValueError("Invalid operation! Cannot divide")

    def __pow__(self, other):
        if not isinstance(other, int):
            raise ValueError("Can only raise to integer power")
        if other < 0:
            raise ValueError("Can only raise to nonnegative power")
        if self.nargin != self.nargout:
            raise ShapeError("Can only raise square operators to a power")
        if other == 0:
            return IdentityOperator(self.nargin)
        if other == 1:
            return self
        return self * self ** (other - 1)

H property

The adjoint operator

T property

The transpose operator

matvec(x)

Matrix-vector multiplication.

The matvec property encapsulates the matvec routine specified at construct time, to ensure the consistency of the input and output arrays with the operator's shape.

Source code in brahmap/base/linop.py

def matvec(self, x):
    """
    Matrix-vector multiplication.

    The matvec property encapsulates the `matvec` routine specified at
    construct time, to ensure the consistency of the input and output
    arrays with the operator's shape.
    """
    x = np.asanyarray(x)
    M, N = self.shape

    # check input data consistency
    N = int(N)
    try:
        x = x.reshape(N)
    except ValueError:
        msg = (
            "The size of the input array is incompatible with the "
            "operator dimensions\n"
            f"size of the input array: {len(x)}\n"
            f"shape of the operator: {self.shape}"
        )
        raise ValueError(msg)

    y = self.__matvec(x)

    # check output data consistency
    M = int(M)
    try:
        y = y.reshape(M)
    except ValueError:
        msg = (
            "The size of the output array is incompatible with the "
            "operator dimensions\n"
            f"size of the output array: {len(y)}\n"
            f"shape of the operator: {self.shape}"
        )
        raise ValueError(msg)

    return y

to_array() -> np.ndarray

Returns the dense form of the linear operator as a 2D NumPy array

Warning

This method first allocates a NumPy array of shape self.shape and data-type self.dtype, and then fills them with numbers. As such it can occupy an enormous amount of memory. Don't use it unless you understand the risk!

Source code in brahmap/base/linop.py

def to_array(self) -> np.ndarray:
    """Returns the dense form of the linear operator as a 2D NumPy array

    !!! Warning

        This method first allocates a NumPy array of shape `self.shape`
        and data-type `self.dtype`, and then fills them with numbers. As
        such it can occupy an enormous amount of memory. Don't use it
        unless you understand the risk!
    """
    n, m = self.shape
    H = np.empty((n, m), dtype=self.dtype)
    ej = np.zeros(m, dtype=self.dtype)
    for j in range(m):
        ej[j] = 1.0
        H[:, j] = self * ej
        ej[j] = 0.0
    return H