(x, stats) = cgne(A, b::AbstractVector{FC};
                  N=I, ldiv::Bool=false,
                  λ::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 consistent linear system

Ax + √λs = b

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 of the second kind

(AAᴴ + λI) y = b

but is more stable. When λ = 0, this method solves the minimum-norm problem

min ‖x‖₂ s.t. Ax = b.

When λ > 0, it solves the problem

min ‖(x,s)‖₂  s.t. Ax + √λs = b.

CGNE produces monotonic errors ‖x-x*‖₂ but not residuals ‖r‖₂. It is formally equivalent to CRAIG, though can be slightly less accurate, but simpler to implement. Only the x-part of the solution is returned.

Input arguments

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

Keyword arguments

• N: linear operator that models a Hermitian positive-definite matrix of size n used for preconditioning;
• ldiv: define whether the preconditioner uses ldiv! or mul!;
• λ: 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

• J. E. Craig, The N-step iteration procedures, Journal of Mathematics and Physics, 34(1), pp. 64–73, 1955.
• J. E. Craig, Iterations Procedures for Simultaneous Equations, Ph.D. Thesis, Department of Electrical Engineering, MIT, 1954.

solver = cgne!(solver::CgneSolver, A, b; kwargs...)

where kwargs are keyword arguments of cgne.

See CgneSolver for more details about the solver.

(x, stats) = crmr(A, b::AbstractVector{FC};
                  N=I, ldiv::Bool=false,
                  λ::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 consistent linear system

Ax + √λs = b

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

(AAᴴ + λI) y = b

but is more stable. When λ = 0, this method solves the minimum-norm problem

min ‖x‖₂  s.t.  x ∈ argmin ‖Ax - b‖₂.

When λ > 0, this method solves the problem

min ‖(x,s)‖₂  s.t. Ax + √λs = b.

CRMR produces monotonic residuals ‖r‖₂. It is formally equivalent to CRAIG-MR, though can be slightly less accurate, but simpler to implement. Only the x-part of the solution is returned.

Input arguments

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

Keyword arguments

• N: linear operator that models a Hermitian positive-definite matrix of size n used for preconditioning;
• ldiv: define whether the preconditioner uses ldiv! or mul!;
• λ: 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

• D. Orban and M. Arioli, Iterative Solution of Symmetric Quasi-Definite Linear Systems, Volume 3 of Spotlights. SIAM, Philadelphia, PA, 2017.
• D. Orban, The Projected Golub-Kahan Process for Constrained Linear Least-Squares Problems. Cahier du GERAD G-2014-15, 2014.

solver = crmr!(solver::CrmrSolver, A, b; kwargs...)

where kwargs are keyword arguments of crmr.

See CrmrSolver for more details about the solver.

(x, y, stats) = lnlq(A, b::AbstractVector{FC};
                     M=I, N=I, ldiv::Bool=false,
                     transfer_to_craig::Bool=true,
                     sqd::Bool=false, λ::T=zero(T),
                     σ::T=zero(T), utolx::T=√eps(T),
                     utoly::T=√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}.

Find the least-norm solution of the consistent linear system

Ax + λ²y = b

of size m × n using the LNLQ method, where λ ≥ 0 is a regularization parameter.

For a system in the form Ax = b, LNLQ method is equivalent to applying SYMMLQ to AAᴴy = b and recovering x = Aᴴy but is more stable. Note that y are the Lagrange multipliers of the least-norm problem

minimize ‖x‖  s.t.  Ax = b.

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

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

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

min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.

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

If λ = 0, LNLQ solves the symmetric and indefinite system

[ -F   Aᴴ ] [ x ]   [ 0 ]
[  A   0  ] [ y ] = [ b ].

The system above represents the optimality conditions of

minimize ‖x‖²_F  s.t.  Ax = b.

In this case, M can still be specified and indicates the weighted norm in which residuals are measured.

In this implementation, both the x and y-parts of the solution are returned.

utolx and utoly are tolerances on the upper bound of the distance to the solution ‖x-x‖ and ‖y-y‖, respectively. The bound is valid if λ>0 or σ>0 where σ should be strictly smaller than the smallest positive singular value. For instance σ:=(1-1e-7)σₘᵢₙ .

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!;
• transfer_to_craig: transfer from the LNLQ point to the CRAIG 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⁻⁷)σₘᵢₙ;
• utolx: tolerance on the upper bound on the distance to the solution ‖x-x*‖;
• utoly: tolerance on the upper bound on the distance to the solution ‖y-y*‖;
• 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;
• y: a dense vector of length m;
• stats: statistics collected on the run in a LNLQStats structure.

Reference

• R. Estrin, D. Orban, M.A. Saunders, LNLQ: An Iterative Method for Least-Norm Problems with an Error Minimization Property, SIAM Journal on Matrix Analysis and Applications, 40(3), pp. 1102–1124, 2019.

solver = lnlq!(solver::LnlqSolver, A, b; kwargs...)

where kwargs are keyword arguments of lnlq.

See LnlqSolver for more details about the solver.

(x, y, stats) = craig(A, b::AbstractVector{FC};
                      M=I, N=I, ldiv::Bool=false,
                      transfer_to_lsqr::Bool=false, sqd::Bool=false,
                      λ::T=zero(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}.

Find the least-norm solution of the consistent linear system

Ax + λ²y = b

of size m × n using the Golub-Kahan implementation of Craig's method, where λ ≥ 0 is a regularization parameter. This method is equivalent to CGNE but is more stable.

For a system in the form Ax = b, Craig's method is equivalent to applying CG to AAᴴy = b and recovering x = Aᴴy. Note that y are the Lagrange multipliers of the least-norm problem

minimize ‖x‖  s.t.  Ax = b.

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

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

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

min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.

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

If λ = 0, CRAIG solves the symmetric and indefinite system

[ -F   Aᴴ ] [ x ]   [ 0 ]
[  A   0  ] [ y ] = [ b ].

The system above represents the optimality conditions of

minimize ‖x‖²_F  s.t.  Ax = b.

In this case, M can still be specified and indicates the weighted norm in which residuals are measured.

In this implementation, both the x and y-parts of the solution are returned.

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!;
• 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;
• 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;
• y: a dense vector of length m;
• stats: statistics collected on the run in a SimpleStats structure.

References

• J. E. Craig, The N-step iteration procedures, Journal of Mathematics and Physics, 34(1-4), pp. 64–73, 1955.
• C. C. Paige and M. A. Saunders, LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares, ACM Transactions on Mathematical Software, 8(1), pp. 43–71, 1982.
• M. A. Saunders, Solutions of Sparse Rectangular Systems Using LSQR and CRAIG, BIT Numerical Mathematics, 35(4), pp. 588–604, 1995.

solver = craig!(solver::CraigSolver, A, b; kwargs...)

where kwargs are keyword arguments of craig.

See CraigSolver for more details about the solver.

(x, y, stats) = craigmr(A, b::AbstractVector{FC};
                        M=I, N=I, ldiv::Bool=false,
                        sqd::Bool=false, λ::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 consistent linear system

Ax + λ²y = b

of size m × n using the CRAIGMR method, where λ ≥ 0 is a regularization parameter. This method is equivalent to applying the Conjugate Residuals method to the normal equations of the second kind

(AAᴴ + λ²I) y = b

but is more stable. When λ = 0, this method solves the minimum-norm problem

min ‖x‖  s.t.  x ∈ argmin ‖Ax - b‖.

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

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

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

min ‖x‖²_F + λ²‖y‖²_E  s.t.  Ax + λ²Ey = b.

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

If λ = 0, CRAIGMR solves the symmetric and indefinite system

[ -F   Aᴴ ] [ x ]   [ 0 ]
[  A   0  ] [ y ] = [ b ].

The system above represents the optimality conditions of

min ‖x‖²_F  s.t.  Ax = b.

In this case, M can still be specified and indicates the weighted norm in which residuals are measured.

CRAIGMR produces monotonic residuals ‖r‖₂. It is formally equivalent to CRMR, though can be slightly more accurate, and intricate to implement. Both the x- and y-parts of the solution are returned.

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!;
• sqd: if true, set λ=1 for Hermitian quasi-definite systems;
• λ: 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;
• y: a dense vector of length m;
• stats: statistics collected on the run in a SimpleStats structure.

References

• D. Orban and M. Arioli. Iterative Solution of Symmetric Quasi-Definite Linear Systems, Volume 3 of Spotlights. SIAM, Philadelphia, PA, 2017.
• D. Orban, The Projected Golub-Kahan Process for Constrained, Linear Least-Squares Problems. Cahier du GERAD G-2014-15, 2014.

solver = craigmr!(solver::CraigmrSolver, A, b; kwargs...)

where kwargs are keyword arguments of craigmr.

See CraigmrSolver for more details about the solver.

(x, stats) = usymlq(A, b::AbstractVector{FC}, c::AbstractVector{FC};
                    transfer_to_usymcg::Bool=true, 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) = usymlq(A, b, c, x0::AbstractVector; kwargs...)

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

USYMLQ determines the least-norm solution of the consistent linear system Ax = b of size m × n.

USYMLQ 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 error norm ‖x - x*‖ monotonously decreases in USYMLQ. It's considered as a generalization of SYMMLQ.

It can also be applied to under-determined and over-determined problems. In all cases, problems must be consistent.

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

• transfer_to_usymcg: transfer from the USYMLQ point to the USYMCG point, when it exists. The transfer is based on the residual norm;
• 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

• M. A. Saunders, H. D. Simon, and E. L. Yip, Two Conjugate-Gradient-Type Methods for Unsymmetric Linear Equations, SIAM Journal on Numerical Analysis, 25(4), pp. 927–940, 1988.
• A. Buttari, D. Orban, D. Ruiz and D. Titley-Peloquin, A tridiagonalization method for symmetric saddle-point and quasi-definite systems, SIAM Journal on Scientific Computing, 41(5), pp. 409–432, 2019.
• A. Montoison and D. Orban, BiLQ: An Iterative Method for Nonsymmetric Linear Systems with a Quasi-Minimum Error Property, SIAM Journal on Matrix Analysis and Applications, 41(3), pp. 1145–1166, 2020.

solver = usymlq!(solver::UsymlqSolver, A, b, c; kwargs...)
solver = usymlq!(solver::UsymlqSolver, A, b, c, x0; kwargs...)

where kwargs are keyword arguments of usymlq.

See UsymlqSolver for more details about the solver.