Re: LU Decomposition w/o Pivoting
- To: mathgroup at smc.vnet.net
- Subject: [mg91734] Re: LU Decomposition w/o Pivoting
- From: Benjamin Reedlunn <breedlun at umich.edu>
- Date: Sun, 7 Sep 2008 05:32:40 -0400 (EDT)
- References: <g9t6m5$j89$1@smc.vnet.net> <48C2AF95.5080201@gmail.com>
On Sep 6, 2008, at 12:28 PM, Jean-Marc Gulliet wrote:
> Benjamin Reedlunn wrote:
>
>> I am trying to calculate the signs of the eigenvalues of a bunch
>> of moderately large matrices (thousand by thousand), which takes a
>> long time. It has come to my attention that the diagonal terms in
>> a LU decomposition have the same signs as the eigenvalues. So I'm
>> hoping that the LU decomposition is faster than calculating the
>> actual eigenvalues. I tried testing this out using mathematica's
>> built in function LUDecomposition[] and after a large amount of
>> monkeying around, I realized the signs of the diagonal terms only
>> match the eigenvalues if there is no pivoting employed.
>> Unfortunately, I cannot figure out how to turn off pivoting in
>> LUDecomposition[]. While I am sure I could write my own LU
>> decomposition function, I'm also sure it would not be nearly as
>> elegant and computationally efficient as the built in function.
>> So does anyone know how to turn off LUDecomposition[] or have
>> optimized code to perform a LU decomposition?
>
> FWIW,
>
> You might be able to use the second element returned by
> LUDecomposition for it is a vector specifying the rows used for
> pivoting.
>
> In[1]:= m = Table[N[1/(i + j + 1)], {i, 5}, {j, 5}];
> {lu, p, c} = LUDecomposition[m];
> lu // Diagonal
> p
> Eigenvalues[m]
>
> Out[3]= {0.333333, 0.0178571, -0.00111111, -0.0000340136, \
> -8.58929*10^-7}
>
> Out[4]= {1, 5, 2, 3, 4}
>
> Out[5]= {0.83379, 0.0430979, 0.00129982, 0.0000229962, 1.79889*10^-7}
>
> HTH,
> -- Jean-Marc
>
>
>
Thanks Jean-Marc for the suggestion, but I don't know if the pivoting
list returned by LUDecomposition[] is useful to me. The list is
generating during the construction of the LU decomposition and the
rows chosen to pivot totally change the composition of the L and U
matrices. The only way I can see using the pivoting list is by going
backwards through the algorithm for LUDecomposition[] and undo each
row swap.
I modified your code a bit below. I showed how I can reconstruct the
original matrix by multiplying Inverse[P].L.U, but I can't think of a
simple way to use the pivot list (or pivot matrix) to get the LU
decomposition that would have been calculated without pivoting.
In[315]:= n = 3;
m = Table[RandomInteger[10], {i, n}, {j, n}]
{lu, p, c} = LUDecomposition[m];
lu // Diagonal
Eigenvalues[m] // N
Out[316]= ( {
{2, 10, 9},
{2, 3, 10},
{2, 5, 3}
} )
Out[318]= {2, -5, 47/5}
Out[319]= {13.4853, -3.48528, -2.}
In[320]:= P = Table[0, {i, n}, {j, n}];
Do[P[[i, p[[i]]]] = 1, {i, n}];
L = (lu*SparseArray[{i_, j_} /; j < i -> 1, {n, n}] +
IdentityMatrix[n]);
U = lu SparseArray[{i_, j_} /; j >= i -> 1, {n, n}];
Inverse[P].L.U === m
Out[324]= True
In[356]:= Inverse[P].lu // Diagonal
lu.Inverse[P] // Diagonal
U.Inverse[P] // Diagonal
P.U // Diagonal
U.P // Diagonal
Out[356]= {2, 7/5, -6}
Out[357]= {2, -6, 7/5}
Out[358]= {2, -6, 0}
Out[359]= {2, 0, -6}
Out[360]= {2, -6, 0}
-Ben Reedlunn