diff --git a/spectral-clustering.cabal b/spectral-clustering.cabal --- a/spectral-clustering.cabal +++ b/spectral-clustering.cabal @@ -1,6 +1,6 @@ cabal-version: >=1.10 name: spectral-clustering -version: 0.2.2.3 +version: 0.3.0.0 license: GPL-3 license-file: LICENSE copyright: 2019 Gregory W. Schwartz @@ -28,6 +28,7 @@ ghc-options: -O2 build-depends: base >=4.7 && <5, + containers >=0.5.11.0, clustering >=0.4.0, eigen ==3.3.4.1, hmatrix >=0.19.0.0, diff --git a/src/Math/Clustering/Spectral/Dense.hs b/src/Math/Clustering/Spectral/Dense.hs --- a/src/Math/Clustering/Spectral/Dense.hs +++ b/src/Math/Clustering/Spectral/Dense.hs @@ -4,6 +4,7 @@ Collects the functions pertaining to spectral clustering. -} + {-# LANGUAGE BangPatterns #-} module Math.Clustering.Spectral.Dense @@ -27,10 +28,11 @@ -- Remote import Data.Bool (bool) import Data.Function (on) -import Data.List (sortBy) +import Data.List (sortBy, maximumBy, transpose) import Data.Maybe (fromMaybe) import Safe (headMay) import qualified AI.Clustering.KMeans as K +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 @@ -59,6 +61,10 @@ -- Spectral Neighborhood Blocking for Entity Resolution", 2011. newtype B = B { unB :: H.Matrix Double } deriving (Show) +-- | Assign close to 0 as 0. +epsilonZero :: Double -> Double +epsilonZero x = if abs x < 1e-12 then 0 else x + -- | Map hmatrix with indices. cimap :: (Int -> Int -> Double -> Double) -> H.Matrix Double -> H.Matrix Double cimap f mat = H.assoc (H.size mat) 0 @@ -144,7 +150,8 @@ spectral n e b | e < 1 = error "Less than 1 eigenvector chosen for clustering." | n < 1 = error "N < 1, cannot go before first eigenvector." - | otherwise = secondLeft n e . unC . bdToC b . bToD $ b + | otherwise = + fmap (H.cmap 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 @@ -169,22 +176,50 @@ spectralClusterK e k (B b) | H.rows b < 1 = H.fromList [] | H.rows b == 1 = H.fromList [0] - | otherwise = kmeansVec k . spectral 1 e $ B b + | otherwise = kmeansVec k . spectral 2 e $ B b -- | Executes kmeans to cluster a vector. kmeansVec :: Int -> [H.Vector Double] -> LabelVector -kmeansVec k = V.convert - . U.map fromIntegral - . K.membership - . (\x -> K.kmeansBy k x id K.defaultKMeansOpts) +kmeansVec k = consensusKmeans 100 . V.fromList - . fmap V.convert + . fmap U.convert . H.toRows . H.fromColumns . fmap H.normalize -- Normalize within eigenvectors (columns). . H.toColumns . H.fromRows + where +-- | Consensus kmeans. +consensusKmeans :: Int -> V.Vector (U.Vector Double) -> LabelVector +consensusKmeans x vs = H.fromList + . 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 + -- | Get the cosine similarity between two rows using B2. getSimilarityFromB2 :: B2 -> Int -> Int -> Double getSimilarityFromB2 (B2 b2) i j = @@ -201,7 +236,7 @@ | H.rows mat < 1 = H.fromList [] | H.rows mat == 1 = H.fromList [0] | otherwise = kmeansVec k - . spectralNorm 1 e + . spectralNorm 2 e $ mat -- | Returns the eigenvector with the second smallest eigenvalue of the @@ -212,8 +247,10 @@ spectralClusterNorm mat | H.rows mat < 1 = H.fromList [] | H.rows mat == 1 = H.fromList [0] - | otherwise = - H.cmap (bool 0 1 . (>= 0)) . mconcat . spectralNorm 2 1 $ mat + | otherwise = H.cmap (bool 0 1 . (>= 0)) + . mconcat + . spectralNorm 2 1 + $ mat -- | Returns the eigenvectors with the Nth smallest eigenvalue and on of the -- symmetric normalized Laplacian L. Computes real symmetric part of L, so @@ -224,18 +261,26 @@ | e < 1 = error "Less than 1 eigenvector chosen for clustering." | n < 1 = error "N < 1, cannot go before first eigenvector." | otherwise = H.toRows + . H.cmap epsilonZero -- Correct accuracy. . flip (H.??) (H.All, H.TakeLast e) . flip (H.??) (H.All, H.DropLast (n - 1)) . snd . H.eigSH $ lNorm where - lNorm = H.sym $ i - mconcat [invD, mat, invD] + lNorm = H.trustSym $ i - mconcat [invD, mat, invD] invD = H.diag . H.cmap (\x -> if x == 0 then x else x ** (- 1 / 2)) . getDegreeVector $ mat i = H.ident . H.rows $ mat + +-- | Sort an matrix by values in a vector. +sortMatrixByVec :: H.Vector Double -> H.Matrix Double -> H.Matrix Double +sortMatrixByVec xs mat = mat H.¿ sortedIdx + where + sortedIdx = + reverse . fmap fst . sortBy (compare `on` snd) . zip [0..] . H.toList $ xs -- | Obtain the signed degree matrix. getDegreeMatrix :: AdjacencyMatrix -> H.Matrix Double diff --git a/src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs b/src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs --- a/src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs +++ b/src/Math/Clustering/Spectral/Eigen/AdjacencyMatrix.hs @@ -20,12 +20,13 @@ import Control.Monad (replicateM) import Data.Bool (bool) import Data.Function (on) -import Data.List (sortBy) +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 @@ -38,6 +39,10 @@ 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 @@ -46,7 +51,7 @@ spectralClusterKNorm e k mat | S.rows mat < 1 = S.fromDenseList [[]] | S.rows mat == 1 = S.fromDenseList [[0]] - | otherwise = kmeansVec k . spectralNorm 1 e $ mat + | 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 @@ -71,7 +76,7 @@ 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 = secondLeft n e lNorm + | otherwise = S._map epsilonZero $ secondLeft n e lNorm where lNorm = i + (S.transpose invRootD * (mat * invRootD)) invRootD = S.diagRow 0 @@ -125,11 +130,7 @@ -- | Executes kmeans to cluster a one dimensional vector. kmeansVec :: Int -> S.SparseMatrixXd -> LabelVector -kmeansVec k = S.fromDenseList - . fmap ((:[]) . fromIntegral) - . U.toList - . K.membership - . (\x -> K.kmeansBy k x id K.defaultKMeansOpts) +kmeansVec k = consensusKmeans 100 . V.fromList . fmap U.fromList . concatMap S.toDenseList @@ -137,6 +138,36 @@ . 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 diff --git a/src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs b/src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs --- a/src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs +++ b/src/Math/Clustering/Spectral/Eigen/FeatureMatrix.hs @@ -22,10 +22,12 @@ -- Remote import Data.Bool (bool) import Data.Function (on) -import Data.List (sortBy) +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 @@ -49,6 +51,10 @@ -- 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) = @@ -113,7 +119,7 @@ -- 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 = secondLeft n e . unC . bdToC b . bToD $ b +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 @@ -139,11 +145,7 @@ spectralClusterK e k (B b) | S.rows b < 1 = S.fromDenseList [[]] | S.rows b == 1 = S.fromDenseList [[0]] - | otherwise = S.fromDenseList - . fmap ((:[]) . fromIntegral) - . U.toList - . K.membership - . (\x -> K.kmeansBy k x id K.defaultKMeansOpts) + | otherwise = consensusKmeans 100 . V.fromList . fmap U.fromList . concatMap S.toDenseList @@ -151,8 +153,38 @@ . S.fromCols . fmap normNormalize . S.getCols - . spectral 1 e + . 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 diff --git a/src/Math/Clustering/Spectral/Sparse.hs b/src/Math/Clustering/Spectral/Sparse.hs --- a/src/Math/Clustering/Spectral/Sparse.hs +++ b/src/Math/Clustering/Spectral/Sparse.hs @@ -27,15 +27,17 @@ import Data.Bool (bool) import Data.Maybe (fromMaybe) import Data.Function (on) -import Data.List (sortBy, foldl1') +import Data.List (sortBy, foldl1', maximumBy, transpose) import Safe (headMay) import qualified AI.Clustering.KMeans as K +import qualified Data.Map.Strict as Map import qualified Data.Sparse.Common as S 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 Debug.Trace -- Local @@ -115,8 +117,10 @@ secondLeft n e m = fmap (S.sparsifySV . S.fromListDenseSV e . drop (n - 1) . H.toList) . H.toColumns + . (\x -> traceShow x x) . (\(!x, _, _) -> x) . SVD.sparseSvd (e + (n - 1)) + . (\x -> traceShow x x) . H.mkCSR . fmap (\(!i, !j, !x) -> ((i, j), x)) . S.toListSM @@ -136,7 +140,8 @@ spectral n e b | e < 1 = error "Less than 1 eigenvector chosen for clustering." | n < 1 = error "N < 1, cannot go before first eigenvector." - | otherwise = secondLeft n e . unC . bdToC b . bToD $ b + | otherwise = + fmap S.sparsifySV . 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 @@ -164,16 +169,11 @@ spectralClusterK e k (B b) | S.nrows b < 1 = S.zeroSV 0 | S.nrows b == 1 = S.zeroSV 1 - | otherwise = kmeansVec k . spectral 1 e $ B b + | otherwise = kmeansVec k . spectral 2 e $ B b -- | Executes kmeans to cluster a vector. kmeansVec :: Int -> [S.SpVector Double] -> LabelVector -kmeansVec k = S.sparsifySV - . S.vr - . fmap fromIntegral - . U.toList - . K.membership - . (\x -> K.kmeansBy k x id K.defaultKMeansOpts) +kmeansVec k = consensusKmeans 100 . V.fromList . fmap (U.fromList . S.toDenseListSV) . S.toRowsL @@ -183,6 +183,37 @@ . S.transpose . S.fromColsL +-- | Consensus kmeans. +consensusKmeans :: Int -> V.Vector (U.Vector Double) -> LabelVector +consensusKmeans x vs = S.sparsifySV + . S.vr + . 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 + -- | Get the cosine similarity between two rows using B2. getSimilarityFromB2 :: B2 -> Int -> Int -> Double getSimilarityFromB2 (B2 b2) i j = @@ -195,7 +226,7 @@ -- 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.SpVector Double] -spectralNorm n e mat = secondLeft n e lNorm +spectralNorm n e mat = fmap S.sparsifySV $ secondLeft n e lNorm where lNorm = i S.^+^ (S.transpose invRootD S.#~# (mat S.#~# invRootD)) invRootD = S.diagonalSM @@ -211,7 +242,7 @@ -- 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 = kmeansVec k . spectralNorm 1 e +spectralClusterKNorm e k = kmeansVec k . spectralNorm 2 e -- | Returns the eigenvector with the second smallest eigenvalue of the -- symmetric normalized Laplacian L. Computes real symmetric part of L, so