lsmr#

mrinufft.extras.lsmr(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 LSMR algorithm [1]_.

See the reference LSMR software [2]_ for further details.

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

Notes

  • LSMR is generally more stable than LSQR, notably in term of image residual norm \(\|A^H(Ax-b)\|\), and is similar to the MINRES algorithm for least squares problems.

  • It usually converges faster than LSQR and can stop in fewer iterations.

References

D. C.-L. Fong and M. A. Saunders, “LSMR: An iterative algorithm for sparse least-squares problems”, SIAM J. Sci. Comput., vol. 33, pp. 2950-2971, 2011.

LSMR Software, https://web.stanford.edu/group/SOL/software/lsmr/

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)