[Haskell-cafe] Re: recursive matrix computations

[Wilhelm B. Kloke]

Andrew U. Frank <frank at geoinfo.tuwien.ac.at> schrieb:
>  
> it appears possible to define matrix operations like LU-decomposition, SVD,
> inverse etc. in a recursive fashion (see transpose in the prelude).
>
> is this possible? has anybody code in haskell which does this (not asking
> for high performance here, more instructive value).

I would like to see this, too. Myself I did some experiments in trying
to implement the transformation to upper diagonal form by rational Givens
transforms. The following code fragment does this recursively:

-- rationalGivens :: Fractional a => a -> [a] -> [[a]] -> [[a]]
-- marginal cases first
rationalGivens qq [x] ((pp:[]):[]) = ( pp + qq * x * x : [] ) : []
rationalGivens _ [] [[]] = [[]]
rationalGivens qq xy pmat | qq == 0 = pmat
rationalGivens qq (x:y) (ppr:pmat) | x == 0 = ppr:rationalGivens qq y pmat
rationalGivens qq (x:y) [[]] = ( qq * x * x : (1/x).*y ) : []
-- main case
rationalGivens qq (x:y) ((pp:r):pmat) = let
        pp' = pp + qq * x * x
        qq' = pp * qq / pp'
        s = x * qq / pp'
        y' = y - x .* r
        r' = r + s .* y'
        in ( pp' : r' ) : rationalGivens qq' y' pmat

-- rationalGivens 1 [2,1] [[1,0],[1]] == [[5.0,0.4],[1.2]]

Arrays are double lists in this code,
q a scale factor (starting with 1)
xy a row vector to be added to the u.t. matrix pmat.

The diagonal elements of pmat contain the square of the scale factor of
the row, i.E. the matrix [[5.0,0.4],[1.2]] is to be interpreted as the product

        (sqrt(5)      0   ) ( 1  0.4 )
        (   0    sqrt(1.2)) ( 0   1  )

-- 
Dipl.-Math. Wilhelm Bernhard Kloke
Institut fuer Arbeitsphysiologie an der Universitaet Dortmund
Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257