[Haskell-cafe] Re: recursive matrix computations
Wilhelm B. Kloke
wb at arb-phys.uni-dortmund.de
Wed Apr 19 07:31:25 EDT 2006
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],] == [[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
More information about the Haskell-Cafe