spectral-clustering 0.3.1.0 → 0.3.1.1
raw patch · 3 files changed
+1/−416 lines, 3 filesdep −eigen
Dependencies removed: eigen
Files
- spectral-clustering.cabal +1/−4
- src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs +0/−213
- src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs +0/−199
spectral-clustering.cabal view
@@ -1,6 +1,6 @@ cabal-version: >=1.10 name: spectral-clustering-version: 0.3.1.0+version: 0.3.1.1 license: GPL-3 license-file: LICENSE copyright: 2019 Gregory W. Schwartz@@ -21,8 +21,6 @@ exposed-modules: Math.Clustering.Spectral.Dense Math.Clustering.Spectral.Sparse- Math.Clustering.Spectral.Eigen.FeatureMatrix- Math.Clustering.Spectral.Eigen.AdjacencyMatrix hs-source-dirs: src default-language: Haskell2010 ghc-options: -O2@@ -30,7 +28,6 @@ base >=4.7 && <5, containers >=0.5.11.0, clustering >=0.4.0,- eigen ==3.3.4.1, hmatrix >=0.19.0.0, hmatrix-svdlibc >=0.5.0.1, mwc-random >=0.13.6.0,
− src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs
@@ -1,213 +0,0 @@-{- Math.Clustering.Spectral.Eigen.AdjacencyMatrix-Gregory W. Schwartz--Collects the functions pertaining to spectral clustering.--}--{-# LANGUAGE BangPatterns #-}--module Math.Clustering.Spectral.Eigen.AdjacencyMatrix- ( spectralClusterKNorm- , spectralClusterNorm- , spectralNorm- , getDegreeMatrix- , secondLeft- , AdjacencyMatrix (..)- , LabelVector (..)- ) where---- Remote-import Control.Monad (replicateM)-import Data.Bool (bool)-import Data.Function (on)-import Data.List (sortBy, maximumBy, transpose)-import Data.Maybe (fromMaybe)-import Safe (headMay)-import System.Random.MWC (createSystemRandom, uniform)-import qualified AI.Clustering.KMeans as K-import qualified Data.Eigen.SparseMatrix as S-import qualified Data.Map.Strict as Map-import qualified Data.Vector as V-import qualified Data.Vector.Unboxed as U-import qualified Numeric.LinearAlgebra as H-import qualified Numeric.LinearAlgebra.Devel as H-import qualified Numeric.LinearAlgebra.SVD.SVDLIBC as SVD-import qualified Statistics.Quantile as Stat---- Local--type LabelVector = S.SparseMatrixXd-type AdjacencyMatrix = S.SparseMatrixXd---- | Assign close to 0 as 0.-epsilonZero :: Double -> Double-epsilonZero x = if abs x < 1e-12 then 0 else x---- | Returns the clustering of the eigenvectors with the second smallest--- eigenvalues of the symmetric normalized Laplacian L. Computes real symmetric--- part of L, so ensure the input is real and symmetric. Diagonal should be 0s--- for adjacency matrix. Clusters the eigenvector using kmeans into k groups.-spectralClusterKNorm :: Int -> Int -> AdjacencyMatrix -> LabelVector-spectralClusterKNorm e k mat- | S.rows mat < 1 = S.fromDenseList [[]]- | S.rows mat == 1 = S.fromDenseList [[0]]- | otherwise = kmeansVec k . spectralNorm 2 e $ mat---- | Returns the clustering of the eigenvectors with the second smallest--- eigenvalues of the symmetric normalized Laplacian L. Computes real symmetric--- part of L, so ensure the input is real and symmetric. Diagonal should be 0s--- for adjacency matrix. Clusters the eigenvector by sign.-spectralClusterNorm :: AdjacencyMatrix -> LabelVector-spectralClusterNorm mat- | S.rows mat < 1 = S.fromDenseList [[]]- | S.rows mat == 1 = S.fromDenseList [[0]]- | otherwise = S.fromDenseList- . (fmap . fmap) (bool 0 1 . (>= 0))- . S.toDenseList- . spectralNorm 2 1- $ mat---- | Returns the eigenvector with the second smallest eigenvalue (or N start)--- and E on of the symmetric normalized Laplacian L. Computes real symmetric--- part of L, so ensure the input is real and symmetric. Diagonal should be 0s--- for adjacency matrix. Uses I + Lnorm instead of I - Lnorm to find second--- largest singular value instead of second smallest for Lnorm.-spectralNorm :: Int -> Int -> AdjacencyMatrix -> S.SparseMatrixXd-spectralNorm n e mat- | e < 1 = error "Less than 1 eigenvector chosen for clustering."- | n < 1 = error "N < 1, cannot go before first eigenvector."- | otherwise = S._map epsilonZero $ secondLeft n e lNorm- where- lNorm = i + (S.transpose invRootD * (mat * invRootD))- invRootD = S.diagRow 0- . S._map (\x -> if x == 0 then x else x ** (- 1 / 2))- . getDegreeVector- $ mat- i = S.ident . S.rows $ mat---- | Second largest eigenvector. Unused, untested, don't use.-secondLargest :: S.SparseMatrixXd -> IO S.SparseMatrixXd-secondLargest a = do- (first, firstVal) <- powerIt a-- let firstScale = S.scale firstVal (first * S.transpose first)- b = a - firstScale-- fmap fst $ powerIt b---- | Rayleigh quotient. Takes a matrix and an eigenvector to find the--- corresponding eigenvalue. Unused, untested, don't use.-rayQuot :: S.SparseMatrixXd -> S.SparseMatrixXd -> Double-rayQuot a x = ((S.transpose (a * x) * x) S.! (0, 0))- / ((S.transpose x * x) S.! (0, 0))---- | Power iteration. Unused, untested, don't use.-powerIt :: S.SparseMatrixXd -> IO (S.SparseMatrixXd, Double)-powerIt a = do- g <- createSystemRandom- start <-- fmap (S.fromDenseList . fmap (: [])) . replicateM (S.rows a) . uniform $ g- let go- :: Int- -> Double- -> S.SparseMatrixXd -- Eigenvector guess.- -> Double -- Eigenvalue guess.- -> (S.SparseMatrixXd, Double)- go !i !e !b !lambda =- if (abs (lambda' - lambda) < e) || (i < 0)- then (b', lambda')- else go (i - 1) e b' lambda'- where- absMat :: S.SparseMatrixXd -> S.SparseMatrixXd- absMat = S._map abs- b' :: S.SparseMatrixXd- b' = S.scale (1 / maxElement ab) ab- ab :: S.SparseMatrixXd- ab = a * b- maxElement = maximum . fmap (\(_, _, !x) -> abs x) . S.toList- lambda' = S.norm ab / S.norm b- return $ go 1000 0.000001 start (1 / 0)---- | Executes kmeans to cluster a one dimensional vector.-kmeansVec :: Int -> S.SparseMatrixXd -> LabelVector-kmeansVec k = consensusKmeans 100- . V.fromList- . fmap U.fromList- . concatMap S.toDenseList- . S.getRows- . S.fromCols- . fmap normNormalize- . S.getCols---- | Consensus kmeans.-consensusKmeans :: Int -> V.Vector (U.Vector Double) -> LabelVector-consensusKmeans x vs = S.fromDenseList- . fmap ((:[]) . fromIntegral . mostCommon)- . transpose- . fmap kmeansFunc- $ [1 .. fromIntegral x]- where- kmeansFunc run =- (\xs -> if headMay xs == Just 1 then fmap (bool 0 1 . (== 0)) xs else xs)- . U.toList- . K.membership- . K.kmeansBy 2 vs id- $ K.defaultKMeansOpts- { K.kmeansMethod = K.Forgy- , K.kmeansClusters = False- , K.kmeansSeed = U.fromList [run]- }---- | Get the most common element of a list.-mostCommon :: (Ord a) => [a] -> a-mostCommon [] = error "Cannot find most common element of empty list."-mostCommon [x] = x-mostCommon xs = fst- . maximumBy (compare `on` snd)- . Map.toAscList- . Map.fromListWith (+)- . flip zip [1,1..]- $ xs---- | Normalize by the norm of a vector.-normNormalize :: S.SparseMatrixXd -> S.SparseMatrixXd-normNormalize xs = S._map (/ norm) xs- where- norm = S.norm xs---- | Obtain the second largest value singular vector (or Nth) and E on of a--- sparse matrix.-secondLeft :: Int -> Int -> S.SparseMatrixXd -> S.SparseMatrixXd-secondLeft n e m = S.transpose- . S.fromDenseList- . fmap H.toList- . drop (n - 1)- . H.toRows- . (\(!x, _, _) -> x)- . SVD.sparseSvd (e + (n - 1))- . H.mkCSR- . fmap (\(!i, !j, !x) -> ((i, j), x))- . S.toList- $ m---- | Obtain the second largest value singular vector of a sparse matrix.-denseSecondLeft :: S.SparseMatrixXd -> S.SparseMatrixXd-denseSecondLeft m = S.fromDenseList- . fmap (:[])- . H.toList- . (!! 2)- . H.toColumns- . (\(!x, _, _) -> x)- . H.svd- . H.assoc (S.rows m, S.cols m) 0- . fmap (\(!i, !j, !x) -> ((i, j), x))- . S.toList- $ m---- | Obtain the signed degree matrix. Faster for columns.-getDegreeMatrix :: AdjacencyMatrix -> S.SparseMatrixXd-getDegreeMatrix = S.diagRow 0 . getDegreeVector---- | Obtain the signed degree vector. Faster for columns.-getDegreeVector :: AdjacencyMatrix -> S.SparseMatrixXd-getDegreeVector = S.getColSums . S._map abs
− src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs
@@ -1,199 +0,0 @@-{- Math.Clustering.Spectral.Eigen.FeatureMatrix-Gregory W. Schwartz--Collects the functions pertaining to sparse spectral clustering.--}--{-# LANGUAGE BangPatterns #-}--module Math.Clustering.Spectral.Eigen.FeatureMatrix- ( B (..)- , B1 (..)- , B2 (..)- , LabelVector (..)- , spectral- , spectralCluster- , spectralClusterK- , getB- , b1ToB2- , getSimilarityFromB2- ) where---- Remote-import Data.Bool (bool)-import Data.Function (on)-import Data.List (sortBy, maximumBy, transpose)-import Data.Maybe (fromMaybe)-import Safe (headMay)-import qualified AI.Clustering.KMeans as K-import qualified Data.Eigen.SparseMatrix as S-import qualified Data.Map.Strict as Map-import qualified Data.Vector as V-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Unboxed as U-import qualified Numeric.LinearAlgebra as H-import qualified Numeric.LinearAlgebra.Devel as H-import qualified Numeric.LinearAlgebra.SVD.SVDLIBC as SVD---- Local---- | Output vector containing cluster assignment (0 or 1).-type LabelVector = S.SparseMatrixXd--- | B1 observation by feature matrix.-newtype B1 = B1 { unB1 :: S.SparseMatrixXd } deriving (Show)--- | B2 term frequency-inverse document frequency matrix of B1.-newtype B2 = B2 { unB2 :: S.SparseMatrixXd } deriving (Show)--- | Diagonal matrix from \(diag(B(B^{T}1))\).-newtype D = D { unD :: S.SparseMatrixXd } deriving (Show)--- | Matrix from \(D^{-1/2}B}\).-newtype C = C { unC :: S.SparseMatrixXd } deriving (Show)--- | Normed rows of B2. For a complete explanation, see Shu et al., "Efficient--- Spectral Neighborhood Blocking for Entity Resolution", 2011.-newtype B = B { unB :: S.SparseMatrixXd } deriving (Show)---- | Assign close to 0 as 0.-epsilonZero :: Double -> Double-epsilonZero x = if abs x < 1e-12 then 0 else x---- | Normalize the input matrix by column. Here, columns are features.-b1ToB2 :: B1 -> B2-b1ToB2 (B1 b1) =- B2- . S._imap (\ i j x- -> (log (fromIntegral n / (fromMaybe 0 $ dVec VS.!? j))) * x- )- $ b1- where- dVec :: VS.Vector Double- dVec = maybe (error "Cannot get number of non-zeros.") (VS.map fromIntegral)- . S.innerNNZs- . S.uncompress- $ b1- n = S.rows b1---- | Euclidean norm each row.-b2ToB :: B2 -> B-b2ToB (B2 b2) =- B- . S._imap (\ i j x- -> x / (fromMaybe (error "Norm is 0.") $ eVec VS.!? i)- )- $ b2- where- eVec :: VS.Vector Double- eVec = VS.fromList . fmap S.norm . S.getRows $ b2---- | Get the signed diagonal transformed B matrix.-bToD :: B -> D-bToD (B b) = D . S.diagCol 0 $ (S._map abs b) * ((S._map abs $ S.transpose b) * S.ones n)- where- n = S.rows b---- | Get the matrix C as input for SVD.-bdToC :: B -> D -> C-bdToC (B b) (D d) = C $ (S._map (\x -> x ** (- 1 / 2)) d) * b---- | Obtain the second largest value singular vector (or Nth) and E on of a--- sparse matrix.-secondLeft :: Int -> Int -> S.SparseMatrixXd -> S.SparseMatrixXd-secondLeft n e m = S.transpose- . S.fromDenseList- . fmap H.toList- . drop (n - 1)- . H.toRows- . (\(!x, _, _) -> x)- . SVD.sparseSvd (e + (n - 1))- . H.mkCSR- . fmap (\(!i, !j, !x) -> ((i, j), x))- . S.toList- $ m---- | Get the normalized matrix B from an input matrix where the features are--- columns and rows are observations. Optionally, do not normalize.-getB :: Bool -> S.SparseMatrixXd -> B-getB True = b2ToB . b1ToB2 . B1-getB False = b2ToB . B2---- | Returns the second left singular vector (or Nth) of a sparse spectral--- process. Assumes the columns are features and rows are observations. B is the--- normalized matrix (from getB). See Shu et al., "Efficient Spectral--- Neighborhood Blocking for Entity Resolution", 2011.-spectral :: Int -> Int -> B -> S.SparseMatrixXd-spectral n e b = S._map epsilonZero . secondLeft n e . unC . bdToC b . bToD $ b---- | Returns a vector of cluster labels for two groups by finding the second--- left singular vector of a special normalized matrix. Assumes the columns are--- features and rows are observations. B is the normalized matrix (from getB).--- See Shu et al., "Efficient Spectral Neighborhood Blocking for Entity--- Resolution", 2011.-spectralCluster :: B -> LabelVector-spectralCluster (B b)- | S.rows b < 1 = S.fromDenseList [[]]- | S.rows b == 1 = S.fromDenseList [[0]]- | otherwise = S.fromDenseList- . (fmap . fmap) (bool 0 1 . (>= 0))- . S.toDenseList- . spectral 2 1- $ B b---- | Returns a vector of cluster labels for two groups by finding the largest--- singular vectors and on of a special normalized matrix and running kmeans.--- Assumes the columns are features and rows are observations. B is the--- normalized matrix (from getB). See Shu et al., "Efficient Spectral--- Neighborhood Blocking for Entity Resolution", 2011.-spectralClusterK :: Int -> Int -> B -> LabelVector-spectralClusterK e k (B b)- | S.rows b < 1 = S.fromDenseList [[]]- | S.rows b == 1 = S.fromDenseList [[0]]- | otherwise = consensusKmeans 100- . V.fromList- . fmap U.fromList- . concatMap S.toDenseList- . S.getRows- . S.fromCols- . fmap normNormalize- . S.getCols- . spectral 2 e- $ B b---- | Consensus kmeans.-consensusKmeans :: Int -> V.Vector (U.Vector Double) -> LabelVector-consensusKmeans x vs = S.fromDenseList- . fmap ((:[]) . fromIntegral . mostCommon)- . transpose- . fmap kmeansFunc- $ [1 .. fromIntegral x]- where- kmeansFunc run =- (\xs -> if headMay xs == Just 1 then fmap (bool 0 1 . (== 0)) xs else xs)- . U.toList- . K.membership- . K.kmeansBy 2 vs id- $ K.defaultKMeansOpts- { K.kmeansMethod = K.Forgy- , K.kmeansClusters = False- , K.kmeansSeed = U.fromList [run]- }---- | Get the most common element of a list.-mostCommon :: (Ord a) => [a] -> a-mostCommon [] = error "Cannot find most common element of empty list."-mostCommon [x] = x-mostCommon xs = fst- . maximumBy (compare `on` snd)- . Map.toAscList- . Map.fromListWith (+)- . flip zip [1,1..]- $ xs---- | Normalize by the norm of a vector.-normNormalize :: S.SparseMatrixXd -> S.SparseMatrixXd-normNormalize xs = S._map (/ norm) xs- where- norm = S.norm xs---- | Get the cosine similarity between two rows using B2.-getSimilarityFromB2 :: B2 -> Int -> Int -> Double-getSimilarityFromB2 (B2 b2) i j =- (((S.getRow i b2) * (S.transpose $ S.getRow j b2)) S.! (0, 0))- / (S.norm (S.getRow i b2) * S.norm (S.getRow j b2))