kronecker tools

5 files changed, 83 insertions, 19 deletions

diff --git a/CHANGES.md b/CHANGES.md
index 0ce7699..3b19ec9 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,7 +1,9 @@
 0.14.1.0
 --------
 
-- convolution
+- In Numeric.LinearAlgebra.Util:
+    convolution: corr, conv, corr2, conv2, separable, corrMin
+    kronecker: vec, vech, dup, vtrans
 
 0.14.0.0
 --------
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 0ef5cd7..21d6fa4 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -112,7 +112,6 @@ library
                         Numeric.LinearAlgebra.LAPACK,
                         Numeric.LinearAlgebra.Algorithms,
                         Numeric.LinearAlgebra.Util,
-                        Numeric.LinearAlgebra.Util.Convolution,
                         Graphics.Plot,
                         Data.Packed.ST,
                         Data.Packed.Development
@@ -129,7 +128,8 @@ library
                         Numeric.IO,
                         Numeric.Chain,
                         Numeric.Vector,
-                        Numeric.Matrix
+                        Numeric.Matrix,
+                        Numeric.LinearAlgebra.Util.Convolution
 
     C-sources:          lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
                         lib/Numeric/GSL/gsl-aux.c
diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs
index 79b8774..25eb239 100644
--- a/lib/Numeric/LinearAlgebra/Util.hs
+++ b/lib/Numeric/LinearAlgebra/Util.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE FlexibleContexts #-}
 -----------------------------------------------------------------------------
 {- |
 Module      :  Numeric.LinearAlgebra.Util
@@ -11,6 +12,7 @@ Stability   :  provisional
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Util(
+    -- * Convenience functions for real elements
     disp,
     zeros, ones,
     diagl,
@@ -19,11 +21,24 @@ module Numeric.LinearAlgebra.Util(
     (&),(!), (#),
     rand, randn,
     cross,
-    norm
+    norm,
+    -- * Convolution
+    -- ** 1D
+    corr, conv, corrMin,
+    -- ** 2D
+    corr2, conv2, separable,
+    -- * Tools for the Kronecker product
+    --
+    -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@
+    vec,
+    vech,
+    dup,
+    vtrans
 ) where
 
-import Numeric.LinearAlgebra
+import Numeric.LinearAlgebra hiding (i)
 import System.Random(randomIO)
+import Numeric.LinearAlgebra.Util.Convolution
 
 
 disp :: Int -> Matrix Double -> IO ()
@@ -87,7 +102,7 @@ col :: [Double] -> Matrix Double
 col = asColumn . fromList
 
 cross :: Vector Double -> Vector Double -> Vector Double
--- ^ cross product of dimension 3 real vectors
+-- ^ cross product (for three-element real vectors)
 cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
           | otherwise = error $ "cross ("++show x++") ("++show y++")"
   where
@@ -98,7 +113,40 @@ cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
     z3 = x1*y2-x2*y1
 
 norm :: Vector Double -> Double
--- ^ 2-norm of real vectors
+-- ^ 2-norm of real vector
 norm = pnorm PNorm2
 
+--------------------------------------------------------------------------------
+
+vec :: Element t => Matrix t -> Vector t
+-- ^ stacking of columns
+vec = flatten . trans
+
+
+vech :: Element t => Matrix t -> Vector t
+-- ^ half-vectorization (of the lower triangular part)
+vech m = join . zipWith f [0..] . toColumns $ m
+  where
+    f k v = subVector k (dim v - k) v
+
+
+dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t
+-- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k)
+dup k = trans $ fromRows $ map f es
+  where
+    rs = zip [0..] (toRows (ident (k^(2::Int))))
+    es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ]
+    f (i,j) | i == j = g (k*j + i)
+            | otherwise = g (k*j + i) + g (k*i + j)
+    g j = v
+      where
+        Just v = lookup j rs
+
+
+vtrans :: Element t => Int -> Matrix t -> Matrix t
+-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
+vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m
+           | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows"
+  where
+    (q,r) = divMod (rows m) p
 
diff --git a/lib/Numeric/LinearAlgebra/Util/Convolution.hs b/lib/Numeric/LinearAlgebra/Util/Convolution.hs
index b64b169..32cb188 100644
--- a/lib/Numeric/LinearAlgebra/Util/Convolution.hs
+++ b/lib/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -12,12 +12,8 @@ Stability   :  provisional
 -----------------------------------------------------------------------------
 
 module Numeric.LinearAlgebra.Util.Convolution(
-   -- * 1D
-   corr, conv,
-   -- * 2D
-   corr2, conv2,
-   -- * Misc
-   separable,  corrMin
+   corr, conv, corrMin,
+   corr2, conv2, separable
 ) where
 
 import Numeric.LinearAlgebra
@@ -32,8 +28,8 @@ corr :: Product t => Vector t -- ^ kernel
                   -> Vector t
 {- ^ correlation
 
-@\> (fromList[1,2,3]) (fromList [1..10])
-fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]@
+>>> corr (fromList[1,2,3]) (fromList [1..10])
+fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]
 
 -}
 corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker
@@ -43,8 +39,8 @@ corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker
 conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t
 {- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)
 
-@\> conv (fromList[1,1]) (fromList [-1,1])
-fromList [-1.0,0.0,1.0]@
+>>> conv (fromList[1,1]) (fromList [-1,1])
+fromList [-1.0,0.0,1.0]
 
 -}
 conv ker v = corr ker' v'
@@ -113,6 +109,6 @@ conv2 k m = corr2 (fliprl . flipud $ k) pm
 
 
 separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t
--- ^ 2D process implemented as separated 1D processes by rows and columns.
+-- ^ matrix computation implemented as separated vector operations by rows and columns.
 separable f = fromColumns . map f . toColumns . fromRows . map f . toRows
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 4e9a521..99c0c91 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -42,7 +42,7 @@ import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)
 import Control.Arrow((***))
 import Debug.Trace
 import Control.Monad(when)
-import Numeric.LinearAlgebra.Util.Convolution
+import Numeric.LinearAlgebra.Util hiding (ones,row,col)
 
 import Data.Packed.ST
 
@@ -435,6 +435,23 @@ convolutionTest = utest "convolution" ok
 
 --------------------------------------------------------------------------------
 
+kroneckerTest = utest "kronecker" ok
+  where
+    a,x,b :: Matrix Double
+    a = (3><4) [1..]
+    x = (4><2) [3,5..]
+    b = (2><5) [0,5..]
+    v1 = vec (a <> x <> b)
+    v2 = (trans b `kronecker` a) <> vec x
+    s = trans b <> b
+    v3 = vec s
+    v4 = dup 5 <> vech s
+    ok = v1 == v2 && v3 == v4
+      && vtrans 1 a == trans a
+      && vtrans (rows a) a == asColumn (vec a)
+
+--------------------------------------------------------------------------------
+
 
 
 -- | All tests must pass with a maximum dimension of about 20
@@ -611,6 +628,7 @@ runTests n = do
         , conformTest
         , accumTest
         , convolutionTest
+        , kroneckerTest
         ]
     when (errors c + failures c > 0) exitFailure
     return ()