summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs73
-rw-r--r--lib/Numeric/LinearAlgebra/Interface.hs55
-rw-r--r--lib/Numeric/LinearAlgebra/Linear.hs61
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs4
4 files changed, 67 insertions, 126 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 0f3ff0e..28063a9 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -61,9 +61,10 @@ module Numeric.LinearAlgebra.Algorithms (
61 nullspacePrec, 61 nullspacePrec,
62 nullVector, 62 nullVector,
63 nullspaceSVD, 63 nullspaceSVD,
64-- * Norms
65 Normed(..), NormType(..),
64-- * Misc 66-- * Misc
65 eps, peps, i, 67 eps, peps, i,
66 Normed(..), NormType(..),
67-- * Util 68-- * Util
68 haussholder, 69 haussholder,
69 unpackQR, unpackHess, 70 unpackQR, unpackHess,
@@ -172,9 +173,6 @@ thinSVD = {-# SCC "thinSVD" #-} thinSVD'
172singularValues :: Field t => Matrix t -> Vector Double 173singularValues :: Field t => Matrix t -> Vector Double
173singularValues = {-# SCC "singularValues" #-} sv' 174singularValues = {-# SCC "singularValues" #-} sv'
174 175
175instance (Field t, RealOf t ~ Double) => Norm2 Matrix t where
176 norm2 m = singularValues m @> 0
177
178-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. 176-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
179-- 177--
180-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. 178-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
@@ -546,7 +544,7 @@ epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]]
546geps delta = head [ k | (k,g) <- epslist, g<delta] 544geps delta = head [ k | (k,g) <- epslist, g<delta]
547 545
548expGolub m = iterate msq f !! j 546expGolub m = iterate msq f !! j
549 where j = max 0 $ floor $ log2 $ normInf m 547 where j = max 0 $ floor $ log2 $ pnorm Infinity m
550 log2 x = log x / log 2 548 log2 x = log x / log 2
551 a = m */ fromIntegral ((2::Int)^j) 549 a = m */ fromIntegral ((2::Int)^j)
552 q = geps eps -- 7 steps 550 q = geps eps -- 7 steps
@@ -569,7 +567,7 @@ expGolub m = iterate msq f !! j
569{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, 567{- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan,
570 based on a scaled Pade approximation. 568 based on a scaled Pade approximation.
571-} 569-}
572expm :: (Norm Vector t, Field t) => Matrix t -> Matrix t 570expm :: (Normed (Matrix t), Field t) => Matrix t -> Matrix t
573expm = expGolub 571expm = expGolub
574 572
575-------------------------------------------------------------- 573--------------------------------------------------------------
@@ -585,11 +583,11 @@ It only works with invertible matrices that have a real solution. For diagonaliz
585 [ 2.0, 2.25 583 [ 2.0, 2.25
586 , 0.0, 2.0 ]@ 584 , 0.0, 2.0 ]@
587-} 585-}
588sqrtm :: (Norm Vector t, Field t) => Matrix t -> Matrix t 586sqrtm :: (Normed (Matrix t), Field t) => Matrix t -> Matrix t
589sqrtm = sqrtmInv 587sqrtm = sqrtmInv
590 588
591sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) 589sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x))
592 where fixedPoint (a:b:rest) | norm1 (fst a |-| fst b) < peps = a 590 where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a
593 | otherwise = fixedPoint (b:rest) 591 | otherwise = fixedPoint (b:rest)
594 fixedPoint _ = error "fixedpoint with impossible inputs" 592 fixedPoint _ = error "fixedpoint with impossible inputs"
595 f (y,z) = (0.5 .* (y |+| inv z), 593 f (y,z) = (0.5 .* (y |+| inv z),
@@ -632,33 +630,56 @@ luFact (l_u,perm) | r <= c = (l ,u ,p, s)
632 630
633--------------------------------------------------------------------------- 631---------------------------------------------------------------------------
634 632
635data NormType = Infinity | PNorm1 | PNorm2 633data NormType = Infinity | PNorm1 | PNorm2 | Frobenius
636 634
637-- | Old class
638class Normed t where 635class Normed t where
639 pnorm :: NormType -> t -> Double 636 pnorm :: NormType -> t -> Double
640 637
641instance Norm Vector t => Normed (Vector t) where 638instance Normed (Vector Double) where
642 pnorm PNorm1 = realToFrac . norm1 639 pnorm PNorm1 = norm1
643 pnorm PNorm2 = realToFrac . normFrob 640 pnorm PNorm2 = norm2
644 pnorm Infinity = realToFrac . normInf 641 pnorm Infinity = normInf
642 pnorm Frobenius = normInf
643
644instance Normed (Vector (Complex Double)) where
645 pnorm PNorm1 = realPart . norm1
646 pnorm PNorm2 = realPart . norm2
647 pnorm Infinity = realPart . normInf
648 pnorm Frobenius = realPart . normInf
649
650instance Normed (Vector Float) where
651 pnorm PNorm1 = realToFrac . norm1
652 pnorm PNorm2 = realToFrac . norm2
653 pnorm Infinity = realToFrac . normInf
654 pnorm Frobenius = realToFrac . normInf
655
656instance Normed (Vector (Complex Float)) where
657 pnorm PNorm1 = realToFrac . realPart . norm1
658 pnorm PNorm2 = realToFrac . realPart . norm2
659 pnorm Infinity = realToFrac . realPart . normInf
660 pnorm Frobenius = realToFrac . realPart . normInf
661
645 662
646instance Normed (Matrix Double) where 663instance Normed (Matrix Double) where
647 pnorm PNorm1 = norm1 664 pnorm PNorm1 = maximum . map norm1 . toColumns
648 pnorm PNorm2 = norm2 665 pnorm PNorm2 = (@>0) . singularValues
649 pnorm Infinity = normInf 666 pnorm Infinity = pnorm PNorm1 . trans
667 pnorm Frobenius = norm2 . flatten
650 668
651instance Normed (Matrix (Complex Double)) where 669instance Normed (Matrix (Complex Double)) where
652 pnorm PNorm1 = norm1 670 pnorm PNorm1 = maximum . map (realPart.norm1) . toColumns
653 pnorm PNorm2 = norm2 671 pnorm PNorm2 = (@>0) . singularValues
654 pnorm Infinity = normInf 672 pnorm Infinity = pnorm PNorm1 . trans
673 pnorm Frobenius = realPart . norm2 . flatten
655 674
656instance Normed (Matrix Float) where 675instance Normed (Matrix Float) where
657 pnorm PNorm1 = realToFrac . norm1 676 pnorm PNorm1 = realToFrac . maximum . map norm1 . toColumns
658 pnorm PNorm2 = norm2 . double -- not yet optimized for Float 677 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
659 pnorm Infinity = realToFrac . normInf 678 pnorm Infinity = realToFrac . pnorm PNorm1 . trans
679 pnorm Frobenius = realToFrac . norm2 . flatten
660 680
661instance Normed (Matrix (Complex Float)) where 681instance Normed (Matrix (Complex Float)) where
662 pnorm PNorm1 = realToFrac . norm1 682 pnorm PNorm1 = realToFrac . maximum . map (realPart.norm1) . toColumns
663 pnorm PNorm2 = norm2 . double -- not yet optimized for Float 683 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
664 pnorm Infinity = realToFrac . normInf 684 pnorm Infinity = realToFrac . pnorm PNorm1 . trans
685 pnorm Frobenius = realToFrac . realPart . norm2 . flatten
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
index ec08694..13d175b 100644
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -19,7 +19,7 @@ In the context of the standard numeric operators, one-component vectors and matr
19----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
20 20
21module Numeric.LinearAlgebra.Interface( 21module Numeric.LinearAlgebra.Interface(
22 (<>),(<.>),mulG, Adapt, adaptElements, 22 (<>),(<.>),
23 (<\>), 23 (<\>),
24 (.*),(*/), 24 (.*),(*/),
25 (<|>),(<->), 25 (<|>),(<->),
@@ -29,28 +29,20 @@ import Data.Packed.Vector
29import Data.Packed.Matrix 29import Data.Packed.Matrix
30import Numeric.LinearAlgebra.Algorithms 30import Numeric.LinearAlgebra.Algorithms
31import Numeric.LinearAlgebra.Linear 31import Numeric.LinearAlgebra.Linear
32import Data.Complex
33import Control.Arrow((***))
34
35--import Numeric.GSL.Vector
36 32
37class Mul a b c | a b -> c where 33class Mul a b c | a b -> c where
38 infixl 7 <> 34 infixl 7 <>
39 -- | Matrix-matrix, matrix-vector, and vector-matrix products. 35 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
40 (<>) :: Product t => a t -> b t -> c t 36 (<>) :: Product t => a t -> b t -> c t
41 mulG :: (Element r, Element s, Adapt r s t t, Product t) => a r -> b s -> c t
42 37
43instance Mul Matrix Matrix Matrix where 38instance Mul Matrix Matrix Matrix where
44 (<>) = mXm 39 (<>) = mXm
45 mulG a b = uncurry mXm (curry adapt a b)
46 40
47instance Mul Matrix Vector Vector where 41instance Mul Matrix Vector Vector where
48 (<>) m v = flatten $ m <> (asColumn v) 42 (<>) m v = flatten $ m <> (asColumn v)
49 mulG m v = flatten $ m `mulG` (asColumn v)
50 43
51instance Mul Vector Matrix Vector where 44instance Mul Vector Matrix Vector where
52 (<>) v m = flatten $ (asRow v) <> m 45 (<>) v m = flatten $ (asRow v) <> m
53 mulG v m = flatten $ (asRow v) `mulG` m
54 46
55--------------------------------------------------- 47---------------------------------------------------
56 48
@@ -124,48 +116,3 @@ a <|> b = joinH a b
124-- -- | Vertical concatenation of matrices and vectors. 116-- -- | Vertical concatenation of matrices and vectors.
125-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t 117-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
126a <-> b = joinV a b 118a <-> b = joinV a b
127
128----------------------------------------------------
129
130class Adapt a b c d | a b -> c, a b -> d where
131 adapt :: Container k => (k a, k b) -> (k c, k d)
132
133--instance Adapt a a a a where
134-- adapt = id *** id
135
136instance Adapt Float Float Float Float where
137 adapt = id *** id
138
139instance Adapt Double Double Double Double where
140 adapt = id *** id
141
142instance Adapt (Complex Float) (Complex Float) (Complex Float) (Complex Float) where
143 adapt = id *** id
144
145instance Adapt Float Double Double Double where
146 adapt = double *** id
147
148instance Adapt Double Float Double Double where
149 adapt = id *** double
150
151instance Adapt Float (Complex Float) (Complex Float) (Complex Float) where
152 adapt = complex *** id
153
154instance Adapt (Complex Float) Float (Complex Float) (Complex Float) where
155 adapt = id *** complex
156
157instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt a (Complex Double) (Complex Double) (Complex Double) where
158 adapt = complex.double *** id
159
160instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt (Complex Double) a (Complex Double) (Complex Double) where
161 adapt = id *** complex.double
162
163instance Adapt Double (Complex Float) (Complex Double) (Complex Double) where
164 adapt = complex *** double
165
166instance Adapt (Complex Float) Double (Complex Double) (Complex Double) where
167 adapt = double *** complex
168
169adaptElements:: (Adapt a b c d, Container k) => (k a, k b) -> (k c, k d)
170adaptElements p = adapt p
171
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 6c21a16..4d7f608 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -24,8 +24,6 @@ module Numeric.LinearAlgebra.Linear (
24 Product(..), 24 Product(..),
25 mXm,mXv,vXm, 25 mXm,mXv,vXm,
26 outer, kronecker, 26 outer, kronecker,
27 -- * Norms
28 Norm(..), Norm2(..),
29 -- * Creation of numeric vectors 27 -- * Creation of numeric vectors
30 constant, linspace 28 constant, linspace
31) where 29) where
@@ -42,53 +40,64 @@ import Control.Monad(ap)
42class Num e => Vectors a e where 40class Num e => Vectors a e where
43 -- the C functions sumX are twice as fast as using foldVector 41 -- the C functions sumX are twice as fast as using foldVector
44 vectorSum :: a e -> e 42 vectorSum :: a e -> e
45 euclidean :: a e -> e
46 absSum :: a e -> e 43 absSum :: a e -> e
47 vectorMin :: a e -> e 44 vectorMin :: a e -> e
48 vectorMax :: a e -> e 45 vectorMax :: a e -> e
49 minIdx :: a e -> Int 46 minIdx :: a e -> Int
50 maxIdx :: a e -> Int 47 maxIdx :: a e -> Int
51 dot :: a e -> a e -> e 48 dot :: a e -> a e -> e
49 norm1 :: a e -> e
50 norm2 :: a e -> e
51 normInf :: a e -> e
52
52 53
53instance Vectors Vector Float where 54instance Vectors Vector Float where
54 vectorSum = sumF 55 vectorSum = sumF
55 euclidean = toScalarF Norm2 56 norm2 = toScalarF Norm2
56 absSum = toScalarF AbsSum 57 absSum = toScalarF AbsSum
57 vectorMin = toScalarF Min 58 vectorMin = toScalarF Min
58 vectorMax = toScalarF Max 59 vectorMax = toScalarF Max
59 minIdx = round . toScalarF MinIdx 60 minIdx = round . toScalarF MinIdx
60 maxIdx = round . toScalarF MaxIdx 61 maxIdx = round . toScalarF MaxIdx
61 dot = dotF 62 dot = dotF
63 norm1 = toScalarF AbsSum
64 normInf = vectorMax . vectorMapF Abs
62 65
63instance Vectors Vector Double where 66instance Vectors Vector Double where
64 vectorSum = sumR 67 vectorSum = sumR
65 euclidean = toScalarR Norm2 68 norm2 = toScalarR Norm2
66 absSum = toScalarR AbsSum 69 absSum = toScalarR AbsSum
67 vectorMin = toScalarR Min 70 vectorMin = toScalarR Min
68 vectorMax = toScalarR Max 71 vectorMax = toScalarR Max
69 minIdx = round . toScalarR MinIdx 72 minIdx = round . toScalarR MinIdx
70 maxIdx = round . toScalarR MaxIdx 73 maxIdx = round . toScalarR MaxIdx
71 dot = dotR 74 dot = dotR
75 norm1 = toScalarR AbsSum
76 normInf = vectorMax . vectorMapR Abs
72 77
73instance Vectors Vector (Complex Float) where 78instance Vectors Vector (Complex Float) where
74 vectorSum = sumQ 79 vectorSum = sumQ
75 euclidean = (:+ 0) . toScalarQ Norm2 80 norm2 = (:+ 0) . toScalarQ Norm2
76 absSum = (:+ 0) . toScalarQ AbsSum 81 absSum = (:+ 0) . toScalarQ AbsSum
77 vectorMin = ap (@>) minIdx 82 vectorMin = ap (@>) minIdx
78 vectorMax = ap (@>) maxIdx 83 vectorMax = ap (@>) maxIdx
79 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 84 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
80 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 85 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
81 dot = dotQ 86 dot = dotQ
87 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs
88 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapQ Abs
82 89
83instance Vectors Vector (Complex Double) where 90instance Vectors Vector (Complex Double) where
84 vectorSum = sumC 91 vectorSum = sumC
85 euclidean = (:+ 0) . toScalarC Norm2 92 norm2 = (:+ 0) . toScalarC Norm2
86 absSum = (:+ 0) . toScalarC AbsSum 93 absSum = (:+ 0) . toScalarC AbsSum
87 vectorMin = ap (@>) minIdx 94 vectorMin = ap (@>) minIdx
88 vectorMax = ap (@>) maxIdx 95 vectorMax = ap (@>) maxIdx
89 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 96 minIdx = minIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
90 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 97 maxIdx = maxIdx . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
91 dot = dotC 98 dot = dotC
99 norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs
100 normInf = (:+ 0) . vectorMax . fst . fromComplex . vectorMapC Abs
92 101
93---------------------------------------------------- 102----------------------------------------------------
94 103
@@ -268,39 +277,3 @@ kronecker a b = fromBlocks
268 $ flatten a `outer` flatten b 277 $ flatten a `outer` flatten b
269 278
270-------------------------------------------------- 279--------------------------------------------------
271
272-- | simple norms
273class (Element t, RealFloat (RealOf t)) => Norm c t where
274 norm1 :: c t -> RealOf t
275 normInf :: c t -> RealOf t
276 normFrob :: c t -> RealOf t
277
278instance Norm Vector Double where
279 normFrob = toScalarR Norm2
280 norm1 = toScalarR AbsSum
281 normInf = vectorMax . vectorMapR Abs
282
283instance Norm Vector Float where
284 normFrob = toScalarF Norm2
285 norm1 = toScalarF AbsSum
286 normInf = vectorMax . vectorMapF Abs
287
288instance (Norm Vector t, Vectors Vector t, RealElement t
289 , RealOf t ~ t, RealOf (Complex t) ~ t
290 ) => Norm Vector (Complex t) where
291 normFrob = normFrob . asReal
292 norm1 = norm1 . mapVector magnitude
293 normInf = vectorMax . mapVector magnitude
294
295instance Norm Vector t => Norm Matrix t where
296 normFrob = normFrob . flatten
297 norm1 = maximum . map norm1 . toColumns
298 normInf = norm1 . trans
299
300class Norm2 c t where
301 norm2 :: c t -> RealOf t
302
303instance Norm Vector t => Norm2 Vector t where
304 norm2 = normFrob
305
306-- (the instance Norm2 Matrix t requires singular values and is defined later)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index d6bb338..d312e52 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -54,9 +54,9 @@ complex x = complex'' x
54debug x = trace (show x) x 54debug x = trace (show x) x
55 55
56-- relative error 56-- relative error
57--dist :: (Normed t, Num t) => t -> t -> Double 57dist :: (Normed t, Num t) => t -> t -> Double
58dist a b = r 58dist a b = r
59 where norm = normInf 59 where norm = pnorm Infinity
60 na = norm a 60 na = norm a
61 nb = norm b 61 nb = norm b
62 nab = norm (a-b) 62 nab = norm (a-b)