Solve $A^{-1} b$ when one column is replaced

Question

Given square matrix $A_0$, vector $b$, vector $A_0^{-1}b$ and matrices $A_1, A_2, \dots, A_k$, in which each $A_i$ is generated from $A_{i-1}$ by replacing one single column, I would like to find an efficient method to compute all the results of

$$A_1^{-1}b, A_2^{-1}b, \dots, A_k^{-1}b$$

Any suggestions are welcome.

Answer 1

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:

1. Find an LU factorization of $A_{1}$,

$PA_{1}=LU$

Using this factorization, you can easily solve $A_{1}x=b$.

2. 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$.

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