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