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).

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).

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