lsqr#

mrinufft.extras.lsqr(operator: FourierOperatorBase, kspace_data: NDArray, damp: float = 0.0, atol: float = 1e-06, btol: float = 1e-06, conlim: float = 100000000.0, max_iter: int = 100, x0: NDArray | None = None, x_init: NDArray | None = None, callback: Callable | None = None, progressbar: bool | tqdm = True)[source]#

Solve a regularized linear least-squares problem with the LSQR algorithm [1]_.

This implementation follows the numerically stable formulation described in [2]_ and is adapted from the SciPy implementation [3]_.

Solves problems of the form

\[\arg\min \|A x - b\|_2^2 + \gamma^2 \|x - x0\|_2^2\]

Stop iterating if:

  • numerical convergence is reached: \(\|Ax-b\| \leq atol \|A\| * \|x\| + btol * \|b\|\)

  • estimation of the conditioning of the problem diverge: cond(A)>=conlim

  • Maximum number of iteration reached.

Parameters:

  • nufft (FourierOperatorBase) – The NUFFT operator representing the forward model.

  • kspace_data (NDArray) – The right-hand side vector. Shape is typically (n_batchs, n_coils, n_samples).

  • damp (float, optional) – Damping (regularization) parameter. Default is 0.0 (no regularization).

  • x0 (NDArray or None, optional) – Damping vector. If None, uses zero. Shape is typically (n_batchs, n_coils or 1, *nufft.shape).

  • x_init (NDArray or None, optional) – Initial guess vector. If ommitted, default to x0. Must have same shape as x0.

  • callback (Callable, optional) – If provided, a callback function will be called at the end of each iteration with the current estimate. It should have the following signature callback(operator, kspace_data, damp, x0)

  • max_iter (int, optional) – Maximum number of iterations. Default is 100.

  • progressbar (bool, optional) – If True (default) display a progress bar to track iterations.

atolfloat, optional

Stopping tolerance on the absolute error. Default is 1e-6.

btolfloat, optional

Stopping tolerance on the relative error. Default is 1e-6.

conlimfloat, optional

Limit on condition number. Iteration stops if condition exceeds this value. Default is 1e8.

Returns:

Solution vector with shape (n_batchs, n_coils or 1, *nufft.shape), dtype and device matching input.

Return type:

NDArray

References

Paige, C. C., & Saunders, M. A. (1982). LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software, 8(1), 43-71.

S.-C. Choi, “Iterative Methods for Singular Linear Equations and Least-Squares Problems”, Dissertation, http://www.stanford.edu/group/SOL/dissertations/sou-cheng-choi-thesis.pdf

scipy/scipy

Note

This function uses numpy for all CPU arrays, and cupy for all on-gpu array. It will convert all its array argument to the respective array library. The outputs will be converted back to the original array module and device.

See also

This function is part of the optimizer registry, get_optimizer(). You can find other registered functions in this registry below:

Key

Function

"lsqr"

lsqr(operator, kspace_data, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, max_iter=100, x0=None, x_init=None, callback=None, progressbar=True)

"lsmr"

lsmr(operator, kspace_data, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, max_iter=100, x0=None, x_init=None, callback=None, progressbar=True)

"cg"

cg(operator, kspace_data, damp=0.0, x0=None, x_init=None, max_iter=10, tol=0.0001, progressbar=True, callback=None)