Message Boards Message Boards

0
|
7651 Views
|
5 Replies
|
4 Total Likes
View groups...
Share
Share this post:

problem with LUDecomposition and SparseArray

Posted 10 years ago

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

POSTED BY: Luis Ledesma
5 Replies

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]
]];

piv mat

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

(l.u).Inverse[kl] // MatrixForm
POSTED BY: Shenghui Yang

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

permutation matrix

POSTED BY: Shenghui Yang
Posted 10 years ago

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

POSTED BY: Luis Ledesma
Posted 10 years ago

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}}
POSTED BY: Bill Simpson
Posted 10 years ago

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

POSTED BY: Luis Ledesma
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract