diff --git a/modularity.cabal b/modularity.cabal
--- a/modularity.cabal
+++ b/modularity.cabal
@@ -1,9 +1,9 @@
 cabal-version: >=1.10
 name: modularity
-version: 0.2.0.4
+version: 0.2.1.0
 license: GPL-3
 license-file: LICENSE
-copyright: 2018 Gregory W. Schwartz
+copyright: 2019 Gregory W. Schwartz
 maintainer: gsch@pennmedicine.upenn.edu
 author: Gregory W. Schwartz
 homepage: http://github.com/GregorySchwartz/modularity#readme
@@ -31,5 +31,5 @@
         eigen ==3.3.4.1,
         hmatrix >=0.19.0.0,
         sparse-linear-algebra >=0.3.1,
-        spectral-clustering >=0.2.1.1,
+        spectral-clustering >=0.2.2.0,
         vector >=0.12.0.1
diff --git a/src/Math/Modularity/Dense.hs b/src/Math/Modularity/Dense.hs
--- a/src/Math/Modularity/Dense.hs
+++ b/src/Math/Modularity/Dense.hs
@@ -7,11 +7,13 @@
 
 module Math.Modularity.Dense
     ( getModularity
+    , getBModularity
     , Q (..)
     ) where
 
 -- Remote
 import Data.Bool (bool)
+import Math.Clustering.Spectral.Dense (B (..), getB)
 import qualified Data.Vector as V
 import qualified Data.Vector.Storable as VS
 import qualified Numeric.LinearAlgebra as H
@@ -41,3 +43,23 @@
     d = H.vector . fmap H.sumElements . H.toRows $ mat
     s = bool (-1) 1 . (== 0) . H.atIndex moduleVec
     k = H.atIndex d
+
+-- | Find modularity from a vector of community labels (0 or 1) corresponding to
+-- rows in the normalized matrix B. See Shu et al., "Efficient Spectral
+-- Neighborhood Blocking for Entity Resolution", 2011.
+-- L = sum_i^n sum_j^n A(i,j) - n = 1^TA1 - n = (B^T1)^T(B^T1) - n.
+getBModularity :: LabelVector -> B -> Q
+getBModularity moduleVec (B b) = Q . sum . fmap inner $ [first, second]
+  where
+    inner v = (a v v / l) - ((a v (ones n) / l) ** 2)
+    first  = H.fromColumns [moduleVec]
+    second = H.fromColumns [H.cmap (bool 1 0 . (== 1)) moduleVec]
+    l    = a (ones n) (ones n)
+    a :: H.Matrix Double -> H.Matrix Double -> Double
+    a oneL oneR = ( flip H.atIndex (0, 0)
+                  $ (H.tr (partA oneL)) H.<> (partA oneR)
+                  )
+                - (H.sumElements oneL)
+    partA one = (H.tr b) <> one
+    n    = H.rows b
+    ones x = (x H.>< 1) [1,1..]