Conjugate gradient squared method

In numerical linear algebra, the conjugate gradient squared method (CGS) is an iterative algorithm for solving systems of linear equations of the form $$A{\bold x} = {\bold b}$$, particularly in cases where computing the transpose $$A^T$$ is impractical. The CGS method was developed as an improvement to the biconjugate gradient method.

Background
A system of linear equations $$A{\bold x} = {\bold b}$$ consists of a known matrix $$A$$ and a known vector $${\bold b}$$. To solve the system is to find the value of the unknown vector $${\bold x}$$. A direct method for solving a system of linear equations is to take the inverse of the matrix $$A$$, then calculate $$\bold x = A^{-1}\bold b$$. However, computing the inverse is computationally expensive. Hence, iterative methods are commonly used. Iterative methods begin with a guess $$\bold x^{(0)}$$, and on each iteration the guess is improved. Once the difference between successive guesses is sufficiently small, the method has converged to a solution.

As with the conjugate gradient method, biconjugate gradient method, and similar iterative methods for solving systems of linear equations, the CGS method can be used to find solutions to multi-variable optimisation problems, such as power-flow analysis, hyperparameter optimisation, and facial recognition.

Algorithm
The algorithm is as follows:


 * 1) Choose an initial guess $${\bold x}^{(0)}$$
 * 2) Compute the residual $${\bold r}^{(0)} = {\bold b} - A{\bold x}^{(0)}$$
 * 3) Choose $$\tilde {\bold r}^{(0)} = {\bold r}^{(0)}$$
 * 4) For $$i = 1, 2, 3, \dots$$ do:
 * 5) $$\rho^{(i-1)} = \tilde {\bold r}^{T(i-1)}{\bold r}^{(i-1)}$$
 * 6) If $$\rho^{(i-1)} = 0$$, the method fails.
 * 7) If $$i=1$$:
 * 8) $${\bold p}^{(1)} = {\bold u}^{(1)} = {\bold r}^{(0)}$$
 * 9) Else:
 * 10) $$\beta^{(i-1)} = \rho^{(i-1)}/\rho^{(i-2)}$$
 * 11) $${\bold u}^{(i)} = {\bold r}^{(i-1)} + \beta_{i-1}{\bold q}^{(i-1)}$$
 * 12) $${\bold p}^{(i)} = {\bold u}^{(i)} + \beta^{(i-1)}({\bold q}^{(i-1)} + \beta^{(i-1)}{\bold p}^{(i-1)})$$
 * 13) Solve $$M\hat {\bold p}={\bold p}^{(i)}$$, where $$M$$ is a pre-conditioner.
 * 14) $$\hat {\bold v} = A\hat {\bold p}$$
 * 15) $$\alpha^{(i)} = \rho^{(i-1)} / \tilde {\bold r}^T \hat {\bold v}$$
 * 16) $${\bold q}^{(i)} = {\bold u}^{(i)} - \alpha^{(i)}\hat {\bold v}$$
 * 17) Solve $$M\hat {\bold u} = {\bold u}^{(i)} + {\bold q}^{(i)}$$
 * 18) $${\bold x}^{(i)} = {\bold x}^{(i-1)} + \alpha^{(i)} \hat {\bold u}$$
 * 19) $$\hat {\bold q} = A\hat {\bold u}$$
 * 20) $${\bold r}^{(i)} = {\bold r}^{(i-1)} - \alpha^{(i)}\hat {\bold q}$$
 * 21) Check for convergence: if there is convergence, end the loop and return the result