Sometimes it helps to inspect intermediate results. With recursion,
though, it can be a bit tricky for a casual observer to see the
intermediate results. With that in mind, here's what I am seeing for
your example:
a1=: (calcU calcL) saveAA
a2=: (calcU calcL) a1
a3=: (calcU calcL) a2
A4=: ((({.@[ ,: ]) ,&.:(|."1) a3"_) calcL) a2
A5=: ((({.@[ ,: ]) ,&.:(|."1) A4"_) calcL) a1
A6=: ((({.@[ ,: ]) ,&.:(|."1) A5"_) calcL) saveAA
a1, a2 and a3 are progressively smaller square matrices (2x2, 1x1, 0x0)
A4, A5 and A6 are progressively larger matrices which are twice as
tall as wide. If you could compute them in reverse order it might have
made sense to make it twice as wide as tall (with intermediate lu side
by side instead of interleave stacked)?
A6 is the same as lumain saveAA
I should go back and re-read km's implementation. But I will note that
you can cut code size slightly using some cross hooks:
lumain =: (((,:~ {.)~ ,&.:(|."1) $:@calcU) calcL)^:(*@#)
lu =: [: (,:~ |:)/ 1 0 2 |: _2 ]\ lumain
Anyways, I think your O(n^3) space is largely because all intermediate
values from what I have characterized as a (calcU calcL) hook are
"pre"-computed and placed on the stack before proceeding with further
computations.
Thanks,
--
Raul
On Sun, Jan 12, 2014 at 10:10 PM, Henry Rich <henryhr...@nc.rr.com> wrote:
> calcL =: (% {.)@:({."1)
> calcU =: (}.@[ - {.@[ *"1 0 ])&:(}."1)
> lumain =: ((({.@[ ,: ]) ,&.:(|."1) $:@calcU) calcL)^:(*@#)
> lu =: [: (|:@] ,: [)/ 1 0 2 |: _2 ]\ lumain
> NB. Half this code is handling joining ragged lists.
> NB. Is there a better way?
>
> saveAA =: 3 3 $ 2 1 4 _4 _1 _11 2 4 _2
> lu saveAA
>
> 1 0 0
> _2 1 0
> 1 3 1
>
> 2 1 4
> 0 1 _3
> 0 0 3
>
> I suspect that a vectorized explicit version is a better way to go. This
> version has memory requirements of O(n^3).
>
> Henry Rich
>
>
> On 1/12/2014 9:00 PM, km wrote:
>>
>> Verb LU below produces the matrices L and U of the LU decomposition of a
>> square matrix A. L is lower triangular, U is upper triangular, and A is L
>> +/ . * U .
>>
>> Should one attempt a tacit version?
>>
>> eye =: =@i.@] NB. eye 3 is a 3 by 3 identity matrix
>>
>> rop =: 3 : 0 NB. row op: subtract c times row i0 from row i1
>> :
>> 'i1 c i0' =. x
>> ( (i1 { y) - c * i0 { y ) i1 } y
>> )
>>
>> LU =: 3 : 0 NB. square matrices L and U for y -: L +/ . * U
>> m =. # y
>> L =. eye(m)
>> U =. y
>> for_j. i. <: m do.
>> p =. (< j , j) { U
>> for_i. j + >: i. <: m - j do.
>> c =. p %~ (< i , j) { U
>> L =. c (< i , j) } L
>> U =. (i, c, j) rop U
>> end.
>> end.
>> L ,: U
>> )
>>
>> saveAA
>> 2 1 4
>> _4 _1 _11
>> 2 4 _2
>>
>> LU saveAA
>> 1 0 0
>> _2 1 0
>> 1 3 1
>>
>> 2 1 4
>> 0 1 _3
>> 0 0 3
>>
>> saveAA -: +/ . */ LU saveAA
>> 1
>>
>> --Kip Murray
>>
>> Sent from my iPad
>> ----------------------------------------------------------------------
>> For information about J forums see http://www.jsoftware.com/forums.htm
>>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
