LUDecomposition

*To*: mathgroup at smc.vnet.net*Subject*: [mg55074] LUDecomposition*From*: DrBob <drbob at bigfoot.com>*Date*: Fri, 11 Mar 2005 04:20:59 -0500 (EST)*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

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 -- DrBob at bigfoot.com