3
$\begingroup$

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.

Improve this question
$\endgroup$
  • $\begingroup$ How do you get a solution of $A_0^{-1}b$? Do you perform LU-decomposition on $A_0$? For example, see scicomp.stackexchange.com/questions/25063/… $\endgroup$ – Anton Menshov♦ Apr 12 '18 at 17:40
  • $\begingroup$ I use the iterative method: $x_{j+1} = b - (A_0-I)x_{j}$. $\endgroup$ – John Smith Apr 12 '18 at 17:48
  • $\begingroup$ @RodrigodeAzevedo, By viewing the replacement of one single column as a rank-one update of $A$, i.e., $A' = A+uv^T$. I notice that, after applying Sherman–Morrison method, I still need to compute $A^{-1}u$. It appears to me that the cost won't be reduced. $\endgroup$ – John Smith Apr 12 '18 at 17:52
  • $\begingroup$ @RodrigodeAzevedo, he is certainly not finding $A^{-1}$ explicitly. I'm not sure how can you use Sherman-Morrison for the iterative procedure... $\endgroup$ – Anton Menshov♦ Apr 12 '18 at 18:09
  • 1
    $\begingroup$ If you're using a Krylov subspace method, then in the process of computing a solution $A_0^{-1}b$ you are also building up an approximation of $A_0^{-1}$ in the Krylov subspace (e.g. see any reference on the Arnoldi process). This approximation is almost certainly going to be a great preconditioner for $A_1$. In fact I'm fairly sure that you can prove something about the quality of this preconditioner. $\endgroup$ – Richard Zhang Apr 12 '18 at 18:39
6
$\begingroup$

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

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

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

Improve this answer
$\endgroup$
  • $\begingroup$ Thanks, Brian. It seems that solving $A_1 u =v$ for $u$ has the same computational complexity as directly solving $A_1 x =b$ for $x$. What's the main advantage of your method? $\endgroup$ – John Smith Apr 14 '18 at 1:31
  • $\begingroup$ The important point is that with a direct factorization, it takes $O(n^{3})$ time to factor $A_{1}$, and then each additional solve takes only $O(n^{2})$ time. This advantage disappears if you use an iterative method to solve the equations. $\endgroup$ – Brian Borchers Apr 14 '18 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.