diff --git a/cpp/graph_spectrum.hpp b/cpp/graph_spectrum.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4be14e52ec06aff130a54c02533c228d011857a6
--- /dev/null
+++ b/cpp/graph_spectrum.hpp
@@ -0,0 +1,65 @@
+#include "graph.hpp"
+#include <Eigen/Dense>
+#include <Eigen/Eigenvalues>
+
+using MatrixType =
+    Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
+
+// A: Adjacency matrix
+//    lambda_max <= d_max
+//
+// L: Laplacian
+//    L = D - A
+///
+// P: Random walk matrix
+//    lambda_max = 1
+
+class GraphSpectrum {
+  public:
+    GraphSpectrum(const Graph& g) : graph(g) {}
+    ~GraphSpectrum() {}
+
+    std::vector<float> computeAdjacencySpectrum() const {
+        // matrix stored as std::vector<std::vector<bool>>
+        auto& badj = graph.getBooleanAdj();
+
+        // Convert it to MatrixType
+        auto n = badj.size();
+        MatrixType m(n, n);
+        for (auto i = 0u; i < n; ++i)
+            for (auto j = 0u; j < n; ++j)
+                m(i, j) = badj[i][j] ? 1.0f : 0.0f;
+
+        return getEigenvalues_(m);
+    }
+
+    std::vector<float> computeLaplacianSpectrum() const {
+        // matrix stored as std::vector<std::vector<bool>>
+        auto& badj = graph.getBooleanAdj();
+        auto& adj = graph.getAdj();
+
+        // - A
+        auto n = badj.size();
+        MatrixType m(n, n);
+        for (auto i = 0u; i < n; ++i)
+            for (auto j = 0u; j < n; ++j)
+                m(i, j) = badj[i][j] ? -1.0f : 0.0f;
+
+        // + D
+        for (auto i = 0u; i < n; ++i)
+            m(i, i) = float(adj[i].size());
+
+        return getEigenvalues_(m);
+    }
+
+  private:
+    const Graph& graph;
+
+    std::vector<float> getEigenvalues_(const MatrixType& m) const {
+        Eigen::SelfAdjointEigenSolver<MatrixType> es(
+            m, Eigen::DecompositionOptions::EigenvaluesOnly);
+        auto ev = es.eigenvalues();
+        return std::vector<float>(ev.data(), ev.data() + ev.rows() * ev.cols());
+    }
+};
+