brahmap.base.BlockLinearOperator

Bases: LinearOperator

A block-linear operator where each block refers to a LinearOperator.

The input parameter blocks must be a list of lists describing the blocks row-wise.

If there is only one block row, it should be specified as [[b1, b2, ..., bn]], not as [b1, b2, ..., bn].

If the overall linear operator is symmetric, only its upper triangle need be specified, e.g., [[A,B,C], [D,E], [F]], and the blocks on the diagonal must be square and symmetric.

Parameters:

Name	Type	Description	Default
blocks	list[list[LinearOperator]]	A nested list of linear operators defining the block structure	required
symmetric	bool	Indicates if the linear operator is symmetric, by default False	False
**kwargs	Any	Extra keyword arguments	{}

Attributes:

Name	Type	Description
blocks	list[list[LinearOperator]]	The list of blocks defining the block 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/blkop.py

class BlockLinearOperator(LinearOperator):
    """
    A block-linear operator where each block refers to a
    [`LinearOperator`][..LinearOperator].

    The input parameter `blocks` must be a list of lists describing the blocks
    row-wise.

    If there is only one block row, it should be specified as
    `[[b1, b2, ..., bn]]`, not as `[b1, b2, ..., bn]`.

    If the overall linear operator is symmetric, only its upper triangle
    need be specified, e.g., `[[A,B,C], [D,E], [F]]`, and the blocks on the
    diagonal must be square and symmetric.

    Parameters
    ----------
    blocks : List[List[LinearOperator]]
        A nested list of linear operators defining the block structure
    symmetric : bool, optional
        Indicates if the linear operator is symmetric, by default `False`
    **kwargs : Any
        Extra keyword arguments

    Attributes
    ----------
    blocks : list[list[LinearOperator]]
        The list of blocks defining the block operator
    """

    def __init__(
        self,
        blocks: List[List[LinearOperator]],
        symmetric: bool = False,
        **kwargs: Any,
    ) -> None:
        # If building a symmetric operator, fill in the blanks.
        # They're just references to existing objects.

        try:
            for block_row in blocks:
                for block_col in block_row:
                    __ = block_col.shape
        except (TypeError, AttributeError):
            raise ValueError("blocks should be a nested list of operators")

        if symmetric:
            nrow = len(blocks)
            ncol = len(blocks[0])
            if nrow != ncol:
                raise ShapeError("Inconsistent shape.")

            for block_row in blocks:
                if not block_row[0].symmetric:
                    raise ValueError("Blocks on diagonal must be symmetric.")

            self._blocks = blocks[:]
            for i in range(1, nrow):
                for j in range(i - 1, -1, -1):
                    self._blocks[i].insert(0, self._blocks[j][i].T)

        else:
            self._blocks = blocks

        log = kwargs.get("logger", null_log)
        log.debug("Building new BlockLinearOperator")

        nargins = [[blk.shape[-1] for blk in row] for row in self._blocks]
        log.debug("nargins = " + repr(nargins))
        nargins_by_row = [nargin[0] for nargin in nargins]
        if min(nargins_by_row) != max(nargins_by_row):
            raise ShapeError("Inconsistent block shapes")

        nargouts = [[blk.shape[0] for blk in row] for row in self._blocks]
        log.debug("nargouts = " + repr(nargouts))
        for row in nargouts:
            if min(row) != max(row):
                raise ShapeError("Inconsistent block shapes")

        nargin = sum(nargins[0])
        nargout = sum(out[0] for out in nargouts)

        # Create blocks of transpose operator.
        blocksT = [[blk.T for blk in row] for row in zip(*self._blocks)]

        def blk_matvec(x, blks):
            nargins = [[blk.shape[-1] for blk in blkrow] for blkrow in blks]
            nargouts = [[blk.shape[0] for blk in blkrow] for blkrow in blks]
            nargin = sum(nargins[0])
            nargout = sum(out[0] for out in nargouts)
            nx = len(x)
            self.logger.debug("Multiplying with a vector of size %d" % nx)
            self.logger.debug("nargin=%d, nargout=%d" % (nargin, nargout))
            if nx != nargin:
                raise ShapeError("Multiplying with vector of wrong shape.")

            result_type = np.result_type(self.dtype, x.dtype)
            y = np.zeros(nargout, dtype=result_type)

            nblk_row = len(blks)
            nblk_col = len(blks[0])

            row_start = col_start = 0
            for row in range(nblk_row):
                row_end = row_start + nargouts[row][0]
                yout = y[row_start:row_end]
                for col in range(nblk_col):
                    col_end = col_start + nargins[0][col]
                    xin = x[col_start:col_end]
                    B = blks[row][col]
                    yout[:] += B * xin
                    col_start = col_end
                row_start = row_end
                col_start = 0

            return y

        flat_blocks = [blk for row in blocks for blk in row]
        blk_dtypes = [blk.dtype for blk in flat_blocks]
        op_dtype = np.result_type(*blk_dtypes)

        super(BlockLinearOperator, self).__init__(
            nargin,
            nargout,
            symmetric=symmetric,
            matvec=lambda x: blk_matvec(x, self._blocks),
            rmatvec=lambda x: blk_matvec(x, blocksT),
            dtype=op_dtype,
            **kwargs,
        )

        cast(Any, self.H)._blocks = blocksT

    @property
    def blocks(self) -> List[List[LinearOperator]]:
        """The list of blocks defining the block operator.

        Returns
        -------
        list[list[LinearOperator]]
            A nested list of linear operators representing the block
            structure of the operator
        """
        return self._blocks

    def __getitem__(self, indices):
        blks = np.matrix(cast(Any, self._blocks), dtype=object)[indices]
        # If indexing narrowed it down to a single block, return it.
        if isinstance(blks, BaseLinearOperator):
            return blks
        # Otherwise, we have a matrix of blocks.
        return BlockLinearOperator(blks.tolist(), symmetric=False)

    def __contains__(self, op):
        flat_blocks = [blk for row in self.blocks for blk in row]
        return op in flat_blocks

    def __iter__(self):
        for block in self._blocks:
            yield block

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

blocks: List[List[LinearOperator]] property

The list of blocks defining the block operator.

Returns:

Type	Description
list[list[LinearOperator]]	A nested list of linear operators representing the block structure of the 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