CGLS

Krylov.cglsFunction
(x, stats) = cgls(A, b::AbstractVector{FC};
                  M=I, ldiv::Bool=false, radius::T=zero(T),
                  λ::T=zero(T), atol::T=√eps(T), rtol::T=√eps(T),
                  itmax::Int=0, timemax::Float64=Inf,
                  verbose::Int=0, history::Bool=false,
                  callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

Solve the regularized linear least-squares problem

minimize ‖b - Ax‖₂² + λ‖x‖₂²

of size m × n using the Conjugate Gradient (CG) method, where λ ≥ 0 is a regularization parameter. This method is equivalent to applying CG to the normal equations

(AᴴA + λI) x = Aᴴb

but is more stable.

CGLS produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᴴr‖₂. It is formally equivalent to LSQR, though can be slightly less accurate, but simpler to implement.

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m.

Keyword arguments

  • M: linear operator that models a Hermitian positive-definite matrix of size n used for preconditioning;
  • ldiv: define whether the preconditioner uses ldiv! or mul!;
  • radius: add the trust-region constraint ‖x‖ ≤ radius if radius > 0. Useful to compute a step in a trust-region method for optimization;
  • λ: regularization parameter;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;
  • stats: statistics collected on the run in a SimpleStats structure.

References

source
Krylov.cgls!Function
solver = cgls!(solver::CglsSolver, A, b; kwargs...)

where kwargs are keyword arguments of cgls.

See CglsSolver for more details about the solver.

source

CRLS

Krylov.crlsFunction
(x, stats) = crls(A, b::AbstractVector{FC};
                  M=I, ldiv::Bool=false, radius::T=zero(T),
                  λ::T=zero(T), atol::T=√eps(T), rtol::T=√eps(T),
                  itmax::Int=0, timemax::Float64=Inf,
                  verbose::Int=0, history::Bool=false,
                  callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

Solve the linear least-squares problem

minimize ‖b - Ax‖₂² + λ‖x‖₂²

of size m × n using the Conjugate Residuals (CR) method. This method is equivalent to applying MINRES to the normal equations

(AᴴA + λI) x = Aᴴb.

This implementation recurs the residual r := b - Ax.

CRLS produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr‖₂. It is formally equivalent to LSMR, though can be substantially less accurate, but simpler to implement.

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m.

Keyword arguments

  • M: linear operator that models a Hermitian positive-definite matrix of size n used for preconditioning;
  • ldiv: define whether the preconditioner uses ldiv! or mul!;
  • radius: add the trust-region constraint ‖x‖ ≤ radius if radius > 0. Useful to compute a step in a trust-region method for optimization;
  • λ: regularization parameter;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;
  • stats: statistics collected on the run in a SimpleStats structure.

Reference

  • D. C.-L. Fong, Minimum-Residual Methods for Sparse, Least-Squares using Golubg-Kahan Bidiagonalization, Ph.D. Thesis, Stanford University, 2011.
source
Krylov.crls!Function
solver = crls!(solver::CrlsSolver, A, b; kwargs...)

where kwargs are keyword arguments of crls.

See CrlsSolver for more details about the solver.

source

LSLQ

Krylov.lslqFunction
(x, stats) = lslq(A, b::AbstractVector{FC};
                  M=I, N=I, ldiv::Bool=false,
                  window::Int=5, transfer_to_lsqr::Bool=false,
                  sqd::Bool=false, λ::T=zero(T),
                  σ::T=zero(T), etol::T=√eps(T),
                  utol::T=√eps(T), btol::T=√eps(T),
                  conlim::T=1/√eps(T), atol::T=√eps(T),
                  rtol::T=√eps(T), itmax::Int=0,
                  timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                  callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

Solve the regularized linear least-squares problem

minimize ‖b - Ax‖₂² + λ²‖x‖₂²

of size m × n using the LSLQ method, where λ ≥ 0 is a regularization parameter. LSLQ is formally equivalent to applying SYMMLQ to the normal equations

(AᴴA + λ²I) x = Aᴴb

but is more stable.

If λ > 0, we solve the symmetric and quasi-definite system

[ E      A ] [ r ]   [ b ]
[ Aᴴ  -λ²F ] [ x ] = [ 0 ],

where E and F are symmetric and positive definite. Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators. If sqd=true, λ is set to the common value 1.

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.

For a symmetric and positive definite matrix K, the K-norm of a vector x is ‖x‖²_K = xᴴKx. LSLQ is then equivalent to applying SYMMLQ to (AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b with r = E⁻¹(b - Ax).

If λ = 0, we solve the symmetric and indefinite system

[ E    A ] [ r ]   [ b ]
[ Aᴴ   0 ] [ x ] = [ 0 ].

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹.

In this case, N can still be specified and indicates the weighted norm in which x and Aᴴr should be measured. r can be recovered by computing E⁻¹(b - Ax).

Main features

  • the solution estimate is updated along orthogonal directions
  • the norm of the solution estimate ‖xᴸₖ‖₂ is increasing
  • the error ‖eₖ‖₂ := ‖xᴸₖ - x*‖₂ is decreasing
  • it is possible to transition cheaply from the LSLQ iterate to the LSQR iterate if there is an advantage (there always is in terms of error)
  • if A is rank deficient, identify the minimum least-squares solution

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m.

Keyword arguments

  • M: linear operator that models a Hermitian positive-definite matrix of size m used for centered preconditioning of the augmented system;
  • N: linear operator that models a Hermitian positive-definite matrix of size n used for centered preconditioning of the augmented system;
  • ldiv: define whether the preconditioners use ldiv! or mul!;
  • window: number of iterations used to accumulate a lower bound on the error;
  • transfer_to_lsqr: transfer from the LSLQ point to the LSQR point, when it exists. The transfer is based on the residual norm;
  • sqd: if true, set λ=1 for Hermitian quasi-definite systems;
  • λ: regularization parameter;
  • σ: strict lower bound on the smallest positive singular value σₘᵢₙ such as σ = (1-10⁻⁷)σₘᵢₙ;
  • etol: stopping tolerance based on the lower bound on the error;
  • utol: stopping tolerance based on the upper bound on the error;
  • btol: stopping tolerance used to detect zero-residual problems;
  • conlim: limit on the estimated condition number of A beyond which the solution will be abandoned;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;

  • stats: statistics collected on the run in a LSLQStats structure.

  • stats.err_lbnds is a vector of lower bounds on the LQ error–-the vector is empty if window is set to zero

  • stats.err_ubnds_lq is a vector of upper bounds on the LQ error–-the vector is empty if σ == 0 is left at zero

  • stats.err_ubnds_cg is a vector of upper bounds on the CG error–-the vector is empty if σ == 0 is left at zero

  • stats.error_with_bnd is a boolean indicating whether there was an error in the upper bounds computation (cancellation errors, too large σ ...)

Stopping conditions

The iterations stop as soon as one of the following conditions holds true:

  • the optimality residual is sufficiently small (stats.status = "found approximate minimum least-squares solution") in the sense that either
    • ‖Aᴴr‖ / (‖A‖ ‖r‖) ≤ atol, or
    • 1 + ‖Aᴴr‖ / (‖A‖ ‖r‖) ≤ 1
  • an approximate zero-residual solution has been found (stats.status = "found approximate zero-residual solution") in the sense that either
    • ‖r‖ / ‖b‖ ≤ btol + atol ‖A‖ * ‖xᴸ‖ / ‖b‖, or
    • 1 + ‖r‖ / ‖b‖ ≤ 1
  • the estimated condition number of A is too large in the sense that either
    • 1/cond(A) ≤ 1/conlim (stats.status = "condition number exceeds tolerance"), or
    • 1 + 1/cond(A) ≤ 1 (stats.status = "condition number seems too large for this machine")
  • the lower bound on the LQ forward error is less than etol * ‖xᴸ‖
  • the upper bound on the CG forward error is less than utol * ‖xᶜ‖

References

source
Krylov.lslq!Function
solver = lslq!(solver::LslqSolver, A, b; kwargs...)

where kwargs are keyword arguments of lslq.

See LslqSolver for more details about the solver.

source

LSQR

Krylov.lsqrFunction
(x, stats) = lsqr(A, b::AbstractVector{FC};
                  M=I, N=I, ldiv::Bool=false,
                  window::Int=5, sqd::Bool=false, λ::T=zero(T),
                  radius::T=zero(T), etol::T=√eps(T),
                  axtol::T=√eps(T), btol::T=√eps(T),
                  conlim::T=1/√eps(T), atol::T=zero(T),
                  rtol::T=zero(T), itmax::Int=0,
                  timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                  callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

Solve the regularized linear least-squares problem

minimize ‖b - Ax‖₂² + λ²‖x‖₂²

of size m × n using the LSQR method, where λ ≥ 0 is a regularization parameter. LSQR is formally equivalent to applying CG to the normal equations

(AᴴA + λ²I) x = Aᴴb

(and therefore to CGLS) but is more stable.

LSQR produces monotonic residuals ‖r‖₂ but not optimality residuals ‖Aᴴr‖₂. It is formally equivalent to CGLS, though can be slightly more accurate.

If λ > 0, LSQR solves the symmetric and quasi-definite system

[ E      A ] [ r ]   [ b ]
[ Aᴴ  -λ²F ] [ x ] = [ 0 ],

where E and F are symmetric and positive definite. Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators. If sqd=true, λ is set to the common value 1.

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.

For a symmetric and positive definite matrix K, the K-norm of a vector x is ‖x‖²_K = xᴴKx. LSQR is then equivalent to applying CG to (AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b with r = E⁻¹(b - Ax).

If λ = 0, we solve the symmetric and indefinite system

[ E    A ] [ r ]   [ b ]
[ Aᴴ   0 ] [ x ] = [ 0 ].

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹.

In this case, N can still be specified and indicates the weighted norm in which x and Aᴴr should be measured. r can be recovered by computing E⁻¹(b - Ax).

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m.

Keyword arguments

  • M: linear operator that models a Hermitian positive-definite matrix of size m used for centered preconditioning of the augmented system;
  • N: linear operator that models a Hermitian positive-definite matrix of size n used for centered preconditioning of the augmented system;
  • ldiv: define whether the preconditioners use ldiv! or mul!;
  • window: number of iterations used to accumulate a lower bound on the error;
  • sqd: if true, set λ=1 for Hermitian quasi-definite systems;
  • λ: regularization parameter;
  • radius: add the trust-region constraint ‖x‖ ≤ radius if radius > 0. Useful to compute a step in a trust-region method for optimization;
  • etol: stopping tolerance based on the lower bound on the error;
  • axtol: tolerance on the backward error;
  • btol: stopping tolerance used to detect zero-residual problems;
  • conlim: limit on the estimated condition number of A beyond which the solution will be abandoned;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;
  • stats: statistics collected on the run in a SimpleStats structure.

Reference

source
Krylov.lsqr!Function
solver = lsqr!(solver::LsqrSolver, A, b; kwargs...)

where kwargs are keyword arguments of lsqr.

See LsqrSolver for more details about the solver.

source

LSMR

Krylov.lsmrFunction
(x, stats) = lsmr(A, b::AbstractVector{FC};
                  M=I, N=I, ldiv::Bool=false,
                  window::Int=5, sqd::Bool=false, λ::T=zero(T),
                  radius::T=zero(T), etol::T=√eps(T),
                  axtol::T=√eps(T), btol::T=√eps(T),
                  conlim::T=1/√eps(T), atol::T=zero(T),
                  rtol::T=zero(T), itmax::Int=0,
                  timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                  callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

Solve the regularized linear least-squares problem

minimize ‖b - Ax‖₂² + λ²‖x‖₂²

of size m × n using the LSMR method, where λ ≥ 0 is a regularization parameter. LSMR is formally equivalent to applying MINRES to the normal equations

(AᴴA + λ²I) x = Aᴴb

(and therefore to CRLS) but is more stable.

LSMR produces monotonic residuals ‖r‖₂ and optimality residuals ‖Aᴴr‖₂. It is formally equivalent to CRLS, though can be substantially more accurate.

LSMR can be also used to find a null vector of a singular matrix A by solving the problem min ‖Aᴴx - b‖ with any nonzero vector b. At a minimizer, the residual vector r = b - Aᴴx will satisfy Ar = 0.

If λ > 0, we solve the symmetric and quasi-definite system

[ E      A ] [ r ]   [ b ]
[ Aᴴ  -λ²F ] [ x ] = [ 0 ],

where E and F are symmetric and positive definite. Preconditioners M = E⁻¹ ≻ 0 and N = F⁻¹ ≻ 0 may be provided in the form of linear operators. If sqd=true, λ is set to the common value 1.

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹ + λ²‖x‖²_F.

For a symmetric and positive definite matrix K, the K-norm of a vector x is ‖x‖²_K = xᴴKx. LSMR is then equivalent to applying MINRES to (AᴴE⁻¹A + λ²F)x = AᴴE⁻¹b with r = E⁻¹(b - Ax).

If λ = 0, we solve the symmetric and indefinite system

[ E    A ] [ r ]   [ b ]
[ Aᴴ   0 ] [ x ] = [ 0 ].

The system above represents the optimality conditions of

minimize ‖b - Ax‖²_E⁻¹.

In this case, N can still be specified and indicates the weighted norm in which x and Aᴴr should be measured. r can be recovered by computing E⁻¹(b - Ax).

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m.

Keyword arguments

  • M: linear operator that models a Hermitian positive-definite matrix of size m used for centered preconditioning of the augmented system;
  • N: linear operator that models a Hermitian positive-definite matrix of size n used for centered preconditioning of the augmented system;
  • ldiv: define whether the preconditioners use ldiv! or mul!;
  • window: number of iterations used to accumulate a lower bound on the error;
  • sqd: if true, set λ=1 for Hermitian quasi-definite systems;
  • λ: regularization parameter;
  • radius: add the trust-region constraint ‖x‖ ≤ radius if radius > 0. Useful to compute a step in a trust-region method for optimization;
  • etol: stopping tolerance based on the lower bound on the error;
  • axtol: tolerance on the backward error;
  • btol: stopping tolerance used to detect zero-residual problems;
  • conlim: limit on the estimated condition number of A beyond which the solution will be abandoned;
  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;
  • stats: statistics collected on the run in a LsmrStats structure.

Reference

source
Krylov.lsmr!Function
solver = lsmr!(solver::LsmrSolver, A, b; kwargs...)

where kwargs are keyword arguments of lsmr.

See LsmrSolver for more details about the solver.

source

USYMQR

Krylov.usymqrFunction
(x, stats) = usymqr(A, b::AbstractVector{FC}, c::AbstractVector{FC};
                    atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
                    timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                    callback=solver->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(x, stats) = usymqr(A, b, c, x0::AbstractVector; kwargs...)

USYMQR can be warm-started from an initial guess x0 where kwargs are the same keyword arguments as above.

USYMQR solves the linear least-squares problem min ‖b - Ax‖² of size m × n. USYMQR solves Ax = b if it is consistent.

USYMQR is based on the orthogonal tridiagonalization process and requires two initial nonzero vectors b and c. The vector c is only used to initialize the process and a default value can be b or Aᴴb depending on the shape of A. The residual norm ‖b - Ax‖ monotonously decreases in USYMQR. It's considered as a generalization of MINRES.

It can also be applied to under-determined and over-determined problems. USYMQR finds the minimum-norm solution if problems are inconsistent.

Input arguments

  • A: a linear operator that models a matrix of dimension m × n;
  • b: a vector of length m;
  • c: a vector of length n.

Optional argument

  • x0: a vector of length n that represents an initial guess of the solution x.

Keyword arguments

  • atol: absolute stopping tolerance based on the residual norm;
  • rtol: relative stopping tolerance based on the residual norm;
  • itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
  • timemax: the time limit in seconds;
  • verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
  • history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
  • callback: function or functor called as callback(solver) that returns true if the Krylov method should terminate, and false otherwise;
  • iostream: stream to which output is logged.

Output arguments

  • x: a dense vector of length n;
  • stats: statistics collected on the run in a SimpleStats structure.

References

source
Krylov.usymqr!Function
solver = usymqr!(solver::UsymqrSolver, A, b, c; kwargs...)
solver = usymqr!(solver::UsymqrSolver, A, b, c, x0; kwargs...)

where kwargs are keyword arguments of usymqr.

See UsymqrSolver for more details about the solver.

source