5
$\begingroup$

Ok, I've been working on a matrix solver for a while (no, its not going well, thanks for asking), and just discovered something that worrys me alot;

The problem at hand is LU Decomposition of a square matrix, for example

$\begin{matrix} 1&2&3\\ 2&-1&1\\ 3&4&-1\\ \end{matrix}$

WolframAlpha's answer is completely different from the example source

Wolfram

$\begin{matrix} 1&2&3\\ 3&-2&-10\\ 2&2.5&20\\ \end{matrix}$

Source

$\begin{matrix} 1&2&3\\ 2&-5&-5\\ 3&0.4&-8\\ \end{matrix}$

What's going on?

  • 0
    I think your example source has some problem. Compare equation (3) with equation (6). $A$ is supposed to equal $LU$...2011-04-17
  • 0
    The Wolfram one also ends up being trying to take the LU decomposition of a row permutation of $A$. (It is the one for with the bottom two rows switched.)2011-04-17
  • 1
    @Willie: The example source is correct. Equation (6) you're referring to is showing how to store the $LU$ decomposition efficiently in one matrix.2011-04-17
  • 0
    @Theo: ah, thanks, I didn't pay enough attention.2011-04-17
  • 2
    Also, I think the question statement has a typo. Entry $(2,2)$ of the Wolfram compact representation should be $-2$ instead of $-1$.2011-04-17
  • 0
    FWIW, you might profit from having a look at [chapter 3 of Golub/Van Loan](http://books.google.com/books?id=mlOa7wPX6OYC&pg=PA87) while you're attempting to implement Gaussian elimination.2011-04-18
  • 0
    @JM already have it open on the desk! but thanks.2011-04-18

1 Answers 1

9

It looks like Wolfram Alpha is doing an LU decomposition with partial pivoting. So, it's actually solving a system like $$\newcommand{\bm}[1]{\mathbf{#1}} \bm{P} \bm{A} = \bm{L} \bm{U} $$ where $\bm{P}$ is a permutation matrix.

Curiously, it appears that Alpha simply isn't listing the permutation used in solving the system. In this case it is $$ \bm{P} = \left( \begin{matrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{matrix} \right) . $$

This is a common technique to enhance numerical stability.

You can recover $\bm{A}$ via $$ \bm{A} = \bm{P}^{-1} \bm{L} \bm{U} = \bm{P L U} $$ since in this case $\bm{P} = \bm{P}^{-1}$.

  • 1
    So both are correct but WA is swapping around A before decomposition?2011-04-17
  • 0
    I didn't check closely the other solution, but the Wolfram one is correct when the permutation is included.2011-04-18