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