[Math] LU factorization problem – Writing a code, don’t understand partial pivoting

gaussian eliminationlinear algebramatricesmatrix decomposition

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

Best Answer

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

A question in the comments:

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.

Related Question