-norm, +rcond, +kronecker

2 files changed, 39 insertions, 4 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 0683956..52f9b6f 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -22,7 +22,7 @@ module Numeric.LinearAlgebra.Algorithms (
 -- * Linear Systems
     linearSolve,
     inv, pinv,
-    pinvTol, det, rank,
+    pinvTol, det, rank, rcond,
 -- * Matrix factorizations
 -- ** Singular value decomposition
     svd,
@@ -244,10 +244,10 @@ pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (liftVector
 pnormCM Infinity m = vectorMax $ liftMatrix (liftVector magnitude) m `mXv` constant 1 (cols m)
 --pnormCM _ _ = error "p norm not yet defined"
 
+-- | Objects which have a p-norm.
+-- Using it you can define convenient shortcuts: @norm2 = pnorm PNorm2@, @frobenius = norm2 . flatten@, etc.
 class Normed t where
     pnorm :: NormType -> t -> Double
-    norm :: t -> Double
-    norm = pnorm PNorm2
 
 instance Normed (Vector Double) where
     pnorm = pnormRV
@@ -356,6 +356,12 @@ uH (pq, tau) = (p,h)
 
 --------------------------------------------------------------------------
 
+-- | Reciprocal of the 2-norm condition number of a matrix, computed from the SVD.
+rcond :: GenMat t => Matrix t -> Double
+rcond m = last s / head s
+    where (_,s',_) = svd m
+          s = toList s'
+
 -- | Number of linearly independent rows or columns.
 rank :: GenMat t => Matrix t -> Int
 rank m | pnorm PNorm1 m < eps = 0
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
index 3017936..94f6958 100644
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ b/lib/Numeric/LinearAlgebra/Linear.hs
@@ -16,7 +16,7 @@ Basic optimized operations on vectors and matrices.
 
 module Numeric.LinearAlgebra.Linear (
     Linear(..),
-    multiply, dot, outer
+    multiply, dot, outer, kronecker
 ) where
 
 
@@ -100,3 +100,32 @@ dot u v = dat (multiply r c) `at` 0
 -}
 outer :: (Field t) => Vector t -> Vector t -> Matrix t
 outer u v = asColumn u `multiply` asRow v
+
+{- | Kronecker product of two matrices.
+
+@m1=(2><3)
+ [ 1.0,  2.0, 0.0
+ , 0.0, -1.0, 3.0 ]
+m2=(4><3)
+ [  1.0,  2.0,  3.0
+ ,  4.0,  5.0,  6.0
+ ,  7.0,  8.0,  9.0
+ , 10.0, 11.0, 12.0 ]@
+
+@\> kronecker m1 m2
+(8><9)
+ [  1.0,  2.0,  3.0,   2.0,   4.0,   6.0,  0.0,  0.0,  0.0
+ ,  4.0,  5.0,  6.0,   8.0,  10.0,  12.0,  0.0,  0.0,  0.0
+ ,  7.0,  8.0,  9.0,  14.0,  16.0,  18.0,  0.0,  0.0,  0.0
+ , 10.0, 11.0, 12.0,  20.0,  22.0,  24.0,  0.0,  0.0,  0.0
+ ,  0.0,  0.0,  0.0,  -1.0,  -2.0,  -3.0,  3.0,  6.0,  9.0
+ ,  0.0,  0.0,  0.0,  -4.0,  -5.0,  -6.0, 12.0, 15.0, 18.0
+ ,  0.0,  0.0,  0.0,  -7.0,  -8.0,  -9.0, 21.0, 24.0, 27.0
+ ,  0.0,  0.0,  0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
+-}
+kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t
+kronecker a b = fromBlocks
+              . partit (cols a)
+              . map (reshape (cols b))
+              . toRows
+              $ flatten a `outer` flatten b