problem with LUDecomposition and SparseArray

Questions
,
1

Hi everybody, I have the following problem, given a matrix lk

kl = {{1, 0, 0, 0, 0, 0},  {-1, -6, -2, 0, 0, 0}, {0, -2, -8, -2, 0, 0},   {0, 0, -2, 8, 6, 0},  {0, 0, 0, 6, -4, -8},  {0, 0, 0, 0, 0, 1}}

and do this

{lu, p, c} = LUDecomposition[kl]

L is obtained

l = lu SparseArray[{i_, j_} /; j < i -> 1, {6, 6}] + IdentityMatrix[6]

U is obtained

u = lu SparseArray[{i_, j_} /; j >= i -> 1, {6, 6}]

The problem arises when multiplying l.u

l.u

the result is not the matrix lk of above,Any idea why this happens?, thanks in advance

2

The rows are out of order

In[1]:= kl = {
   {1, 0, 0, 0, 0, 0},
   {-1, -6, -2, 0, 0,  0},
   {0, -2, -8, -2, 0, 0},
   {0, 0, -2, 8, 6, 0},
   {0, 0, 0, 6, -4, -8},
   {0, 0, 0, 0, 0, 1}};
{lu, p, c} = LUDecomposition[kl];
l = lu SparseArray[{i_, j_} /; j < i -> 1, {6, 6}] + IdentityMatrix[6];
u = lu SparseArray[{i_, j_} /; j >= i -> 1, {6, 6}];
l.u

Out[5]= {
  {1, 0, 0, 0, 0, 0},
  {0, -2, -8, -2, 0, 0},
  {0, 0, -2, 8, 6, 0},
  {0, 0, 0, 6, -4, -8},
  {-1, -6, -2, 0, 0, 0},
  {0, 0, 0, 0, 0, 1}}
3

LU decomp does not guarantee a unique solution. Click link here. The permutation matrix is possible: PA = LU if P is not identity matrix. The “A” is your “kl” in the problem.

Check:

(l.u).Inverse[kl] // MatrixForm

ScreenShot2014-06-12at2.23.56AM.png
ScreenShot2014-06-12at2.23.56AM.png324×160 4.71 KB

4

Bill Simpson thanks for the help, but I do not understand what you mean with The rows are out of order, is bad
built my matrix?,anxiously await your response in order to continue using this command

Shenghui Yang thanks for the link, is very explicit in the decomposition, and to say of the programs, just what I need, but that can also be done in Mathematica using the commands described above? I hope you can guide me on this topic

5

Shenghui Yang thanks for the link, is very explicit in the decomposition, and to say of the programs, just what I need, but that can also be done in Mathematica using the commands described above? I hope you can guide me on this topic

6

The pivoting matrix is treated at early stage of LU decomp/LUD, which in general is about gaussian elimination/GE. GE is the very fundamental O(n^2) scheme behind the LUD. This actually require some readings and exercise to understand. Here are two links can help you to some extend:

  1. partial pivoting
  2. full and partial pivoting

In Mathematica, the pivoting matrix is given by a permutation of the list [ 1,2,...,n ] , with n = dimension of the square matrix. In your case, kl is 6 by 6 so the n is 6. In your demonstation, you can see that p = [1,3,4,5,2,6] , which is a permutation of [1,2,3,4,5,6]. Because we are using the row pivoting, each number in p indicates the where the identity, aka. 1, is in each row. For instance, 3 is the second item in the list p, so there is a 1 at row 2 colume 3. Therefore you can use the following code to generate the pivoting matrix by yourself

pvtMat = SparseArray[
   Thread[
     MapThread[{#1, #2} &, {Range[6], p}] -> ConstantArray[1, 6]
]];

ScreenShot2014-06-13at1.23.49AM.png
ScreenShot2014-06-13at1.23.49AM.png524×252 9.99 KB

The matrix above should be identical the one generate by the following command in my last reply

(l.u).Inverse[kl] // MatrixForm