brahmap.base.BlockDiagonalLinearOperator

Bases: LinearOperator

Base class for a block-diagonal linear operator

Parameters:

Name	Type	Description	Default
block_list	list[LinearOperator]	A flat list of linear operators representing the individual diagonal blocks	required
**kwargs	Any	Extra keyword arguments	{}

Attributes:

Name	Type	Description
block_list	list[LinearOperator]	A flat list of linear operators representing the individual diagonal blocks
num_blocks	int	The total number of diagonal blocks in the operator
row_size	NDArray[integer]	Array containing the number of rows for each block
col_size	NDArray[integer]	Array containing the number of columns for each block

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/blkop.py

class BlockDiagonalLinearOperator(LinearOperator):
    """Base class for a block-diagonal linear operator

    Parameters
    ----------
    block_list : List[LinearOperator]
        A flat list of linear operators representing the individual diagonal blocks
    **kwargs : Any
        Extra keyword arguments

    Attributes
    ----------
    block_list : list[LinearOperator]
        A flat list of linear operators representing the individual diagonal blocks
    num_blocks : int
        The total number of diagonal blocks in the operator
    row_size : npt.NDArray[np.integer]
        Array containing the number of rows for each block
    col_size : npt.NDArray[np.integer]
        Array containing the number of columns for each block
    """

    def __init__(
        self,
        block_list: List[LinearOperator],
        **kwargs: Any,
    ) -> None:
        try:
            for block in block_list:
                __, __ = block.shape
        except (TypeError, AttributeError):
            raise ValueError("The `block_list` must be a flat list of linearoperators")

        self.__row_size = np.asarray(
            [block.shape[0] for block in block_list], dtype=int
        )
        self.__col_size = np.asarray(
            [block.shape[-1] for block in block_list], dtype=int
        )

        nargin = sum(self.__col_size)
        nargout = sum(self.__row_size)
        symmetric = reduce(
            lambda x, y: x and y, [block.symmetric for block in block_list]
        )
        dtype = np.result_type(*[block.dtype for block in block_list])

        self.__block_list = block_list

        # transpose operator
        blocks_list_transposed = [block.T for block in block_list]

        matvec = partial(
            self._mult,
            block_list=self.block_list,
            dtype=dtype,
        )
        rmatvec = partial(
            self._mult,
            block_list=blocks_list_transposed,
            dtype=dtype,
        )

        super(BlockDiagonalLinearOperator, self).__init__(
            nargin=nargin,
            nargout=nargout,
            symmetric=symmetric,
            matvec=matvec,
            rmatvec=rmatvec,
            dtype=dtype,
            **kwargs,
        )

    @property
    def block_list(self) -> List[LinearOperator]:
        """A list of linear operators representing the individual diagonal blocks.

        Returns
        -------
        list[LinearOperator]
            A list of linear operators representing the individual diagonal blocks
        """
        return self.__block_list

    @property
    def num_blocks(self) -> int:
        """The total number of diagonal blocks in the operator.

        Returns
        -------
        int
            The total number of diagonal blocks in the operator
        """
        return len(self.block_list)

    @property
    def row_size(self) -> npt.NDArray[np.integer]:
        """Array containing the number of rows for each block.

        Returns
        -------
        npt.NDArray[np.integer]
            Array containing the number of rows for each block
        """
        return self.__row_size

    @property
    def col_size(self) -> npt.NDArray[np.integer]:
        """Array containing the number of columns for each block.

        Returns
        -------
        npt.NDArray[np.integer]
            Array containing the number of columns for each block
        """
        return self.__col_size

    def __getitem__(self, idx):
        block_range = self.block_list[idx]
        if isinstance(idx, slice):
            return BlockDiagonalLinearOperator(
                block_list=block_range,
            )
        else:
            return block_range

    def _mult(
        self,
        vec: npt.NDArray[np.number],
        block_list: List,
        dtype,
    ) -> npt.NDArray[np.number]:
        nrows = sum(block.shape[0] for block in block_list)

        prod = np.zeros(nrows, dtype=dtype)

        start_row_idx = 0
        start_col_idx = 0
        for idx, block in enumerate(block_list):
            end_row_idx = start_row_idx + block.shape[0]
            end_col_idx = start_col_idx + block.shape[1]

            prod[start_row_idx:end_row_idx] = block * vec[start_col_idx:end_col_idx]

            start_row_idx = end_row_idx
            start_col_idx = end_col_idx

        return prod

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

block_list: List[LinearOperator] property

A list of linear operators representing the individual diagonal blocks.

Returns:

Type	Description
list[LinearOperator]	A list of linear operators representing the individual diagonal blocks

num_blocks: int property

The total number of diagonal blocks in the operator.

Returns:

Type	Description
int	The total number of diagonal blocks in the operator

row_size: npt.NDArray[np.integer] property

Array containing the number of rows for each block.

Returns:

Type	Description
NDArray[integer]	Array containing the number of rows for each block

col_size: npt.NDArray[np.integer] property

Array containing the number of columns for each block.

Returns:

Type	Description
NDArray[integer]	Array containing the number of columns for each block

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