diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -109,3 +109,7 @@
 
 0.2.5.4:
 		updated for hmatrix 0.18
+
+0.3:
+		changed PCA to use SVD as suggested by Pavol Klacansky 
+		(issue #3)
diff --git a/hstatistics.cabal b/hstatistics.cabal
--- a/hstatistics.cabal
+++ b/hstatistics.cabal
@@ -1,8 +1,8 @@
 Name:               hstatistics
-Version:            0.2.5.4
+Version:            0.3
 License:            BSD3
 License-file:       LICENSE
-Copyright:          (c) A.V.H. McPhail 2010, 2011, 2012, 2013, 2014, 2016
+Copyright:          (c) A.V.H. McPhail 2010--2014, 2016, 2017
 Author:             Vivian McPhail
 Maintainer:         haskell.vivian.mcphail <at> gmail <dot> com
 Stability:          provisional
@@ -30,8 +30,8 @@
     Build-Depends:      base >= 4 && < 5,
                         array, random,
                         vector,
-                        hmatrix >= 0.17,
-                        hmatrix-gsl-stats >= 0.4
+                        hmatrix >= 0.18,
+                        hmatrix-gsl-stats >= 0.4.1.6
 
     Extensions:         
 
@@ -46,12 +46,9 @@
     other-modules:      
     C-sources:          
 
-    ghc-prof-options:   -auto
-
     ghc-options:        -Wall -fno-warn-missing-signatures
                               -fno-warn-orphans
                               -fno-warn-unused-binds
-                        -O2
 
 source-repository head
     type:     git
diff --git a/lib/Numeric/Statistics/PCA.hs b/lib/Numeric/Statistics/PCA.hs
--- a/lib/Numeric/Statistics/PCA.hs
+++ b/lib/Numeric/Statistics/PCA.hs
@@ -1,7 +1,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Statistics.PCA
--- Copyright   :  (c) A. V. H. McPhail 2010, 2014
+-- Copyright   :  (c) A. V. H. McPhail 2010, 2014, 2017
 -- License     :  BSD3
 --
 -- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
@@ -27,38 +27,47 @@
 
 import Numeric.GSL.Statistics
 
-import Numeric.Statistics
+--import Numeric.Statistics
 
 -----------------------------------------------------------------------------
 
 -- | find the principal components of multidimensional data greater than
 --    the threshhold
-pca :: I.Array Int (Vector Double)    -- the data
+pca :: I.Array Int (Vector Double)     -- the data
     -> Double                         -- eigenvalue threshold
-    -> Matrix Double
+    -> (Vector Double, Matrix Double) -- Eignevalues, Principal components
 pca d q = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
-              cv = covarianceMatrix d'
-              (val',vec') = eigSH $ trustSym cv -- the covariance matrix is real symmetric
-              val = toList val'
-              vec = toColumns vec'
-              v' = zip val vec
+              d'' = fromColumns $ I.elems d'
+              (_,vec',uni') = svd d''
+              vec = toList vec'
+              uni = toColumns uni'
+              v' = zip vec uni
               v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
-          in fromColumns $ snd $ unzip v
+              (eigs,vs) = unzip v
+          in (fromList eigs,fromColumns vs) 
 
 -- | find N greatest principal components of multidimensional data
 --    according to size of the eigenvalue
 pcaN :: I.Array Int (Vector Double)    -- the data
      -> Int                            -- number of components to return
-     -> Matrix Double
+     -> (Vector Double, Matrix Double) -- Eignevalues, Principal components
 pcaN d n = let d' = fmap (\x -> x - (scalar $ mean x)) d -- remove the mean from each dimension
-               cv = covarianceMatrix d'
-               (val',vec') = eigSH $ trustSym cv  -- the covariance matrix is real symmetric
-               val = toList val'
-               vec = toColumns vec'
-               v' = zip val vec
+               d'' = fromColumns $ I.elems d'
+               (_,vec',uni') = svd d''
+               vec = toList vec'
+               uni = toColumns uni'
+               v' = zip vec uni
                v = take n $ reverse $ sortBy (comparing fst) v'
-           in fromColumns $ snd $ unzip v
+               (eigs,vs) = unzip v
+           in (fromList eigs,fromColumns vs) 
 
+v1 = fromList [1,2,3,4,5,6::Double]
+v2 = fromList [2,3,4,5,6,7::Double]
+v3 = fromList [3,4,5,6,7,8::Double]
+
+a = fromColumns [v1,v2,v3]
+b = I.listArray (1,3::Int) [v1,v2,v3] :: I.Array Int (Vector Double)
+                
 -- | perform a PCA transform of the original data (remove mean)
 -- |     Final = M^T Data^T
 pcaTransform :: I.Array Int (Vector Double)    -- ^ the data
@@ -74,11 +83,11 @@
           -> I.Array Int (Vector Double)      -- ^ the reduced data
 pcaReduce d q = let u = fmap (scalar . mean) d
                     d' = zipWith (-) (I.elems d) (I.elems u)
-                    cv = covarianceMatrix $ I.listArray (I.bounds d) d'
-                    (val',vec') = eigSH $ trustSym cv -- the covariance matrix is real symmetric
-                    val = toList val'
-                    vec = toColumns vec'
-                    v' = zip val vec
+                    d'' = fromColumns d'
+                    (_,vec',uni') = svd d''
+                    vec = toList vec'
+                    uni = toColumns uni'
+                    v' = zip vec uni
                     v = filter (\(x,_) -> x > q) v'  -- keep only eigens > than parameter
                     m = fromColumns $ snd $ unzip v
                  in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u) 
@@ -89,11 +98,11 @@
            -> I.Array Int (Vector Double)      -- ^ the reduced data, with n principal components
 pcaReduceN d n = let u = fmap (scalar . mean) d
                      d' = zipWith (-) (I.elems d) (I.elems u)
-                     cv = covarianceMatrix $ I.listArray (I.bounds d) d'
-                     (val',vec') = eigSH $ trustSym cv -- the covariance matrix is real symmetric
-                     val = toList val'
-                     vec = toColumns vec'
-                     v' = zip val vec
+                     d'' = fromColumns d'
+                     (_,vec',uni') = svd d''
+                     vec = toList vec'
+                     uni = toColumns uni'
+                     v' = zip vec uni
                      v = take n $ reverse $ sortBy (comparing fst) v'
                      m = fromColumns $ snd $ unzip v
                   in I.listArray (I.bounds d) $ zipWith (+) (toRows $ m <> (tr' m) <> fromRows d') (I.elems u)