# Solve rank one update to LU using plain vanilla LU routine

I have a large number of systems of the form: $$A_ix=b,$$ where each $A_i,i>0$ is a rank one update of $A_{i-1}$ and the $A_i$ are dense matrices.

I was wondering whether it is possible to use the classical LU algorithm to solve these (rank one updated) versions of the initial problem more efficiently?

P.S. Right now I'm using eigen c++'s PartialPivLU implementation.

• Does the rank one update have some special form (such as replacement of a column of A), or is it a completely general rank one update $A_{i}=A_{i-1}+uv^{T}$? – Brian Borchers Nov 14 '15 at 3:57
• How many matrices are there? – Kirill Nov 14 '15 at 7:11
• @BrianBorchers: each update replaces a row of $A_i$. Kirill: about 25. – user189035 Nov 14 '15 at 9:16
• A technique like the eta factorization used in the simplex method for column replacements could be devised for your row updates. – Brian Borchers Nov 14 '15 at 14:37

The eta factorization of the basis is a technique widely used in the simplex method for LP to handle rank one column updates. The same idea can be modified to handle rank one row updates,

I'll assume that you start with an LU factorization of $A_{0}$ and can solve systems of equations of the form $A_{0}x=b$. We can also solve systems of the form $yA_{0}=c$ where $c$ is a given row vector and $y$ is the row vector that you are solving for. If $PA_{0}=LU$, then we can solve $yA_{0}=c$ by

1. Solve $wU=c$.
2. Solve $vL=w$.
3. Let $y=vP$.

Suppose that in the first iteration, row $k$ of $A_{0}$ is to be replaced by a row vector $c_{1}$. Solve

$yA_{0}=c_{1}$

Let $E_{1}$ be the identity matrix with $y$ in place of row row $k$. Then

$E_{1}A_{0}=A_{1}$.

Now, we can solve a system of equations $A_{1}x=b$ by

1. Solve $E_{1}u=b$.
2. Solve $A_{0}x=u$.

Note that because $E_{1}$ is nearly the identity matrix except for row $k$, it's trivial to solve $E_{1}u=b$.

You can repeat this process to get

$E_{2}E_{1}A_{0}=A_{2}$

$\ldots$

$E_{k} \cdots E_{1}A_{0}=A_{k}$.

In the simplex method it is typical to do this for something like 30 iterations before refactoring $A_{30}$ from scratch.

• A related technique computes $A_{k}^{-1}=A_{0}^{-1}E_{1}^{-1} \cdots E_{k}^{-1}$. This product form of the inverse can be convenient if you have lots of right hand sides $b$, but $A_{0}^{-1}$ is likely to be dense. – Brian Borchers Nov 14 '15 at 18:49