       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)
• 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

```

• Prev by Date: Re: SQLSelect
• Next by Date: Re: Canonical order...