Iterative Solvers

 cgls(A, b; shift, abstol, reltol, maxiter, x0) -> (x, info)

Solve Ax = b or minimize norm(Ax-b) using CGLS, the conjugate gradient method for unsymmetric linear equations and least squares problems. A can be specified either as a rectangular matrix or as a linear operator, as defined in the LinearMaps package. It is desirable that eltype(A) == eltype(b), otherwise errors may result or additional allocations may occur in operator-vector products.

The keyword argument shift specifies a regularization parameter as shift = s. If s = 0 (default), then CGLS is Hestenes and Stiefel's specialized form of the conjugate-gradient method for least-squares problems. If s ≠ 0, the system (A'*A + s*I)*b = A'*b is solved.

An absolute tolerance abstol and a relative tolerance reltol can be specified for stopping the iterative process (default: abstol = 0, reltol = 1.e-6).

The maximum number of iterations can be specified using maxiter (default: maxiter = max(size(A),20)).

An initial guess for the solution can be specified using the keyword argument vector x0 (default: x0 = missing).

The resulting named tuple info contains (flag, resNE, iter), with convergence related information, as follows:

 `info.flag`  - convergence flag with values:  
                1, if convergence occured; 
                2, if the maximum number of iterations has been reached without convergence;
                3, if the matrix `A'*A + s*I` seems to be singular or indefinite;
                4, if instability seems likely meaning `(A'*A + s*I)` indefinite and `norm(x)` decreased;  

 `info.resNE` - the relative residual for the normal equations `norm(A'*b - (A'*A + s*I)*x)/norm(A'*b)`;  

 `info.iter`  - the iteration number at which `x` was computed.

This function is a translation of the MATLAB implementation of CGLS, the conjugate gradient method for nonsymmetric linear equations and least squares problems https://web.stanford.edu/group/SOL/software/cgls/. The author of the code is Michael Saunders, with contributions from Per Christian Hansen, Folkert Bleichrodt and Christopher Fougner.

Note: Two alternative solvers lsqr and lsmr, available in the IterativeSolvers package, can also be employed. For example, the following call to lsqr can be alternatively used:

  using IterativeSolvers
  lsqr(A, b; kwargs...) -> x[, history]

where kwargs contains solver-specific keyword arguments. A similar call to lsmr can be used.

lyapci(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the continuous Lyapunov matrix equation

            A*X + X*A' + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

lyapci(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized continuous Lyapunov equation

AXE' + EXA' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

plyapci(A, E, B; cyclic, abstol, reltol, maxiter, nshifts, shifts, loginf) -> (U, info)
plyapci(A, B; cyclic, abstol, reltol, maxiter, nshifts, shifts, loginf) -> (U, info)

Compute a low-rank factor U of the solution X = UU' of the generalized continuous Lyapunov equation

  AXE' + EXA' + BB' = 0,

where A and E are square real matrices and B is a real matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with negative real parts. E = I is assumed in the second call.

The named tuple info contains information related to the execution of the LR-ADI algorithm as follows: info.niter contains the number of performed iterations; info.res_fact contains the norm of the residual factor; info.res contains, if loginf = true, the vector of normalized residual norms (normalized with respect to the norm of the initial approximation); info.rc contains, if loginf = true, the vector of norms of relative changes in building the solution; info.used_shift contains the vector of used shifts.

The keyword argument abstol (default: abstol = 1e-12) is the tolerance used for convergence test on the normalized residuals, while the keyword argument reltol (default: reltol = 0) is the tolerance for the relative changes of the solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 100). The keyword argument nshifts specifies the desired number of shifts to be used in an iteration cycle (default: nshifts = 6). The keyword argument shifts can be used to provide a pre-calculated set of complex conjugated shifts to be used to start the iterations (default: shifts = missing). With the keyword argument loginf = true, the normalized residual values and the norms of increments of the solution can be saved as outputs in the resulting info structure (default: loginf = false).

The low-rank ADI (LR-ADI) method with enhancements proposed in [1] is implemented, based on MATLAB codes of the free software described in [2]. If cyclic = true, the cyclic low-rank method of [3] is used, with the pre-calculated shifts provided in the keyword argument shifts.

References

[1] P. Kürschner. Efficient Low-Rank Solution of Large-Scale Matrix Equations. Dissertation, Otto-von-Guericke-Universität, Magdeburg, Germany, 2016. Shaker Verlag,

[2] P. Benner, M. Köhler, and J. Saak. “Matrix Equations, Sparse Solvers: M-M.E.S.S.-2.0.1— Philosophy, Features, and Application for (Parametric) Model Order Reduction.” In Model Reduction of Complex Dynamical Systems, Eds. P. Benner et.al., 171:369–92, Springer, 2021.

[3] T. Penzl, A cyclic low-rank Smith method for large sparse Lyapunov equations, SIAM J. Sci. Comput. 21 (4) (1999) 1401–1418.

lyapdi(A, C; abstol, reltol, maxiter) -> (X,info)

Compute for a square A and a hermitian/symmetric C a solution X of the discrete Lyapunov matrix equation

            AXA' - X + C = 0.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

lyapdi(A, E, C; abstol, reltol, maxiter) -> (X,info)

Compute X, the solution of the generalized discrete Lyapunov equation

AXA' - EXE' + C = 0,

where A and E are square real or complex matrices and C is a square matrix. A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

plyapdi(A, E, B; cyclic, abstol, reltol, maxiter, nshifts, shifts, loginf) -> (U, info)
plyapdi(A, B; cyclic, abstol, reltol, maxiter, nshifts, shifts, loginf) -> (U, info)

Compute a low-rank factor U of the solution X = UU' of the generalized discrete Lyapunov equation

  AXA' - EXE' + BB' = 0,

where A and E are square real matrices and B is a real matrix with the same number of rows as A. The pencil A - λE must have only eigenvalues with moduli less than one. E = I is assumed in the second call.

The named tuple info contains information related to the execution of the LR-ADI algorithm as follows: info.niter contains the number of performed iterations; info.res_fact contains the norm of the residual factor; info.res contains, if loginf = true, the vector of normalized residual norms (normalized with respect to the norm of the initial approximation); info.rc contains, if loginf = true, the vector of norms of relative changes in building the solution; info.used_shift contains the vector of used shifts.

The keyword argument abstol (default: abstol = 1e-12) is the tolerance used for convergence test on the normalized residuals, while the keyword argument reltol (default: reltol = 0) is the tolerance for the relative changes of the solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 100). The keyword argument nshifts specifies the desired number of shifts to be used in an iteration cycle (default: nshifts = 6). The keyword argument shifts can be used to provide a pre-calculated set of complex conjugated shifts to be used to start the iterations (default: shifts = missing). With the keyword argument loginf = true, the normalized residual values and the norms of increments of the solution can be saved as outputs in the resulting info structure (default: loginf = false).

The low-rank ADI (LR-ADI) method with enhancements proposed in [1] is adapted to the discrete case, inspired by the MATLAB codes of the free software described in [2]. If cyclic = true, the cyclic low-rank method of [3] is adapted, with the pre-calculated shifts provided in the keyword argument shifts.

References

[1] P. Kürschner. Efficient Low-Rank Solution of Large-Scale Matrix Equations. Dissertation, Otto-von-Guericke-Universität, Magdeburg, Germany, 2016. Shaker Verlag,

[2] P. Benner, M. Köhler, and J. Saak. “Matrix Equations, Sparse Solvers: M-M.E.S.S.-2.0.1— Philosophy, Features, and Application for (Parametric) Model Order Reduction.” In Model Reduction of Complex Dynamical Systems, Eds. P. Benner et.al., 171:369–92, Springer, 2021.

[3] T. Penzl, A cyclic low-rank Smith method for large sparse Lyapunov equations, SIAM J. Sci. Comput. 21 (4) (1999) 1401–1418.

tlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous T-Lyapunov matrix equation

            A*X +isig*transpose(X)*transpose(A) = C   if adj = false,

or

            A*transpose(X) + isig*X*transpose(A) = C   if adj = true,

where for isig = 1, C is a symmetric matrix and for isig = -1, C is a skew-symmetric matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

hlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the continuous H-Lyapunov matrix equation

            A*X + isig*X'*A' = C   if adj = false,

or

            A*X' + isig*X*A' = C   if adj = true,

where for isig = 1, C is a hermitian matrix and for isig = -1, C is a skew-hermitian matrix.

For a matrix A, a least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

tulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a symmetric Q an upper triangular solution X of the continuous T-Lyapunov matrix equation

  transpose(U)*X + transpose(X)*U = Q   if adj = false,

or

  U*transpose(X) + X*transpose(U) = Q   if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

hulyapci(U, Q; adj = false, abstol, reltol, maxiter) -> (X,info)

Compute for an upper triangular U and a hermitian Q an upper triangular solution X of the continuous H-Lyapunov matrix equation

            U'*X + X'*U = Q   if adj = false,

or

            U*X' + X*U' = Q    if adj = true.

For a n×n upper triangular matrix U, a least-squares upper-triangular solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator L:X -> Y, which maps upper triangular matrices X into upper triangular matrices Y, and the associated matrix M = Matrix(L) is $n(n+1)/2 \times n(n+1)/2$. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

X = sylvci(A,B,C)

Solve the continuous Sylvester matrix equation

            AX + XB = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

X = sylvdi(A,B,C)

Solve the discrete Sylvester matrix equation

            AXB + X = C ,

where A and B are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

X = gsylvi(A,B,C,D,E)

Solve the generalized Sylvester matrix equation

AXB + CXD = E ,

where A, B, C and D are square matrices.

A least-squares solution X is determined using a conjugate gradient based iterative method applied to a suitably defined Lyapunov linear operator L:X -> Y such that L(X) = C or norm(L(X) - C) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

 gtsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized T-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*transpose(X)*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8

 ghsylvi(A, B, C, D, E; mx, nx, abstol, reltol, maxiter) -> (X,info)

Compute a solution X of the generalized H-Sylvester matrix equation

  ∑ A_i*X*B_i + ∑ C_j*X'*D_j = E,

where A_i and C_j are matrices having the same row dimension equal to the row dimension of E and B_i and D_j are matrices having the same column dimension equal to the column dimension of E. A_i and B_i are contained in the k-vectors of matrices A and B, respectively, and C_j and D_j are contained in the l-vectors of matrices C and D, respectively. Any of the component matrices can be given as an UniformScaling. The keyword parameters mx and nx can be used to specify the row and column dimensions of X, if they cannot be inferred from the data.

A least-squares solution X is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Sylvester linear operator L:X -> Y such that L(X) = E or norm(L(X) - E) is minimized. The keyword arguments abstol (default: abstol = 0) and reltol (default: reltol = sqrt(eps())) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument maxiter can be used to set the maximum number of iterations (default: maxiter = 1000).

Note: For the derivation of the adjoint equation see reference [1], which also served as motivation to implement a general linear matrix equation solver in Julia.

[1] Uhlig, F., Xu, A.B. Iterative optimal solutions of linear matrix equations for hyperspectral and multispectral image fusing. Calcolo 60, 26 (2023). https://doi.org/10.1007/s10092-023-00514-8