# Solving $Ax=b$: Projection onto subspace with a canonical basis of largest error

The goal is to solve the linear system $$Ax = b$$, where $$A$$ is symmetric and positive definite (SPD). Consider the one-dimensional projection method given by equation (1):

$$x_{k+1} = \operatorname{argmin}_{x \in x_k + \text{span }{e_{i_k}}} || x - x^* ||_A, \quad (1)$$

where $$e_{i_k}$$ is the $$i_k$$-th canonical basis and $$x^*$$ is the actual solution. I need to use the following form of $$i_k$$: $$i_k = \operatorname{argmax}_{1 \leq i \leq n} | \langle x_k - x^*, e_i \rangle_A |, \quad (2)$$

Based on my understanding, this should mean update the equation corresponding to the largest error.

I've only dealt with simple 1-D projection of steepest descent so far so I'm not quite sure if I can simply replicate the steps. In steepest descent, the subspace would be span{$$b-Ax_k$$} instead of span{$$e_{i_k}$$} so the steps were a lot easier to understand. However, here the subspace and the choice of $$i_k$$ is confusing me a bit. I tried looking up stuff online but couldn't find anything relevant.

I would really appreciate if someone could help me understand what such a projection would mean and how I could develop an algorithm out of this idea. I would also appreciate any relevant reference materials which include algorithms/theory that deal with this specific type of projection.

Edit: After further digging, I've noticed that this is somewhat related to P5.1 from Yousef Saad's book 'Iterative Methods for Sparse Linear System' 2nd ed. Does this mean this is some additive/multiplicative type projection?

Similar to the steepest descent method, this is yet another example of a 1D projection method. The only difference is that we change a single component of $$x_k$$ to update the solution and this component is determined as that which has the largest projection of the error.

Since the error is unknown, we need to manipulate the expressions a bit to get computable expressions for the index $$i_k$$ and $$x_{k+1}$$. Computing the index $$i_k$$ is fairly easy: $$i_k = \mathop{\arg\max}_{1\leq i\leq n}|\langle x_k-x^*,e_i\rangle_A| =\mathop{\arg\max}_{1\leq i\leq n}|\langle A(x_k-x^*),e_i\rangle| =\mathop{\arg\max}_{1\leq i\leq n}|(r_k)_i|,$$ where $$r_k:=b-Ax_i=-A(x_k-x^*)$$ is the residual vector and $$(r_k)_i$$ is its $$i$$th component. It turns out, that $$i_k$$ is simply the index of the largest component of the residual $$r_k$$ in absolute value.

If $$x_{k+1}\in x_k+\mathrm{span}\{e_{i_k}\}$$, it can be written as $$x_{k+1}=x_k+\alpha_ke_{i_k},$$ where $$\alpha_k$$ is a scalar to be determined from the minimal error property (1). We have $$\begin{split} \|x_{k+1}-x^*\|_A^2&=\langle x_{k+1}-x^*,x_{k+1}-x^*\rangle_A \\&=\langle x_k+\alpha_ke_{i_k}-x^*,x_k+\alpha_ke_{i_k}-x^*\rangle_A \\&=\langle x_k-x^*,x_k-x^*\rangle_A+2\alpha_k\langle x_k-x^*,e_{i_k}\rangle_A+\alpha_k^2\langle e_{i_k},e_{i_k}\rangle_A, \end{split}$$ which is minimized for $$\alpha_k=-\frac{\langle x_k-x^*,e_{i_k}\rangle_A}{\langle e_{i_k},e_{i_k}\rangle_A} =\frac{(r_k)_{i_k}}{(A)_{i_k,i_k}}.$$

We can use the update formula for $$x_{k+1}$$ to determine the update formula for the residual as well: $$r_{k+1}=b-Ax_{k+1}=r_k-\alpha_kAe_{i_k}=r_k-\alpha_k (A)_{:,i_k}.$$

To summarize, the algorithm would go as follows: starting from some $$x_0$$, for $$k=0,1,2,...$$ until some convergence criterions is satisfied, do:

• pick the component $$i_k$$ of the residual vector $$r_k$$ largest in the magnitude,
• add $$\alpha_k=(r_k)_{i_k}/(A)_{i_k,i_k}$$ to the $$i_k$$th component of $$x_k$$ to compute $$x_{k+1}$$,
• subtract $$\alpha_k$$ multiple of the $$i_k$$ column of $$A$$ from $$r_k$$ to obtain $$r_{k+1}$$.

Note that we have $$\begin{split} \langle x_{k+1}-x^*,e_{i_k}\rangle_A &=\langle x_k+\alpha_ke_{i_k}-x^*,e_{i_k}\rangle_A \\&=\langle x_k-x^*,e_{i_k}\rangle_A+\alpha_k\langle e_{i_k},e_{i_k}\rangle_A \\&=-(r_k)_{i_k}+\alpha_k(A)_{i_k,i_k}=0. \end{split}$$ It means that (1) and (2) lead to a method such that the projection of the error of the new iteration is zero in the direction of the vector $$e_{i_k}$$.

I suppose this method is an analog of the Gauss-Seidel method with a "smarter" selection of the direction vector. This classical iterative method can be interpreted (see also the mentioned book by Saad) as a 1D projection method, which, in contrast to this method, sequentially zeroes out the $$i_k$$th component of the residual vector $$r_k$$, where $$i_k=(k\mod n)+1$$.