[Math] LU Factorization algorithm always fails.

matrices

I'm trying to implement LU factorization in openCL but I'm struggling to get my sequential algorithm working properly.

I implemented a sequential algorithm that works perfectly. Next I wanted to implement pivoting to enhance my precision. This is where I'm stuck. I can't figure out how to do it exactly. As a check for my calculations I use this website. From what I understand about pivoting is that each time you will swap rows such that the element with the biggest value in the pivot column is in the pivot position. Great.

So what I do in pseudocode for an MxM matrix, A is the following:

for each pivot index: 0 -> M:
    * Swap the current row with the row that has the biggest value in column M.
    * Swap the same row indices in P.
    * Calculate the new value for each element.
    * Recur

So as an example I took this input matrix:

239.000000      522.000000      287.000000      615.000000
373.000000       70.000000      816.000000      646.000000
766.000000      940.000000      482.000000      397.000000
347.000000       45.000000      909.000000      594.000000

We can clearly see that we have to swap row 0 with row 2, which results in the following matrix:

766.000000      940.000000      482.000000      397.000000
373.000000       70.000000      816.000000      646.000000
239.000000      522.000000      287.000000      615.000000
347.000000       45.000000      909.000000      594.000000

Now, when we do 1 iteration of decomposition we get the following result:

766.000000      940.000000      482.000000      397.000000
  0.000000     -387.728485      581.292419      452.682770
  0.000000      228.710175      136.610962      491.131866
  0.000000     -380.822449      690.652710      414.157959

Up to this point it matches with what the website gives me. So next iteration, we can see that row 2 now has the maximum value in the pivot column (228.710175) so we swap row 1 and 2, resulting in this matrix:

766.000000      940.000000      482.000000      397.000000
  0.000000      228.710175      136.610962      491.131866
  0.000000     -387.728485      581.292419      452.682770
  0.000000     -380.822449      690.652710      414.157959

And this is where it goes wrong. When I look at the result on the website I see this matrix:

 766.000  940.000  482.000  397.000
   0.000 -387.728  581.292  452.683
   0.000    0.000  479.499  758.157
   0.000    0.000    0.000 -219.747

It appears to be that they did not swap any rows in the second iteration. I have tried to test if PA = LU with my solution but it's wrong. Could anyone tell me what I am doing wrong?

Best Answer

You have to switch rows so the biggest absolute value is in the pivot position.

This is done because of two reasons:

  1. You avoid dividing by zero (if all pivot candidates are zero, the matrix is singular and no LU factorization exists).
  2. You enhance the numerical stability of your algorithm. Imagine you have a value 100.1, but your (imaginary) computer can only store integers so it stores 100 (you have a fault of 0.1). If you divide 100/2, you get 50 (it should be 50.05, you have a fault of 0.05). If you divide 100/-20, you get -5 (it should be -5.005, you have a fault of 0.005). As you can see you have smaller faults if you divide by bigger absolute numbers.

In the unlikely case you were unsure where the division takes place, it takes place in your lower triangular matrix. You divide the elements under the pivot by the pivot.

Related Question