The question is from section 2.5 Iterative Improvement of a Solution to Linear Equations in Numerical Recipes book.
When we solve $\mathbf{Ax = b}$ using LU decomposition numerically, the result is not exact due to the round off error accumulation. So we use a trick called interative improvement of the solution to restore the full machine precision. Suppose that a vector $\mathbf{x}$ is the exact solution of the linear set $$\mathbf{Ax = b}$$ You don't however know $\mathbf{x}$. You only know slightly wrong solution $\mathbf{x} + \delta\mathbf{x}$, where $\delta\mathbf{x}$ is the unknown error, when multiplying by $\mathbf{A}$, it becomes $$\mathbf{A\cdot(x + \delta x) = b + \delta b}\tag{1}$$ Then we have $$\mathbf{A\cdot\delta x = \delta b}\tag{2}$$ $\delta b$ is easy to get by equation (1), thus we can iterative improve the result using (1) and (2).
But we neglected the fact that the LU decomposition of $\mathbf{A}$ is itself exact. A different analytical approach starts with some matrix $\mathbf{B_0}$ that is assumed to be an approximate inverse of the matrix $\mathbf{A}$, so that $\mathbf{B_0\cdot A}$ is approximately the identity matrix $\mathbf{1}$. Define the residual matrix $\mathbf{R}$ of $\mathbf{B_0}$ as $$\mathbf{R\equiv 1 - B_0\cdot A}$$ which is supposed to be "small", therefore $$\mathbf{B_0\cdot A = 1 - R}$$ Next consider the following formal manipulation $$\mathbf{A^{-1} = A^{-1}\cdot (B_0^{-1}\cdot B_0) = (A^{-1}\cdot B_0^{-1})\cdot B_0 = {(B_0\cdot A)}^{-1}\cdot B_0 \\ = {(1 - R)}^{-1}\cdot B_0 = (1 + R + R^2 + R^3 + \cdots)\cdot B_0}$$ We can define the nth partial sum of the last expression by $$\mathbf{B_n \equiv (1 + R + \cdots + R^n)\cdot B_0}\tag{3}$$ so that $\mathbf{B_\infty \to A^{-1}}$, if the limit exits.
It now is straight to verify that equation (3) satisfies some interesting recurrence relations. As regards solving $\mathbf{Ax = b}$, where $\mathbf{x}$ and $\mathbf{b}$ are vectors, define $$\mathbf{x_n \equiv B_n\cdot b}\tag{4}$$ It can be show that $$\mathbf{x_{n+1} = x_n + B_0 \cdot (b - A\cdot x_n)}\tag{5}$$ Here is what I cannot understand, how to deduce equation (5) from above equations.
And the book implies that $$\mathbf{B_{2n+1} = 2B_n - B_n\cdot A\cdot B_n} \qquad n = 0, 1, 3, 7 ...$$ can be drived from equation (3), and I don't know how. Can some one make these two deduction more clear, Thank you!
For (5), we can proceed as follows. You can check using (3) that $$ (I-R^{n+1})B_0=(I-R)B_n=B_0AB_n. $$ Then we have $$ \begin{split} x_{n+1}-x_n&=(B_{n+1}-B_n)b=R^{n+1}B_0b \\ &=[I-(I-R^{n+1})]B_0b=B_0b-B_0AB_nb\\ &=B_0(b-Ax_n). \end{split} $$
For the last identity, we have $$ \begin{split} B_nAB_n&=(I+R+\cdots+R^n)B_0A(I+R+\cdots+R^n)B_0 \\&=(I+R+\cdots+R^n)(I-R)(I+R+\cdots+R^n)B_0 \\&=(I+R+\cdots+R^n)(I-R^{n+1})B_0 \\&=(I+R+\cdots+R^n-R^{n+1}-\cdots-R^{2n+1})B_0 \\&=[2(I+R+\cdots+R^n)-(I+R+\cdots+R^{2n+1})]B_0 \\&=2B_n-B_{2n+1}. \end{split} $$ Note that the identity holds for any non-negative integer.