LSMR

template<Number T, DataLocation L = Host>
class LSMR

Iterative solver for least-square problems using the LSMR algorithm [1]. It solves the system of linear equations \( \mathbf{A}\vec{x} = \vec{b} \). If the system is inconsistent, it solves the least-square problem

\[ \min_{\vec{x}} \| \mathbf{A}\vec{x} - \vec{b} \|_2 \]
where A is a real rectangular matrix of dimension m-by-n and b is a real vector of length m. The matrix A may be dense or sparse (usually sparse).

To improve the convergence rate of the solver, you can specify a right preconditioner \( \mathbf{P}^{-1} \), such that now LSMR solves the following least square problem

\[ \min_{\vec{x}} \| (\mathbf{AP})^\dagger \; \mathbf{P} \vec{x} - \vec{b} \|_2 \]

More specifically, the Preconditioner object is a function that computes \( \vec{y} = \mathbf{P}^{-1} \mathbf{P}^{-T} \vec{x}\) for an input vector x and output vector y.

The code was adapted from the SciPy package (https://github.com/scipy/scipy.git).

References:

[1] D. C.-L. Fong and M. Saunders, ‘LSMR: An Iterative Algorithm for Sparse Least-Squares Problems’, SIAM J. Sci. Comput., vol. 33, no. 5, pp. 2950–2971, Jan. 2011, doi: 10.1137/10079687X.

Template Parameters:
  • T – Data type used by the solver. Must be either float or double.

  • L – Location of the data.

Public Functions

LSMR(LSMROptions options = LSMROptions())

Creates a new LSMR solver.

Parameters:

options[in] solver parameters. See LSMROptions for more information.

virtual ~LSMR() = default

Default destructor.

inline LSMRMatResult operator()(MatrixMap<T, L> A, VectorArgs b, std::optional<VectorArgs> x0 = std::nullopt, std::optional<Preconditioner<T, L>> precond = std::nullopt)

Solves the least square problem Ax = b.

Parameters:
  • A[in] coefficient matrix

  • b[in] right-hand side of the equation

  • x0[in] initial guess for the solution (optional)

  • precond[in] a Preconditioner to the system (optional). Suppose that \( \mathbf{L}^{-1}\) is the preconditioner of system, then precond must apply \( \mathbf{L}^{-1} \mathbf{L}^{-T}\) over a vector x.

Returns:

An abstract object representing the solution of the linear system.

LSMRReport report()
Returns:

Some relevant information about the last run. See LSMRReport for more information.

void set_params(LSMROptions options)

Set the parameters of the solver.

Parameters:

options[in] solver parameters. See LSMROptions for more information.

struct LSMROptions

Parameters of the LSMR solver.

Public Members

double damping_factor = 0

Damping factor \( \lambda \) for regularized least-squares problems. If \( \lambda \) is not 0, LSMR solves the regularized least-squares problem::

\[\begin{split} \min_x \| \begin{bmatrix} \vec{b} \\ 0 \end{bmatrix} - \begin{bmatrix} \mathbf{A}\vec{x} \\ \lambda \mathbf{I} \end{bmatrix} \|_2 \end{split}\]

double atol = 1E-6

Stopping criteria. If the system Ax = b seems to be consistent, LSMR terminates when \( \| \vec{r}\|_2 \leq atol * \| \mathbf{A}\vec{x} \|_2 + btol * \| \vec{b}\|\). Otherwise, it terminates when \( \left \| \mathbf{A}^H \vec{r}\|_2 \leq atol * \| \mathbf{A}\vec{r} \right \|_2 \).

double cond_limit = 1E8

LSMR terminates when the estimate for cond(A) exceeds cond_limit.

index_t max_iter = kAuto

LSMR terminates when the number of iterations reaches max_iter. The default is max_iter = min(m, n).

int verbose = 0

Set the verbosity level. 0 - No report. 1 - Report only the residual

\[ \| \mathbf{A}\vec{x} - \vec{b} \|_2 \]
2 - Full iteration log.

struct LSMRReport

Report of the last run of the LSMR solver.

Public Members

std::string stop_condition

Stop condition.

index_t num_iter

Number of iterations.

double cond_A

Estimated condition number of the matrix A

double norm_residual

Norm of the residual \( \| \vec{b} - \mathbf{A}\vec{x}\|_2 \).