brahmap.core.InvNoiseCovLO_Toeplitz01

Bases: InvNoiseCovLinearOperator

A linear operator representing the inverse of a Toeplitz noise covariance matrix \(N^{-1}\)

This inverse operator is based on the PCG based Toeplitz system solver.

The input covariance array must be at least of the size \(n\). The input power spectrum array must be of the size \(2n-2\) or \(2n-1\).

Parameters:

Name	Type	Description	Default
size	int	The size (dimension) of the linear operator	required
input	ArrayLike	The input array or data defining the operator	required
input_type	Literal['covariance', 'power_spectrum']	Specifies whether the input is a covariance array or a power spectrum array, by default "power_spectrum"	'power_spectrum'
precond_op	LinearOperator | Literal[None, 'Strang', 'TChan', 'RChan', 'KK2']	The preconditioner operator, by default None	None
precond_maxiter	int	The maximum number of iterations allowed for the preconditioner, by default 50	50
precond_atol	float	The absolute tolerance setting for the preconditioner solver, by default 1.0e-10	1e-10
precond_callback	Callable | None	An optional callback function executed within the preconditioner step, by default None	None
dtype	DTypeFloat	The data type of the operator, by default np.float64	float64

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.
get_inverse	Returns the inverse of this operator, which is the original

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
size	int	The dimension i.e. the number of rows/columns of the operator
precond_op	LinearOperator | None	The preconditioner operator used to accelerate convergence.
diag	NDArray[number]	The diagonal elements of the inverse noise covariance operator.
previous_num_iterations	int	The number of iterations performed in the last linear solve.

Source code in brahmap/core/noise_ops_toeplitz.py

class InvNoiseCovLO_Toeplitz01(InvNoiseCovLinearOperator):
    """A linear operator representing the inverse of a Toeplitz noise covariance
    matrix $N^{-1}$

    This inverse operator is based on the PCG based Toeplitz system solver.

    The input covariance array must be at least of the size $n$.
    The input power spectrum array must be of the size $2n-2$ or $2n-1$.

    Parameters
    ----------
    size : int
        The size (dimension) of the linear operator
    input : npt.ArrayLike
        The input array or data defining the operator
    input_type : Literal["covariance", "power_spectrum"], optional
        Specifies whether the `input` is a covariance array or a power
        spectrum array, by default `"power_spectrum"`
    precond_op : LinearOperator | Literal[None, "Strang", "TChan", "RChan", "KK2"], optional
        The preconditioner operator, by default `None`
    precond_maxiter : int, optional
        The maximum number of iterations allowed for the preconditioner, by
        default `50`
    precond_atol : float, optional
        The absolute tolerance setting for the preconditioner solver, by
        default `1.0e-10`
    precond_callback : Callable | None, optional
        An optional callback function executed within the preconditioner
        step, by default `None`
    dtype : DTypeFloat, optional
        The data type of the operator, by default `np.float64`
    """

    def __init__(
        self,
        size: int,
        input: npt.ArrayLike,
        input_type: Literal["covariance", "power_spectrum"] = "power_spectrum",
        precond_op: LinearOperator
        | Literal[None, "Strang", "TChan", "RChan", "KK2"] = None,
        precond_maxiter: int = 50,
        precond_atol: float = 1.0e-10,
        precond_callback: Callable | None = None,
        dtype: DTypeFloat = np.float64,
    ) -> None:
        self.__toeplitz_op = NoiseCovLO_Toeplitz01(
            size=size,
            input=input,
            input_type=input_type,
            dtype=dtype,
        )

        self.precond_atol = precond_atol
        self.precond_maxiter = precond_maxiter
        self.precond_callback = precond_callback

        self.__previous_num_iterations = 0

        super().__init__(
            nargin=size,
            matvec=self._mult,
            input_type=input_type,
            dtype=dtype,
        )

        if precond_op is None:
            self.precond_op = None
        elif isinstance(precond_op, LinearOperator) or isinstance(
            precond_op, np.ndarray
        ):
            self.precond_op = precond_op
        elif precond_op in ["Strang", "TChan", "RChan", "KK2"]:
            if input_type == "power_spectrum":
                cov = scipy.fft.ifft(
                    input,
                    workers=MPI_UTILS.nthreads_per_process,
                )[:size]
                cov = cov.real.astype(dtype=dtype, copy=False)  # type: ignore
            else:
                cov = np.asarray(input, dtype=dtype)[:size]

            if precond_op == "Strang":
                temp_size = int(np.floor(cov.size / 2))
                if cov.size % 2 == 0:
                    new_cov = np.concatenate(
                        [cov[:temp_size], cov[1 : temp_size + 1][::-1]]
                    )
                else:
                    new_cov = np.concatenate(
                        [cov[: temp_size + 1], cov[1 : temp_size + 1][::-1]]
                    )
            elif precond_op == "TChan":
                new_cov = np.empty_like(cov)
                new_cov[0] = cov[0]
                n = cov.size
                for idx in range(1, n):
                    new_cov[idx] = ((n - idx) * cov[idx] + idx * cov[n - idx]) / n
            elif precond_op == "RChan":
                new_cov = np.roll(np.flip(cov), 1)
                new_cov += cov
                new_cov[0] = cov[0]
            elif precond_op == "KK2":  # Circulant but not symmetric
                new_cov = np.roll(np.flip(cov), 1)
                new_cov[0] = 0
                new_cov = cov - new_cov

            self.precond_op = InvNoiseCovLO_Circulant(
                size=size,
                input=new_cov,
                input_type="covariance",
                dtype=dtype,
            )
        else:
            raise ValueError("Invalid preconditioner operator provided!")

    @property
    def precond_op(self) -> LinearOperator | None:
        """The preconditioner operator used to accelerate convergence.

        Returns
        -------
        LinearOperator | None
            The preconditioner operator, if set
        """
        return self.__precond_op

    @precond_op.setter
    def precond_op(self, operator: LinearOperator | None) -> None:
        if operator is not None:
            if self.shape != operator.shape:
                raise ValueError(
                    f"The shape of the input operator {operator.shape} is not "
                    f"compatible with the shape of inverse Toeplitz operator {self.shape}"
                )
        self.__precond_op = operator

    @property
    def diag(self) -> npt.NDArray[np.number]:
        """The diagonal elements of the inverse noise covariance operator.

        Returns
        -------
        npt.NDArray[np.number]
            A 1-d array containing the diagonal elements
        """
        try:
            diag_arr = getattr(self, "__diag")
        except AttributeError:
            factor = 1.0
            diag_arr = factor * np.ones(self.size, dtype=self.dtype)
        return diag_arr

    @diag.setter
    def diag(self, diag: npt.NDArray[np.number]) -> None:
        self.__diag = diag

    @property
    def previous_num_iterations(self) -> int:
        """The number of iterations performed in the last linear solve.

        Returns
        -------
        int
            Number of iterations
        """
        return self.__previous_num_iterations

    def get_inverse(self) -> "NoiseCovLO_Toeplitz01":  # type: ignore
        """Returns the inverse of this operator, which is the original
        noise covariance operator.

        Returns
        -------
        NoiseCovLO_Toeplitz01
            The noise covariance operator $N$
        """
        return self.__toeplitz_op

    def __callback(self, x, r, norm_residual) -> None:
        self.__previous_num_iterations += 1
        if self.precond_callback is not None:
            self.precond_callback(x, r, norm_residual)

    def _mult(self, vec: npt.NDArray[np.number]) -> npt.NDArray[np.number]:
        r"""Performs the matrix-vector product $y = N^{-1} v$ by iteratively
        solving $N y = v$ using PCG.

        Parameters
        ----------
        vec : npt.NDArray[np.number]
            The input vector $v$

        Returns
        -------
        npt.NDArray[np.number]
            The resulting vector
        """
        self.__previous_num_iterations = 0

        prod, _ = cg(
            A=self.__toeplitz_op,
            b=vec,
            atol=self.precond_atol,
            maxiter=self.precond_maxiter,
            M=self.precond_op,
            callback=self.__callback,
            parallel=False,
        )

        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

size: int property

The dimension i.e. the number of rows/columns of the operator

Returns:

Type	Description
int	The size of the operator

precond_op: LinearOperator | None property writable

The preconditioner operator used to accelerate convergence.

Returns:

Type	Description
LinearOperator | None	The preconditioner operator, if set

diag: npt.NDArray[np.number] property writable

The diagonal elements of the inverse noise covariance operator.

Returns:

Type	Description
NDArray[number]	A 1-d array containing the diagonal elements

previous_num_iterations: int property

The number of iterations performed in the last linear solve.

Returns:

Type	Description
int	Number of iterations

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

get_inverse() -> NoiseCovLO_Toeplitz01

Returns the inverse of this operator, which is the original noise covariance operator.

Returns:

Type	Description
NoiseCovLO_Toeplitz01	The noise covariance operator \(N\)

Source code in brahmap/core/noise_ops_toeplitz.py

def get_inverse(self) -> "NoiseCovLO_Toeplitz01":  # type: ignore
    """Returns the inverse of this operator, which is the original
    noise covariance operator.

    Returns
    -------
    NoiseCovLO_Toeplitz01
        The noise covariance operator $N$
    """
    return self.__toeplitz_op