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