Re: LUDecomposition

*To*: mathgroup at smc.vnet.net*Subject*: [mg55134] Re: LUDecomposition*From*: Maxim <ab_def at prontomail.com>*Date*: Sun, 13 Mar 2005 04:58:02 -0500 (EST)*References*: <d0rou0$2fv$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

On Fri, 11 Mar 2005 09:38:08 +0000 (UTC), DrBob <drbob at bigfoot.com> wrote: > I don't know if this is new to 5.1.1, and maybe I'm doing something > wrong myself, but this looks like a "problem" to me. (I'll avoid the B > word, for now.) > > I'm attempting to retrieve the appropriate upper and lower triangle > matrices from LUDecomposition and check that the multiplication actually > works, so that lower.upper is a permutation of the original matrix. > > The "errorSQ" function calculates data[[perm]] and lower.upper, and the > difference between the two. A sum-square-difference (it should be zero) > is the output. When the second argument (printFlag) is true, "errorSQ" > also prints data[[perm]] and lower.upper, so that we can compare them > side by side. > > test[f,k,n] calls "errorSQ" for an array of n k-by-k arrays, each of > them determined by f. > > (There's also a check for precisely singular matrices.) > > Clear[errorSQ, upperLU, lowerLU, test] > upperLU[m_?MatrixQ] := Transpose[m]Table[Boole[i <= j], {i, Length@m}, > {j, Length@m}] > lowerLU[m_?MatrixQ] := Transpose[m] - upperLU[m] > + IdentityMatrix@Length@m > errorSQ[data_?MatrixQ, printFlag_:False] := Module[{m, perm, ignored, > lower, upper, error}, > If[Det@data === 0, Infinity, > {m, perm, ignored} = LUDecomposition@data; > {lower, upper} = Through[{lowerLU, upperLU}@m]; > error = #.# &@Flatten[data[[perm]] - lower.upper]; > printFlag && Print /@ MatrixForm /@ {data[[perm]], lower.upper}; > error > ] > ] > test[f_, k_Integer, n_Integer] := errorSQ /@ Array[Array[f, {k, k}] &, > {n}] > > For real-valued matrix input, I always get a difference near zero. For > instance, here's the maximum squared error for 10,000 trials of 4 by 4 > matrices with random real entries between -5 and 5: > > Max@test[Random[Real,{-5,5}]&,4,10000] > > 1.5801870007708393*^-28 > > But when entries are Integer (leaving out the singular cases), the error > is never zero: > > Min@test[Random[Integer, {-5, 5}] &, 4, 10000] > > 30 > > I'm generating matrices and checking results the same way in both cases > (except Integer versus Real), and here's an example with printFlag=True: > > k=4; > errorSQ[Array[Random[Integer, {-5, 5}] &, {k, k}], True] > > MatrixForm[{{-1, 5, 3, -5}, > {0, -5, -4, 2}, {1, -5, 0, > -2}, {-5, -2, -1, 0}}] > MatrixForm[{{-1, 0, -1, 5}, > {-5, -5, -5, 152/5}, > {-3, 20, 0, -(71/15)}, > {5, -10, -16, 0}}] > 452482/225 > > data[[perm]] and lower.upper have different entries, not just the same > entries permuted in some way. > > So... what am I missing? > > Bobby > If all the values in the matrix returned by LUDecomposition are machine-precision numbers, then you have to transpose the matrix, otherwise the transposition is not needed (so, for example, if the output matrix contains an element 1.*x, then Transpose shouldn't be used). If this observation is correct, then Upper and Lower should be defined like this: In[1]:= Upper = If[MatrixQ[#, MachineNumberQ], Transpose[#], #]* Array[Boole[# <= #2]&, Dimensions@ #]&; Lower = If[MatrixQ[#, MachineNumberQ], Transpose[#], #] - Upper[#] + IdentityMatrix[Length@ #]&; In[3]:= A = {{0, 2, 2}, {5, 2, 4}, {5, 4, 0}}; Module[{LU = LUDecomposition[#]}, Lower[LU[[1]]].Upper[LU[[1]]] == #[[LU[[2]]]] ]& /@ {A, N[A], N[A, 20]} Out[4]= {True, True, True} I can't see any rational explanation for this. Maxim Rytin m.r at inbox.ru