matrix

matrix ¶

Convolution matrix construction and inversion.

This module provides functions for building convolution matrices and computing their regularized pseudo-inverses for deconvolution operations.

References

- dcmri (https://github.com/dcmri/dcmri) utils.convmat(), utils.invconvmat()
- Wu et al. (2003). Tracer arrival timing-insensitive technique.
  Magn Reson Med. 50:164-174.

convmat ¶

convmat(h, t, *, order=1)

Construct a convolution matrix from an impulse response.

Creates a lower triangular matrix A such that the convolution f * h can be computed as A @ h when f is known, or the deconvolution can be computed by solving A @ x = y for x.

PARAMETER	DESCRIPTION
`h`	Impulse response function (e.g., AIF). Shape: (n_times,). TYPE: `ndarray`
`t`	Time points in seconds. Shape: (n_times,). Must be monotonically increasing. TYPE: `ndarray`
`order`	Integration order: - 1: First-order (trapezoidal) integration - 2: Second-order (Simpson's rule) integration TYPE: `(1, 2)` DEFAULT: `1`

RETURNS	DESCRIPTION
`ndarray`	Convolution matrix. Shape: (n_times, n_times). Lower triangular Toeplitz-like structure.

Notes

For uniform time grids, the matrix is Toeplitz (constant diagonals). For non-uniform grids, the matrix accounts for varying time steps.

The matrix is constructed such that: (f * h)[i] = sum_j A[i,j] * f[j]

Using first-order (trapezoidal) integration: A[i,j] = dt[j] * (h[i-j] + h[i-j-1]) / 2 for j < i A[i,i] = 0

Examples:

>>> import numpy as np
>>> from osipy.common.convolution import convmat
>>> t = np.linspace(0, 10, 11)
>>> aif = np.exp(-t / 3)
>>> A = convmat(aif, t)
>>> # Convolution: result = A @ signal

References

.. [1] dcmri library: https://github.com/dcmri/dcmri

invconvmat ¶

invconvmat(A, *, method='tsvd', tol=0.1, lambda_reg=None)

Compute regularized pseudo-inverse of a convolution matrix.

Uses SVD decomposition with regularization to compute a stable inverse of the convolution matrix for deconvolution operations.

PARAMETER	DESCRIPTION
`A`	Convolution matrix. Shape: (n, n) or (m, n). TYPE: `ndarray`
`method`	Regularization method: - "tsvd": Truncated SVD. Singular values below tol * max(s) are set to zero. - "tikhonov": Tikhonov regularization. Singular values are damped as s / (s^2 + lambda^2). TYPE: `('tsvd', 'tikhonov')` DEFAULT: `"tsvd"`
`tol`	For TSVD: Relative tolerance for singular value truncation. Singular values s < tol * max(s) are set to zero. For Tikhonov: Used to compute lambda if not provided. TYPE: `float` DEFAULT: `0.1`
`lambda_reg`	Tikhonov regularization parameter. If not provided, computed as lambda = tol * max(singular_values). TYPE: `float` DEFAULT: `None`

RETURNS	DESCRIPTION
`ndarray`	Regularized pseudo-inverse. Shape: (n, m) where A is (m, n).

Notes

TSVD (Truncated Singular Value Decomposition): A_inv = V @ diag(1/s_truncated) @ U.T where s_truncated sets small singular values to 0.

Tikhonov regularization: A_inv = V @ diag(s / (s^2 + lambda^2)) @ U.T This damps rather than truncates small singular values.

The choice of regularization affects noise amplification in deconvolution. TSVD is simpler but can introduce ringing artifacts. Tikhonov provides smoother results but may over-smooth.

Examples:

>>> import numpy as np
>>> from osipy.common.convolution import convmat, invconvmat
>>> t = np.linspace(0, 10, 51)
>>> aif = np.exp(-t / 3)
>>> A = convmat(aif, t)
>>> A_inv = invconvmat(A, method="tsvd", tol=0.1)
>>> # Deconvolution: irf = A_inv @ signal

References

.. [1] dcmri library: https://github.com/dcmri/dcmri .. [2] Hansen PC (1998). Rank-Deficient and Discrete Ill-Posed Problems.

circulant_convmat ¶

circulant_convmat(h, t)

Construct a block-circulant convolution matrix.

Creates a circulant matrix for delay-insensitive deconvolution methods (cSVD/oSVD). The circulant structure allows the AIF to wrap around, making the deconvolution insensitive to tracer arrival time differences.

PARAMETER	DESCRIPTION
`h`	Impulse response function (e.g., AIF). Shape: (n_times,). TYPE: `ndarray`
`t`	Time points in seconds. Shape: (n_times,). TYPE: `ndarray`

RETURNS	DESCRIPTION
`ndarray`	Block-circulant convolution matrix. Shape: (n_times, n_times).

Notes

The circulant matrix has the form: C[i,j] = h[(i-j) mod n]

This is equivalent to circular convolution rather than linear convolution. The benefit is that the deconvolution becomes insensitive to shifts between the AIF and tissue curve.

References

.. [1] Wu et al. (2003). Tracer arrival timing-insensitive technique. Magn Reson Med. 50:164-174.