Skip to content

brahmap.base.InverseLO

Bases: LinearOperator

Constructs the inverse of a linear operator \(A\), represented as another linear operator.

This class wraps the inversion operator, applying an iterative solver provided with the method argument, whenever the inverse operator is multiplied with a vector.

Parameters:

Name Type Description Default
A LinearOperator

The primary linear operator or matrix

required
method Callable

The solver method to use (e.g., cg, pcg)

required
preconditioner LinearOperator | None

An optional preconditioner operator to accelerate convergence, by default None

None

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.

mult

Computes \(y = A^{-1}x\) by solving the linear system \(Ay = x\) for \(y\).

isconverged

Stores the convergence information depending on the exit status of the

Attributes:

Name Type Description
dtype DTypeLike

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 operator

H LinearOperator

The adjoint operator

method Callable

The solver method used to compute the inverse of \(A\).

converged int | None

Provides the solver convergence information.

preconditioner LinearOperator | None

The preconditioner linear operator for the iterative solver.

Source code in brahmap/base/linop.py
class InverseLO(LinearOperator):
    """Constructs the inverse of a linear operator $A$, represented as another linear
    operator.

    This class wraps the inversion operator, applying an iterative
    solver provided with the `method` argument, whenever the inverse
    operator is multiplied with a vector.

    Parameters
    ----------
    A : LinearOperator
        The primary linear operator or matrix
    method : Callable
        The solver method to use (e.g., `cg`, `pcg`)
    preconditioner : LinearOperator | None, optional
        An optional preconditioner operator to accelerate convergence, by default
        `None`
    """

    def __init__(
        self,
        A: LinearOperator,
        method: Callable,
        preconditioner: LinearOperator | None = None,
    ) -> None:
        super(InverseLO, self).__init__(
            nargin=A.shape[0], nargout=A.shape[1], matvec=self.mult, symmetric=True
        )
        self.A = A
        self.__method = method
        self.__preconditioner = preconditioner
        self.__converged = None

    def mult(self, x: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        """Computes $y = A^{-1}x$ by solving the linear system $Ay = x$ for $y$.

        This method uses the iterative solver routine (e.g., `scipy.sparse.linalg.cg`)
        specified during initialization as `method`.

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

        Returns
        -------
        npt.NDArray[np.number]
            The computed solution vector $y$
        """

        if self.method is None:
            raise ValueError(
                f"{self.__class__.__name__}: InverseLO solver method is not specified."
            )
        y, info = self.method(self.A, x, M=self.preconditioner)
        self.isconverged(info)
        return y

    def isconverged(self, info: int) -> None:
        """Stores the convergence information depending on the exit status of the
        solver.

        Parameters
        ----------
        info : int
            The output status code of the solver method
        """
        self.__converged = info

    @property
    def method(self) -> Callable:
        """The solver method used to compute the inverse of $A$."""
        return self.__method

    @property
    def converged(self) -> int | None:
        """Provides the solver convergence information.

        - `0` : Successful exit
        - `>0` : Convergence to tolerance not achieved, number of iterations
        - `<0` : Illegal input or breakdown
        """
        return self.__converged

    @property
    def preconditioner(self) -> LinearOperator | None:
        """The preconditioner linear operator for the iterative solver."""
        return self.__preconditioner

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

method: Callable property

The solver method used to compute the inverse of \(A\).

converged: int | None property

Provides the solver convergence information.

  • 0 : Successful exit
  • >0 : Convergence to tolerance not achieved, number of iterations
  • <0 : Illegal input or breakdown

preconditioner: LinearOperator | None property

The preconditioner linear operator for the iterative solver.

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

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

Computes \(y = A^{-1}x\) by solving the linear system \(Ay = x\) for \(y\).

This method uses the iterative solver routine (e.g., scipy.sparse.linalg.cg) specified during initialization as method.

Parameters:

Name Type Description Default
x NDArray[number]

The input vector \(x\) to be multiplied by the inverse operator

required

Returns:

Type Description
NDArray[number]

The computed solution vector \(y\)

Source code in brahmap/base/linop.py
def mult(self, x: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
    """Computes $y = A^{-1}x$ by solving the linear system $Ay = x$ for $y$.

    This method uses the iterative solver routine (e.g., `scipy.sparse.linalg.cg`)
    specified during initialization as `method`.

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

    Returns
    -------
    npt.NDArray[np.number]
        The computed solution vector $y$
    """

    if self.method is None:
        raise ValueError(
            f"{self.__class__.__name__}: InverseLO solver method is not specified."
        )
    y, info = self.method(self.A, x, M=self.preconditioner)
    self.isconverged(info)
    return y

isconverged(info: int) -> None

Stores the convergence information depending on the exit status of the solver.

Parameters:

Name Type Description Default
info int

The output status code of the solver method

required
Source code in brahmap/base/linop.py
def isconverged(self, info: int) -> None:
    """Stores the convergence information depending on the exit status of the
    solver.

    Parameters
    ----------
    info : int
        The output status code of the solver method
    """
    self.__converged = info