brahmap.base.LinearOperator

Bases: BaseLinearOperator

A concrete linear operator constructed from functional mappings.

This class provides the generic foundation for a linear operator for the matrix-vector multiplication operation. It requires a function matvec for the operation \(x \mapsto A(x)=Ax\), and an optional transposed-matrix-vector function rmatvec for the operation \(x \mapsto A(x)=A^T x\).

For symmetric operators (like \(P^T P\) or block-diagonal preconditioners), 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\), i.e. the number of columns of the operator	required
nargout	int	Size of the output vector \(A(x)\), i.e. the number of rows of the operator	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	{}

Attributes:

Name	Type	Description
dtype	dtype	The data type of the operator
nargin	int	Size of the input vector \(x\), i.e. the number of columns of the operator
nargout	int	Size of the output vector \(A(x)\), i.e. the number of rows of the operator
symmetric	bool	Indicates whether the operator is symmetric or not
shape	tuple[int, int]	A tuple (nargout, nargin) representing the shape of the operator
nMatvec	int	The number of matrix-vector multiplications computed so far
T	LinearOperator	The transpose of this linear operator
H	LinearOperator	The Hermitian adjoint of this linear operator
logger	Logger	The logger instance for this operator

Methods:

Name	Description
reset_counters	Resets matrix-vector product counter to zero.
dot	Numpy-like dot() method.
matvec	Matrix-vector multiplication method.
to_array	Returns the dense form of the linear operator as a 2D NumPy array.

Source code in brahmap/base/linop.py

class LinearOperator(BaseLinearOperator):
    """
    A concrete linear operator constructed from functional mappings.

    This class provides the generic foundation for a linear operator for
    the matrix-vector multiplication operation. It requires a function
    `matvec` for the operation $x \\mapsto A(x)=Ax$, and an optional
    transposed-matrix-vector function `rmatvec` for the operation
    $x \\mapsto A(x)=A^T x$.

    For symmetric operators (like $P^T P$ or block-diagonal
    preconditioners), `rmatvec` is ignored.

    All other keyword arguments are passed directly to the superclass.

    Parameters
    ----------
    nargin : int
        Size of the input vector $x$, i.e. the number of columns of the operator
    nargout : int
        Size of the output vector $A(x)$, i.e. the number of rows of the operator
    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

    Attributes
    ----------
    dtype : np.dtype
        The data type of the operator
    nargin : int
        Size of the input vector $x$, i.e. the number of columns of the operator
    nargout : int
        Size of the output vector $A(x)$, i.e. the number of rows of the operator
    symmetric : bool
        Indicates whether the operator is symmetric or not
    shape : tuple[int, int]
        A tuple `(nargout, nargin)` representing the shape of the operator
    nMatvec : int
        The number of matrix-vector multiplications computed so far
    T : LinearOperator
        The transpose of this linear operator
    H : LinearOperator
        The Hermitian adjoint of this linear operator
    logger : logging.Logger
        The logger instance for this operator
    """

    def __init__(
        self,
        nargin: int,
        nargout: int,
        matvec: Callable,
        rmatvec: Callable | None = None,
        **kwargs,
    ) -> None:
        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

        self.__H: LinearOperator | None = None

        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, LinearOperator):
                    self.__H = adjoint_of
                else:
                    msg = (
                        "kwarg adjoint_of / transpose_of must be of type"
                        " LinearOperator."
                    )
                    msg += " Got " + str(adjoint_of.__class__)
                    msg = f"{self.__class__.__name__}: " + msg
                    raise ValueError(msg)

    @property
    def T(self) -> "LinearOperator":
        """The transpose operator

        Returns
        -------
        LinearOperator
            The transpose of this linear operator
        """
        return cast(LinearOperator, self.__H)

    @property
    def H(self) -> "LinearOperator":
        """The adjoint operator

        Returns
        -------
        LinearOperator
            The Hermitian adjoint of this linear operator
        """
        return cast(LinearOperator, self.__H)

    def matvec(self, x) -> npt.NDArray[np.number]:
        """
        Matrix-vector multiplication method.

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

        Parameters
        ----------
        x : npt.NDArray[np.number]
            The input vector $x$ to be multiplied by the operator

        Returns
        -------
        npt.NDArray[np.number]
            The result of the matrix-vector multiplication $A(x)$
        """
        x = np.asanyarray(x, dtype=self.dtype)
        M, N = self.shape

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

        y = self.__matvec(x)

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

        return y

    def to_array(self) -> npt.NDArray[np.number]:
        """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, for a large linear operator, it can occupy an enormous
            amount of memory and crash your system. Don't use it unless you
            understand the risk!

        Returns
        -------
        npt.NDArray[np.number]
            The dense 2D array representation of the linear operator
        """
        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) -> "LinearOperator":
        # 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) -> "LinearOperator":
        # 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}"
            )
            msg = f"{self.__class__.__name__}: " + msg
            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) -> npt.NDArray[np.number]:
        # 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) -> "LinearOperator | npt.NDArray[np.number]":
        # 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(
                f"{self.__class__.__name__}: Invalid multiplier! Cannot multiply"
            )

    def __rmul__(self, x) -> "LinearOperator | npt.NDArray[np.number]":
        if np.isscalar(x):
            return self.__mul__(x)
        raise ValueError(
            f"{self.__class__.__name__}: Invalid operation! Cannot multiply"
        )

    def __add__(self, other) -> "LinearOperator":
        if not isinstance(other, BaseLinearOperator):
            raise ValueError(
                f"{self.__class__.__name__}: 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}"
            )
            msg = f"{self.__class__.__name__}: " + msg
            raise ShapeError(msg)

        other_op = cast(LinearOperator, other)

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

        def rmatvec(x):
            return self.H(x) + other_op.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) -> "LinearOperator":
        return self * (-1)  # type: ignore

    def __sub__(self, other) -> "LinearOperator":
        if not isinstance(other, BaseLinearOperator):
            raise ValueError(
                f"{self.__class__.__name__}: 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}"
            )
            msg = f"{self.__class__.__name__}: " + msg
            raise ShapeError(msg)

        other_op = cast(LinearOperator, other)

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

        def rmatvec(x):
            return self.H(x) - other_op.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) -> "LinearOperator":
        if isinstance(other, (numbers.Number, np.number)):
            return self * (1.0 / cast(Any, other))
        else:
            raise ValueError(
                f"{self.__class__.__name__}: Invalid operation! Cannot divide"
            )

    def __pow__(self, other) -> "LinearOperator":
        if not isinstance(other, int):
            raise ValueError(
                f"{self.__class__.__name__}: Can only raise to integer power"
            )
        if other < 0:
            raise ValueError(
                f"{self.__class__.__name__}: Can only raise to nonnegative power"
            )
        if self.nargin != self.nargout:
            raise ShapeError(
                f"{self.__class__.__name__}: 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)  # type: ignore

Attributes

dtype: npt.DTypeLike property writable

The data type of the operator.

Returns:

Type	Description
DTypeLike	The NumPy data type of the operator

nargin: int property

Size of the input vector \(x\), i.e. the number of columns of the operator

Returns:

Type	Description
int	The number of input columns

nargout: int property

Size of the output vector \(A(x)\), i.e. the number of rows of the operator

Returns:

Type	Description
int	The number of output rows

symmetric: bool property

Indicates whether the operator is symmetric or not

Returns:

Type	Description
bool	True if symmetric, False otherwise

shape: Tuple[int, int] property

A tuple (nargout, nargin) representing the shape of the operator

Returns:

Type	Description
tuple[int, int]	A tuple (nrows, ncols)

nMatvec: int property

The number of matrix-vector multiplications computed so far

Returns:

Type	Description
int	The number of matrix-vector multiplications performed

T: LinearOperator property

The transpose operator

Returns:

Type	Description
LinearOperator	The transpose of this linear operator

H: LinearOperator property

The adjoint operator

Returns:

Type	Description
LinearOperator	The Hermitian adjoint of this linear operator

Functions

reset_counters() -> None

Resets matrix-vector product counter to zero.

Source code in brahmap/base/linop.py

def reset_counters(self) -> None:
    """Resets matrix-vector product counter to zero."""
    self._nMatvec = 0

dot(x) -> npt.NDArray[np.number]

Numpy-like dot() method.

Parameters:

Name	Type	Description	Default
x	Any	The input vector or object to multiply with.	required

Returns:

Type	Description
NDArray[number]	The result of the dot product.

Source code in brahmap/base/linop.py

def dot(self, x) -> npt.NDArray[np.number]:
    """Numpy-like dot() method.

    Parameters
    ----------
    x : Any
        The input vector or object to multiply with.
    Returns
    -------
    npt.NDArray[np.number]
        The result of the dot product.
    """
    return self.__mul__(x)

matvec(x) -> npt.NDArray[np.number]

Matrix-vector multiplication method.

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

Parameters:

Name	Type	Description	Default
x	NDArray[number]	The input vector \(x\) to be multiplied by the operator	required

Returns:

Type	Description
NDArray[number]	The result of the matrix-vector multiplication \(A(x)\)

Source code in brahmap/base/linop.py

def matvec(self, x) -> npt.NDArray[np.number]:
    """
    Matrix-vector multiplication method.

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

    Parameters
    ----------
    x : npt.NDArray[np.number]
        The input vector $x$ to be multiplied by the operator

    Returns
    -------
    npt.NDArray[np.number]
        The result of the matrix-vector multiplication $A(x)$
    """
    x = np.asanyarray(x, dtype=self.dtype)
    M, N = self.shape

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

    y = self.__matvec(x)

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

    return y

to_array() -> npt.NDArray[np.number]

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, for a large linear operator, it can occupy an enormous amount of memory and crash your system. Don't use it unless you understand the risk!

Returns:

Type	Description
NDArray[number]	The dense 2D array representation of the linear operator

Source code in brahmap/base/linop.py

def to_array(self) -> npt.NDArray[np.number]:
    """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, for a large linear operator, it can occupy an enormous
        amount of memory and crash your system. Don't use it unless you
        understand the risk!

    Returns
    -------
    npt.NDArray[np.number]
        The dense 2D array representation of the linear operator
    """
    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