Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: SQLSelect
  • Next by Date: Re: Canonical order...
  • Previous by thread: id back, please
  • Next by thread: Re: LUDecomposition