hmatrix-nipals 0.1 → 0.2

raw patch · 3 files changed

+98/−27 lines, 3 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• hmatrix-nipals.cabal +16/−10
• src/Numeric/LinearAlgebra/NIPALS.hs +68/−17
• test/tests.hs +14/−0

@@ -1,5 +1,5 @@ Name:                hmatrix-nipals-Version:             0.1+Version:             0.2 Synopsis:   NIPALS method for Principal Components Analysis on large data-sets. @@ -26,16 +26,17 @@   .   NIPALS is not generally recommended because sample matrices where   the largest eigenvalues are close in magnitude will cause NIPALS to-  converge very slowly. In general, Lanczos methods-  <http://en.wikipedia.org/wiki/Lanczos_algorithm> or some other-  truncated singular value decomposition algorithm are preferred to-  NIPALS because of this convergence issue, but these methods often-  require the sample matrix to fit in memory, or store large-  conditioning matrices, which isn't always feasible. However, if you-  know of free and memory-efficient implementations of these more-  sophisticated algorithms, please contact the author with a pointer.+  converge very slowly. For sparse matrices, use Lanczos methods+  <http://en.wikipedia.org/wiki/Lanczos_algorithm>, and for dense+  matrices, random-projection methods+  <http://amath.colorado.edu/faculty/martinss/Pubs/2009_HMT_random_review.pdf>+  can be used. However, these methods are harder to implement in a+  single pass. If you know of a good, single-pass, and+  memory-efficient implementation of either of these methods, please+  contact the author.  Homepage:            http://github.com/alanfalloon/hmatrix-nipals+Bug-reports:         https://github.com/alanfalloon/hmatrix-nipals/issues License:             LGPL-2.1 License-file:        LICENSE Author:              Alan Falloon@@ -44,7 +45,12 @@ Stability:           Experimental Category:            Math Build-type:          Simple-Cabal-version:       >=1.4+Cabal-version:       >=1.6++Source-repository head+  Type:              git+  Location:          git://github.com/alanfalloon/hmatrix-nipals.git+  Branch:            master  Flag test   Description:       Build unit-tests

@@ -4,8 +4,12 @@     ( -- * Simplified Interface       firstPC     , firstPCFromScores+      -- * Monadic interface+    , firstPCFromScoresM+    , residual     ) where +import Data.List (foldl1') import Numeric.LinearAlgebra  -- | Calculate the first principal component of a set of samples.@@ -48,7 +52,7 @@       steps = iterate refine (t0,undefined)       convergence = let scores = map fst steps                         dscores = zipWith diffScores scores $ tail scores-                    in smooth dscores+                    in dscores       (t,p) = let steps' = zip convergence $ tail steps                   steps'' = dropWhile (\(c,_) -> c > threshold) steps'               in snd $ head steps''@@ -59,29 +63,76 @@             p' = toUnit (trans m <> t)             t' = m <> p' -      diffScores ta tb = let xa = ta<.>ta-                             xb = tb<.>tb-                             x = abs $ xa - xb-                         in (x / xb) :: Double+      threshold = 16*eps -      threshold = sqrt (fromIntegral (cols m + rows m)) * eps+diffScores :: Vector Double -> Vector Double -> Double+diffScores ta tb = let xa = ta<.>ta+                       xb = tb<.>tb+                       x = abs $ xa - xb+                   in (x / (eps+xb)) :: Double  toUnit :: Vector Double -> Vector Double toUnit v = if mag <= 0.0 then dim v |> (1 : repeat 0) else scale (1/mag) v     where       mag = norm2 v -smooth :: [Double] -> [Double]-smooth (x0:xs@(x1:_)) = (2*x0+x1)/3 : smooth xs-smooth _ = []- -- | Calculate the first principal component -- calculating the -- samples fresh on every pass. ----- This function calculates the exact same results as 'firstPC', but--- instead of an input 'Matrix', it takes a monad action that yields--- the list of samples, and it guarantees that the list returned by--- the action will be consumed in a single pass. However the action--- may be demanded many times.-firstPCM :: (Product t, Monad m) => m [Vector t] -> m (Vector t, Vector t, [Vector t])-firstPCM = undefined+-- This function calculates the exact same results as+-- 'firstPCFromScores' (minus the residual), but instead of an input+-- 'Matrix', it takes a monad action that yields the list of samples,+-- and it guarantees that the list returned by the action will be+-- consumed in a single pass. However the action may be demanded many+-- times.+--+-- The residual can't be calculated lazily, like it is in+-- 'firstPCFromScores', because the samples would need to be+-- demanded. Instead, to calculate the residual use 'residual'.+--+-- There is no corresponding @firstPCM@ that guesses the initial score+-- vector for you because if you need to use this function instead of+-- 'firstPC', then you really should come up with a reasonable+-- starting point or it will take forever.+firstPCFromScoresM :: Monad m =>+                      m [Vector Double]+                   -> Vector Double+                   -> m (Vector Double, Vector Double)+firstPCFromScoresM samplesM t0 = do+  (t,p) <- refine t0+  iter t0 t p+      where+        refine t = do+          p <- samplesM >>= (\s -> return $ transMult1Pass s t)+          let p' = toUnit p+          t' <- samplesM >>= (\s -> return $ mult1Pass s p')+          return (t',p')+        threshold = 16 * eps+        iter t0 t1 p1 | diff < threshold = return $ (p1,t1)+                      | otherwise = do+                                   (t2,p2) <- refine t1+                                   iter t1 t2 p2+                     where+                       diff = diffScores t0 t1++-- Single pass versions of (trans m <> v) and (m <> v) respectively,+-- where m is a list of rows+transMult1Pass, mult1Pass :: [Vector Double] -> Vector Double -> Vector Double+mult1Pass rows vec = fromList $ map (dot vec) rows+transMult1Pass cols vec = foldl1' add $ zipWith scale (toList vec) cols++-- | Calculate the residuals of a series of samples given a component+-- and score vector.+--+-- > (p,t) <- firstPCFromScoresM samplesM (randomVector 0 Gaussian numSamples)+-- > samples <- samplesM+-- > let r = residual samples p t+residual :: [Vector Double] -- ^ The samples+         -> Vector Double   -- ^ The component (also called the loading)+         -> Vector Double   -- ^ The scores+         -> [Vector Double] -- ^ The residuals for each sample+residual s p t = r+    where+      ts = toList t+      comp = map (`scale` p) ts+      r = zipWith sub s comp

@@ -1,4 +1,5 @@ import Data.List+import Data.Maybe (fromJust) import Numeric.LinearAlgebra import Numeric.LinearAlgebra.NIPALS import Foreign.Storable@@ -15,6 +16,7 @@       ]     , testGroup "Correctness"       [ testProperty "resultCanRecoverInput" prop_recoverInput+      , testProperty "monadSameAnswer" prop_monadSameAnswer       ]     ] @@ -43,6 +45,18 @@       tRelErr = relativeDifference t t'       pRelErr = relativeDifference p p'       th = sqrt $ fromIntegral (cols m * rows m) * eps++prop_monadSameAnswer m = tRelErr <= th .&&. pRelErr <= th .&&. rRelErr <= th+    where+      guess = head $ toColumns m+      samplesM = return $ toRows m+      (p,t,r) = firstPCFromScores m guess+      (p',t') = fromJust $ firstPCFromScoresM samplesM guess+      r' = fromRows $ residual (toRows m) p' t'+      tRelErr = relativeDifference t t'+      pRelErr = relativeDifference p p'+      rRelErr = if rank m > 2 then relativeDifference r r' else 0+      th = 1e-7  relativeDifference :: (Normed c t, Container c t) => c t -> c t -> Double relativeDifference x y = realToFrac (norm (x `sub` y) / (peps + norm x + norm y))