# LU factorization problem - Writing a code, don't understand partial pivoting

I'm trying to write a matlab code for the following question:

The program gets a matrix $A$ (lets say square matrix) and it returns $P,L,U$ such that $PA=PLU$ and $P$ is the permutation matrix, the partial pivoting matrix.

The obvious way to go about this is to: 1) find matrix $P$ via switching rows with the maximum value, and then doing gauss eliminiation, and continuing with the row switching.

2) doing regular LU decomposition of $A$

3) We now have $PA=PLU$

The problem here is that we essential do gauss elimination twice. once to find the $P$ matrix, and then to find $LU$ such that $A=LU$

The second problem, is that the entire point of partial pivoting is so we don't have to find such $LU$ matrices such that $A=LU$ since we don't know how numerically stable it is. So it defeats the point.

My problem is this: after switching rows, I don't know what value to put in matrix $L$.

For example:

$A=\begin{pmatrix} 3 & 1 & 4\\ 1 & 5 & 9 \\ 9 & 2 & 6\end{pmatrix}$

A row change is needed since $9>3$, so we now have:

$P=\begin{pmatrix} 0 & 0 & 1 \\ 0 & 1 & 0 \\ 1 & 0 & 0\end{pmatrix}$

and $PA = \begin{pmatrix} 9 & 2 & 6\\1 & 5 & 9\\ 3 & 1 & 4\end{pmatrix}$

Now I somehow need to find the original $L$ that I would get if I hadn't done the row change, by doing gauss elimination on $PA$ and I don't understand how that's possible.

Take heart, we need to find such matrices $L,U$ such that $A=LU$ without actually doing an lu decomposition of $A$, but of the more numerically stable $PA$ rather. It's kind of hard to explain

• The problem is that you're computing the wrong decomposition. Partial pivoting (at least in any implementation I've seen) produces a factorization where $PA=LU$, not $PA=PLU$. – Alex Becker Mar 24 '14 at 23:34
• I was specifically asked for $PA=PLU$...but even if I wanted to do $PA=LU$ how would I go about it without doing gauss elimination twice? (once to find matrix P, and then once to do LU decomposition of PA) – Oria Gruber Mar 24 '14 at 23:36

If it was indeed $PA = PLU$, then $P$ would be useless. Just premultiply with $P^{-1} = P^T$ and you get $A = LU$. No, it has to be $PA = LU$.

You do it in steps:

1. Find the pivot row and get it to the first position, say by swapping it with the first row. This will give you a permutation $P_1$.

2. Do the first column elimination of $P_1 A$ to obtain $$P_1 A = L_1 A_1,$$ where $L_1$ is lower triangular with a possibly non-zero first column and $1$ on its diagonal, while $A_1$ has zeros below the diagonal in the first column.

3. Repeat on the $A_1[2:n,2:n]$ (i.e., $A_1$ minus its first row and column) to obtain $P_2$, $L_2$, $A_2$, and so on.

4. When done, multiply (or smartly combine) all $P_k$ and $L_k$ (extended to dimension $n$ by identity) to get $P$ and $L$, respectively. You get $U$ by combining the first rows from all $A_k$.

You should probably check out Crout and LUP algorithms, where a similar - but I think a somewhat better - algorithm is described.

A bit on the stability of LU

But you see, $PA = LU$ is the $LU$ decomposition of $PA$. I want the $LU$ decomposition of $A$. But since $A$ is not stable, from what I understood, I need to find $LU$ such that $A = LU$ by doing an LU factorization on $PA$. As in - find the $LU$ factorization of $A$, by running the $LU$ factorization of the stable $PA$. does that not make sense?

Sensible as it may seem, it is not.

Take any $A$ such that $A_{1,1} = 0$ and you will not be able to do LU at all.

Now, take $A_{1,1}>0$, but very small, say $A_{1,1}=10^{-15}$, while all other elements of $A$ are non-zero integers. See what happens.

This is where the instability comes from and the result itself is unstable, so there is no algorithm that can fix it because any algorithm would have to get to that same final solution (LU is unique if $A$ is nonsingular).

By introducing the partial pivoting on the left side of $A$, you are avoiding this effect.

Effectively, these permutations are only permuting equations of your system, nothing more. So your LU of $PA$ is really the Gaussian eliminations algorithm performed on the same equations system, only these equations are in a different order.

Similarly, partial pivoting of columns would be the same as permuting the variables (thus introducing the need to be careful which variable gets which value in the end). The full pivoting would permute both the equations and the variables, for the most stable way to solve the system.

Variable-wise, all of these (when they can be performed) give you the same solution to a system of linear equations. But, in finite arithmetic, some solutions produce more errors than others. So, you'd rather get $x_1 = x_2$ than $x_1 = 10^{17}(1 + 10^{-17}x_2 - 1)$. Notice that the latter would actualy give you $x_1 = 0$ if $x_2 \approx 1$, even though it is (in exact arithmetic) the same solution as $x_1 = x_2$.

That's why $A = LU$ can be a problem, not just they way how you computed it.

• But you see, $PA=LU$ is the $LU$ decomposition of $PA$. I want the $LU$ decomposition of $A$. But since $A$ is not stable, from what I understood, I need to find $LU$ such that $A=LU$ by doing an LU factorization on $PA$. As in - find the $LU$ factorization of $A$, by running the $LU$ factorization of the stable $PA$. does that not make sense? – Oria Gruber Mar 25 '14 at 1:51
• @OriaGruber I've edited to explain the problem of numerical stability here. I hope it is more clear now. – Vedran Šego Mar 25 '14 at 2:13