Mathematica 9 is now available
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

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


  • Prev by Date: Re: Nested iterators in Compile
  • Next by Date: Re: Re: Determinant problem
  • Previous by thread: LUDecomposition
  • Next by thread: Simplfying inside Sqrt