Permalink
Browse files

Singular Value Decomposition working, and reconstruction

  • Loading branch information...
1 parent 695e50d commit 49fc700d08325b63cdbe8d769bcbe6fee5116e89 @markandrus committed Mar 10, 2012
Showing with 216 additions and 25 deletions.
  1. +1 −1 Makefile
  2. +73 −0 Options.hs
  3. +4 −4 PowerMethod.hs
  4. +43 −0 SVD.hs
  5. +16 −1 Utils.hs
  6. +75 −18 main.hs
  7. +3 −0 test.txt
  8. +1 −1 tests.hs
View
@@ -9,7 +9,7 @@ flags = -O2
extra_libs = -L/usr/lib
bins = main tests
-deps = Utils.hs
+deps = Options.hs PowerMethod.hs SVD.hs Utils.hs
all: $(bins)
View
@@ -0,0 +1,73 @@
+module Options
+ ( Options(..)
+ , defaultOptions
+ , options
+ , parseOptions
+ ) where
+
+import System.Environment
+import System.Exit
+import System.IO
+import System.Console.GetOpt
+
+data Options = Options
+ { optVerbose :: Bool
+ , optDecimalPlaces :: Int
+ , optEpsilon :: Double
+ , optDoPowerMethod :: Bool
+ , optDoSVD :: Bool
+ , optDoRecombine :: Bool
+ , optFilePath :: FilePath
+ }
+
+defaultOptions = Options
+ { optVerbose = False
+ , optDecimalPlaces = 3
+ , optEpsilon = 1.25e-3
+ , optDoPowerMethod = False
+ , optDoSVD = False
+ , optDoRecombine = False
+ , optFilePath = ""
+ }
+
+options :: [OptDescr (Options -> IO Options)]
+options =
+ [ Option "d" ["decimal-places"]
+ (ReqArg
+ (\arg opt -> return opt { optDecimalPlaces = read arg :: Int })
+ "INT")
+ "Number of decimal places to display"
+ , Option "e" ["epsilon"]
+ (ReqArg
+ (\arg opt -> return opt { optEpsilon = read arg :: Double })
+ "DOUBLE")
+ "Precision parameter `epsilon`"
+ , Option "P" ["power-method"]
+ (NoArg
+ (\opt -> return opt { optDoPowerMethod = True }))
+ "Compute the eigen-vectors and values of the 2D matrix"
+ , Option "S" ["svd"]
+ (NoArg
+ (\opt -> return opt { optDoSVD = True }))
+ "Compute the singular value decomposition of the 2D matrix"
+ , Option "R" ["recombine"]
+ (NoArg
+ (\opt -> return opt { optDoRecombine = True }))
+ "Compute the Frobenius norm of the difference between the 2D matrix and its reconstruction"
+ , Option "v" ["verbose"]
+ (NoArg
+ (\opt -> return opt { optVerbose = True }))
+ "Enable verbose messages"
+ , Option "h" ["help"]
+ (NoArg
+ (\_ -> do prg <- getProgName
+ hPutStrLn stderr (usageInfo prg options)
+ exitWith ExitSuccess))
+ "Show help"
+ ]
+
+parseOptions :: [String] -> IO Options
+parseOptions args =
+ let (actions, nonOptions, errors) = getOpt RequireOrder options args
+ filePath = if not (null nonOptions) then head nonOptions else "test.txt"
+ in foldl (>>=) (return defaultOptions { optFilePath = filePath }) actions
View
@@ -13,7 +13,7 @@ powerIterations
:: Int -> Matrix Double -> [(Vector Double, Double)]
powerIterations seed arr =
-- NOTE: that `rows arr == cols arr`
- let r = (fromList . take (rows arr) $ repeat 1) `add` (randomVector seed Uniform $ rows arr)
+ let r = (fromList . take (rows arr) $ repeat 1) `add` randomVector seed Uniform (rows arr)
-- Here we construct an infinite list of applications of
--
-- r \leftarrow \frac{ Ar } { ||Ar|| }
@@ -30,7 +30,7 @@ powerIterations seed arr =
-- |eigVecAndVal takes an Int `i` representing which eigenvalue we are looking for (`i` is used
-- to compute `lambda`), a precision parameter `epsilon`, a random seed for `powerIterations` and
-- an array; NOTE: the array should already have eigenvectors 1 through `i` zeroed out
-eigVecAndVal :: Int -> Int -> Double -> Matrix Double -> (Vector Double, Double)
+eigVecAndVal :: Int -> Int -> Epsilon -> Matrix Double -> (Vector Double, Double)
-- NOTE: since we implement no checks as to whether the matrix passed in is actually decomposible,
-- i.e., whether or not `r` will converge to the dominant eigenvector, `eigVecAndVal`
-- occasionally stalls calling `powerIterations`--at which point it may be useful to
@@ -57,13 +57,13 @@ eigVecAndVal seed i epsilon arr =
-- |powerMethod takes a gen to compute random vectors, a precision parameter `epsion`, and an NxN
-- matrix and computes the decomposition
--
--- Av_i = V \Lambda V^T = \lambda_i v_i
+-- A = V \Lambda V^T
--
-- powerMethod returns $(V, \Lambda)$
--
-- NOTE: the output resembles that of HMatrix's built-in `eigSH` for easy checking
powerMethod
- :: StdGen -> Double -> Matrix Double -> (Vector Double, Matrix Double)
+ :: StdGen -> Epsilon -> Matrix Double -> (Vector Double, Matrix Double)
powerMethod g epsilon arr =
-- NOTE: that `rows arr == cols arr`
let is = [1..cols arr]
View
@@ -0,0 +1,43 @@
+module SVD (mySVD) where
+
+import Numeric.LinearAlgebra
+import System.Random
+import Debug.Trace
+
+import PowerMethod
+import Utils
+
+-- |mySVD, given a random number generator, a precision parameter `epsilon`, and a 2D MxN matrix,
+-- computes a three-tuple containing the MxM matrix U, the RxR diagonal S (where R=min(M,N)), and
+-- the NxN matrix V such that the following holds
+--
+-- A = U S V^T
+--
+mySVD :: StdGen -> Epsilon -> Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
+mySVD g epsilon arr =
+ let aT = trans arr
+ aTa = aT <> arr
+ -- NOTE: that the eigenvalues of $A^T A$ and $A A^T$ are the same, therefore, we keep track
+ -- of just `eigVals`
+ (eigVals,bigV) = powerMethod g epsilon aTa
+ --
+ -- S = [ e_1 ... ... ...
+ -- , ... e_2 ... ...
+ -- , ... ... ... ...
+ -- , ... ... ... e_r ]
+ --
+ bigS = fromDiag eigVals
+ -- bigVT = trans bigV
+ -- `S` is a diagonal matrix, therefore
+ --
+ -- S^{-1} = [ 1/s_1 ... ... ...
+ -- , ... 1/s_2 ... ...
+ -- , ... ... ... ...
+ -- , ... ... ... 1/s_R ]
+ --
+ bigSInverse = fromDiag $ scaleRecip 1 eigVals
+ --
+ -- U = A V S^{-1}
+ --
+ bigU = arr <> bigV <> bigSInverse
+ in (bigU, eigVals, bigV)
View
@@ -1,13 +1,19 @@
{-# LANGUAGE FlexibleContexts #-}
module Utils
- ( isSymmetric
+ ( Epsilon
+ , isSymmetric
, vNormalize
, mNormalize
+ , fromDiag
+ , fromDecomposition
+ , fromSVD
) where
import Numeric.LinearAlgebra
+type Epsilon = Double
+
isSymmetric :: (Element a, Eq a) => Matrix a -> Bool
isSymmetric a = toLists a == toLists (trans a)
@@ -16,3 +22,12 @@ vNormalize v = recip (norm2 v) `scale` v
mNormalize :: Matrix Double -> Matrix Double
mNormalize m = recip (det m) `scale` m
+
+fromDiag :: Vector Double -> Matrix Double
+fromDiag v = diagRect 0 v (dim v) (dim v)
+
+fromDecomposition :: (Vector Double, Matrix Double) -> Matrix Double
+fromDecomposition (bigLambda, bigV) = bigV <> (fromDiag bigLambda) <> (trans bigV)
+
+fromSVD :: (Matrix Double, Vector Double, Matrix Double) -> Matrix Double
+fromSVD (bigU, bigS, bigV) = bigU <> (fromDiag bigS) <> (trans bigV)
View
@@ -1,28 +1,85 @@
module Main where
+import Control.Monad
import Data.List
import Numeric.LinearAlgebra
+import System.Environment
+import System.IO
import System.Random
+import Text.Printf
-import Utils
+import Options
import PowerMethod
+import SVD
+import Utils
main :: IO ()
main = do
- g_0 <- newStdGen
- let (seed,g) = next g_0
- epsilon = 0.000125
- n = 10
- arr_0 = reshape n $ randomVector seed Uniform (n*n)
- arr = arr_0 <> (trans arr_0)
- -- arr = ident n
- -- arr = reshape 3 $ fromList [1, 3, (-3), (-3), 7, (-3), (-6), 6, (-2)]
- (lambda, v) = powerMethod g epsilon (arr)
- putStrLn "powerMethod..."
- print arr
- print lambda
- print v
- putStrLn "eigSH..."
- let (lambda', v') = eigSH arr
- print lambda'
- print v'
+ -- Parse options
+ opts <- parseOptions =<< getArgs
+ let Options { optVerbose = verbose
+ , optDecimalPlaces = decimalPlaces
+ , optEpsilon = epsilon
+ , optDoPowerMethod = doPowerMethod
+ , optDoSVD = doSVD
+ , optDoRecombine = doRecombine
+ , optFilePath = filePath
+ } = opts
+ -- Setup matrix and vector display functins
+ mDisp = putStr . disps decimalPlaces
+ vDisp = putStr . vecdisp (disps decimalPlaces)
+ -- Create a random number generator
+ gen = mkStdGen 777
+ -- Read input matrix
+ when verbose (printf "Read 2D matrix from `%s'...\n\n\tA = ...\n\n" filePath)
+ arr <- loadMatrix filePath
+ when verbose (mDisp arr)
+ -- Decompose A using the Power Method
+ when doPowerMethod (do
+ when verbose (hPutStrLn stderr $
+ "\n\nComputing eigen-vectors and -values using the Power Method...\n\n"
+ ++ "\tA = V \\Lambda V^T")
+ -- Compute
+ let d@(bigLambda, bigV) = powerMethod gen epsilon arr
+ -- Print
+ putStrLn "\n# \\Lambda : Vector of eigenvalues"
+ vDisp bigLambda
+ putStrLn "\n# V : Matrix whose columns are eigenvectors"
+ mDisp bigV
+ -- Reconstruct A' from the decomposition of A
+ when doRecombine (do
+ when verbose (hPutStrLn stderr $
+ "\n\nComputing the Frobenius norm of the difference between A and "
+ ++ "A' = V \\Lambda V^T...")
+ -- Compute
+ let arr' = fromDecomposition d
+ -- Print
+ putStrLn "\n# A' : V \\Lambda V^T"
+ mDisp arr'
+ )
+ )
+ -- Do the Singular Value Decomposition
+ when doSVD (do
+ when verbose
+ (hPutStrLn stderr $ "\n\nComputing the Singular Value Decomposition...\n\n"
+ ++ "\tA = U S V^T")
+ -- Compute
+ let d@(bigU, bigS, bigV) = mySVD gen epsilon arr
+ -- Print
+ putStrLn "\n# U"
+ mDisp bigU
+ putStrLn "\n# S : Diagonal of singular values"
+ vDisp bigS
+ putStrLn "\n# V"
+ mDisp bigV
+ -- Reconstruct A' from the SVD of A
+ when doRecombine (do
+ when verbose (hPutStrLn stderr
+ "\n\nComputing the Frobenius norm of the difference vetween A and A' = U S V^T...")
+ -- Compute
+ let arr' = fromSVD d
+ -- Print
+ putStrLn "\n# A' : U S V^T"
+ mDisp arr'
+ )
+ )
View
@@ -0,0 +1,3 @@
+1 0 0
+0 1 0
+0 0 1
View
@@ -9,8 +9,8 @@ import System.Random
import Test.QuickCheck
import Text.Printf
-import Utils
import PowerMethod
+import Utils
compare :: (Eq a, Show a) => a -> a -> Property
compare ans ref = printTestCase message (ref==ans)

0 comments on commit 49fc700

Please sign in to comment.