MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Mapping A Condition Across Two Different Lists
  • Next by Date: Re: Thinking Mathematica: Any suggestions?
  • Previous by thread: Re: LU Decomposition w/o Pivoting
  • Next by thread: How can I do a "grep -v" equivalent in Import[]?