This can be done using a technique called the $\eta$ (eta) factorization. This method is commonly used in implementations of the simplex method for linear programming and can be found in many textbooks on linear programming.
The procedure is as follows:
- Find an LU factorization of $A_{1}$,
$PA_{1}=LU$
Using this factorization, you can easily solve $A_{1}x=b$.
- Suppose that $A_{2}$ is constructed from $A_{1}$ by replacing column
$p$ of $A_{1}$ by a vector $v$. Then we can write
$A_{2}=A_{1}E_{2}$
where $E_{2}$ is a so-called $\eta$ matrix that is an identity matrix with column $p$ replaced by the solution to $A_{1}u=v$. When we multiply $A_{1}$ times the columns of $E_{2}$ taken from $I$, we get the original $A_{1}$ columns copied into $A_{2}$. When we multiply $A_{1}$ times $u$ in the $p$th column of $E_{2}$, we get the desired vector $v$.
- To solve $A_{2}x=b$, we use the following procedure. First, solve $A_{1}w=b$ using the original LU factorization of $A_{1}$. Then, solve $E_{2}x=w$ by first solving for $x_{p}=w_{p}/u_{p}$ and then using back substitution to solve for the remaining entries in $x$. Thus
$A_{2}x=A_{1}E_{2}x=A_{1}w=b$.
This process is continued to get
$A_{k}=A_{1}E_{2}E_{3}E_{4}\cdots E_{k}$.
The worst-case computational complexity of the procedure is as follows. Computing the LU factorization of $A_{1}$ takes $O(n^{3})$ time. Computing an $\eta$ matrix takes $O(n^{2})$ time. The solution of the system $A_{1}w=b$ takes $O(n^{2}$ time. The solution of a system $E_{j}v_{j}=w_{j-1}$ takes $O(n)$ time.
In practice, the matrix is usually refactorized once the number of $\eta$ matrices hits a configurable limit. In linear programming, refactorizing the basis matrix after 30 $\eta$ updates is a common rule of thumb.
Although this method is commonly employed with a direct factorization method to solve $A_{1}x=b$, you could apply the same approach with an iterative method for solving $A_{1}x=b$. Finding each $\eta$ matrix $E_{k}$ would require a full iterative solution plus $k-1$ $\eta$ updates. The solution of $A_{k}x=b$ would require one full iterative solution, plus $k$ $\eta$ updates, each of which takes only $O(n)$ time. This really wouldn't save any computational effort compared to solving $A_{k}x=b$ using the iterative method.