diff --git a/hmatrix-nipals.cabal b/hmatrix-nipals.cabal --- a/hmatrix-nipals.cabal +++ b/hmatrix-nipals.cabal @@ -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 - 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 + , and for dense + matrices, random-projection methods + + 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 diff --git a/src/Numeric/LinearAlgebra/NIPALS.hs b/src/Numeric/LinearAlgebra/NIPALS.hs --- a/src/Numeric/LinearAlgebra/NIPALS.hs +++ b/src/Numeric/LinearAlgebra/NIPALS.hs @@ -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 diff --git a/test/tests.hs b/test/tests.hs --- a/test/tests.hs +++ b/test/tests.hs @@ -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))