fix(#169): replace broken QR with Colt EigenvalueDecomposition

The hand-rolled Householder+implicit-QR routine in GraphSpectrumAnalyzer.computeEigenvalues() threw ArrayIndexOutOfBoundsException at qrTridiagonal:342 on every non-trivial graph (K_3, K_4, P_3, S_4, two disjoint edges), making the entire spectral feature dead code. The bug is a one-off in the inner sweep (offdiag[i + 2] = r) combined with an unpadded length-n offdiag array. Even after padding it produced mathematically wrong eigenvalues on K_3 (returned {2, -1, -3} instead of {2, -1, -1}).

Colt 1.2.0 is already on the classpath (Gvisual/lib/colt-1.2.0.jar) and ships a well-tested EigenvalueDecomposition. Delegate to it for both adjacency and Laplacian spectra. The hand-rolled tridiagonalize / qrTridiagonal helpers and the now-unused MAX_ITERATIONS constant are removed.

Adds 12 JUnit tests in GraphSpectrumAnalyzerTest covering the four issue-#169 reproductions (with closed-form expected spectra), derived metrics, and edge cases. All 12 pass locally.

Closes #169.

Author: sauravbhattacharya001

Gvisual/src/gvisual/GraphSpectrumAnalyzer.java

@@ -2,6 +2,9 @@
 
 import edu.uci.ics.jung.graph.Graph;
 import java.util.*;
+import cern.colt.matrix.DoubleMatrix2D;
+import cern.colt.matrix.impl.DenseDoubleMatrix2D;
+import cern.colt.matrix.linalg.EigenvalueDecomposition;
 
 /**
  * Computes the eigenvalue spectrum of the adjacency and Laplacian matrices
@@ -37,7 +40,6 @@
 public class GraphSpectrumAnalyzer {
 
     private static final double EPSILON = 1e-10;
-    private static final int MAX_ITERATIONS = 500;
 
     private final Graph<String, Edge> graph;
     private double[] adjacencyEigenvalues;
@@ -238,126 +240,35 @@ public String getSummary() {
         return sb.toString();
     }
 
-    // ---- QR Algorithm for symmetric matrix eigenvalues ----
+    // ---- Eigenvalue computation (Colt EigenvalueDecomposition) ----
 
+    /**
+     * Computes the eigenvalues of a real symmetric matrix using Colt's
+     * {@link EigenvalueDecomposition}.
+     *
+     * <p>Previously this class shipped a hand-rolled Householder+implicit-QR
+     * routine that crashed with {@code ArrayIndexOutOfBoundsException} on
+     * every non-trivial graph (issue #169) and produced incorrect eigenvalues
+     * even when it did not crash. Colt is already on the classpath
+     * ({@code lib/colt-1.2.0.jar}), so we delegate to its well-tested
+     * symmetric-matrix path.</p>
+     */
     private double[] computeEigenvalues(double[][] matrix) {
         int n = matrix.length;
         if (n == 0) return new double[0];
         if (n == 1) return new double[]{matrix[0][0]};
 
-        // Copy matrix
-        double[][] M = new double[n][n];
+        DoubleMatrix2D m = new DenseDoubleMatrix2D(n, n);
         for (int i = 0; i < n; i++) {
-            M[i] = Arrays.copyOf(matrix[i], n);
-        }
-
-        // Reduce to tridiagonal form using Householder reflections
-        double[] diag = new double[n];
-        double[] offdiag = new double[n];
-        tridiagonalize(M, diag, offdiag);
-
-        // Apply implicit QR shifts to tridiagonal matrix
-        qrTridiagonal(diag, offdiag, n);
-
-        return diag;
-    }
-
-    /**
-     * Householder tridiagonalization for symmetric matrix.
-     * After this, diag contains diagonal and offdiag contains sub-diagonal.
-     */
-    private void tridiagonalize(double[][] M, double[] diag, double[] offdiag) {
-        int n = M.length;
-        for (int i = n - 1; i > 0; i--) {
-            double scale = 0;
-            for (int k = 0; k < i; k++) scale += Math.abs(M[i][k]);
-
-            if (scale < EPSILON) {
-                offdiag[i] = M[i][i - 1];
-            } else {
-                double h = 0;
-                for (int k = 0; k < i; k++) {
-                    M[i][k] /= scale;
-                    h += M[i][k] * M[i][k];
-                }
-                double f = M[i][i - 1];
-                double g = f >= 0 ? -Math.sqrt(h) : Math.sqrt(h);
-                offdiag[i] = scale * g;
-                h -= f * g;
-                M[i][i - 1] = f - g;
-
-                double[] u = new double[i];
-                for (int j = 0; j < i; j++) {
-                    u[j] = 0;
-                    for (int k = 0; k <= j; k++) u[j] += M[j][k] * M[i][k];
-                    for (int k = j + 1; k < i; k++) u[j] += M[k][j] * M[i][k];
-                    u[j] /= h;
-                }
-
-                double k2 = 0;
-                for (int j = 0; j < i; j++) k2 += M[i][j] * u[j];
-                k2 /= (2.0 * h);
-
-                for (int j = 0; j < i; j++) u[j] -= k2 * M[i][j];
-
-                for (int j = 0; j < i; j++) {
-                    for (int k2b = 0; k2b <= j; k2b++) {
-                        M[j][k2b] -= M[i][j] * u[k2b] + u[j] * M[i][k2b];
-                    }
-                }
-            }
-            diag[i] = M[i][i];
-        }
-        diag[0] = M[0][0];
-        offdiag[0] = 0;
-    }
-
-    /**
-     * Implicit QR algorithm with Wilkinson shifts for tridiagonal symmetric matrix.
-     */
-    private void qrTridiagonal(double[] diag, double[] offdiag, int n) {
-        for (int l = 0; l < n; l++) {
-            int iter = 0;
-            while (true) {
-                int m = l;
-                while (m < n - 1) {
-                    double dd = Math.abs(diag[m]) + Math.abs(diag[m + 1]);
-                    if (Math.abs(offdiag[m + 1]) + dd == dd) break;
-                    m++;
-                }
-                if (m == l) break;
-                if (++iter > MAX_ITERATIONS) break;
-
-                // Wilkinson shift
-                double g = (diag[l + 1] - diag[l]) / (2.0 * offdiag[l + 1]);
-                double r = Math.sqrt(g * g + 1.0);
-                double shift = diag[m] - diag[l]
-                        + offdiag[l + 1] / (g + (g >= 0 ? r : -r));
-
-                double s = 1.0, c = 1.0, p = 0.0;
-                for (int i = m - 1; i >= l; i--) {
-                    double f = s * offdiag[i + 1];
-                    double b = c * offdiag[i + 1];
-                    r = Math.sqrt(f * f + shift * shift);
-                    offdiag[i + 2] = r;
-                    if (Math.abs(r) < EPSILON) {
-                        diag[i + 1] -= p;
-                        offdiag[m + 1 < n ? m + 1 : m] = 0;
-                        break;
-                    }
-                    s = f / r;
-                    c = shift / r;
-                    g = diag[i + 1] - p;
-                    r = (diag[i] - g) * s + 2.0 * c * b;
-                    p = s * r;
-                    diag[i + 1] = g + p;
-                    shift = c * r - b;
-                }
-                diag[l] -= p;
-                offdiag[l + 1] = shift;
-                if (m + 1 < n) offdiag[m + 1] = 0;
+            double[] row = matrix[i];
+            for (int j = 0; j < n; j++) {
+                m.setQuick(i, j, row[j]);
             }
         }
+        EigenvalueDecomposition eig = new EigenvalueDecomposition(m);
+        double[] real = eig.getRealEigenvalues().toArray();
+        // Imaginary parts must be ~0 for symmetric input; we trust Colt here.
+        return real;
     }
 
     private void reverseArray(double[] arr) {

Gvisual/src/test/gvisual/GraphSpectrumAnalyzerTest.java

@@ -0,0 +1,189 @@
+package gvisual;
+
+import edu.uci.ics.jung.graph.Graph;
+import edu.uci.ics.jung.graph.UndirectedSparseGraph;
+import org.junit.Test;
+
+import java.util.Arrays;
+
+import static org.junit.Assert.*;
+
+/**
+ * Tests for {@link GraphSpectrumAnalyzer}.
+ *
+ * <p>Verifies the fix for issue #169: the previous hand-rolled QR routine
+ * threw {@code ArrayIndexOutOfBoundsException} on every non-trivial graph
+ * and returned mathematically wrong eigenvalues even when it did not crash.
+ * After delegating to Colt's {@code EigenvalueDecomposition}, the analyzer
+ * matches closed-form spectra for canonical graph families.</p>
+ */
+public class GraphSpectrumAnalyzerTest {
+
+    private static final double TOL = 1e-6;
+
+    private static Edge addEdge(Graph<String, Edge> g, String u, String v) {
+        if (!g.containsVertex(u)) g.addVertex(u);
+        if (!g.containsVertex(v)) g.addVertex(v);
+        Edge e = new Edge("f", u, v);
+        e.setWeight(1f);
+        g.addEdge(e, u, v);
+        return e;
+    }
+
+    private static Graph<String, Edge> emptyGraph() {
+        return new UndirectedSparseGraph<String, Edge>();
+    }
+
+    private static Graph<String, Edge> completeGraph(int n) {
+        Graph<String, Edge> g = emptyGraph();
+        for (int i = 0; i < n; i++) g.addVertex("v" + i);
+        for (int i = 0; i < n; i++) {
+            for (int j = i + 1; j < n; j++) {
+                addEdge(g, "v" + i, "v" + j);
+            }
+        }
+        return g;
+    }
+
+    /** Path graph on n vertices: v0 - v1 - ... - v(n-1). */
+    private static Graph<String, Edge> pathGraph(int n) {
+        Graph<String, Edge> g = emptyGraph();
+        for (int i = 0; i < n; i++) g.addVertex("v" + i);
+        for (int i = 0; i + 1 < n; i++) addEdge(g, "v" + i, "v" + (i + 1));
+        return g;
+    }
+
+    /** Star graph K_{1, n-1}: hub v0 connected to v1..v(n-1). */
+    private static Graph<String, Edge> starGraph(int n) {
+        Graph<String, Edge> g = emptyGraph();
+        g.addVertex("hub");
+        for (int i = 1; i < n; i++) addEdge(g, "hub", "leaf" + i);
+        return g;
+    }
+
+    // ─── Issue #169 reproductions: these used to throw ──────────────
+
+    /** Triangle K_3. Old code: ArrayIndexOutOfBoundsException at line 342. */
+    @Test
+    public void triangleSpectrumDoesNotCrash() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(completeGraph(3));
+        a.compute();
+        double[] eig = a.getAdjacencyEigenvalues();
+        assertEquals(3, eig.length);
+        // K_3 adjacency spectrum: {2, -1, -1} (descending).
+        assertEquals(2.0, eig[0], TOL);
+        assertEquals(-1.0, eig[1], TOL);
+        assertEquals(-1.0, eig[2], TOL);
+    }
+
+    @Test
+    public void completeK4Spectrum() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(completeGraph(4));
+        a.compute();
+        double[] eig = a.getAdjacencyEigenvalues();
+        // K_n adjacency spectrum: {n-1, -1, -1, ..., -1}.
+        assertEquals(4, eig.length);
+        assertEquals(3.0, eig[0], TOL);
+        for (int i = 1; i < 4; i++) assertEquals(-1.0, eig[i], TOL);
+    }
+
+    @Test
+    public void path3Spectrum() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(pathGraph(3));
+        a.compute();
+        double[] eig = a.getAdjacencyEigenvalues();
+        // P_3 spectrum: { sqrt(2), 0, -sqrt(2) }.
+        assertEquals(3, eig.length);
+        assertEquals(Math.sqrt(2.0), eig[0], TOL);
+        assertEquals(0.0, eig[1], TOL);
+        assertEquals(-Math.sqrt(2.0), eig[2], TOL);
+    }
+
+    @Test
+    public void star4Spectrum() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(starGraph(4));
+        a.compute();
+        double[] eig = a.getAdjacencyEigenvalues();
+        // K_{1,3} spectrum: { sqrt(3), 0, 0, -sqrt(3) }.
+        assertEquals(4, eig.length);
+        assertEquals(Math.sqrt(3.0), eig[0], TOL);
+        assertEquals(0.0, eig[1], TOL);
+        assertEquals(0.0, eig[2], TOL);
+        assertEquals(-Math.sqrt(3.0), eig[3], TOL);
+    }
+
+    // ─── Derived metrics ────────────────────────────────────────────
+
+    @Test
+    public void spectralRadiusMatchesLargestEigenvalue() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(completeGraph(5));
+        a.compute();
+        // K_5 spectral radius = 4.
+        assertEquals(4.0, a.getSpectralRadius(), TOL);
+    }
+
+    @Test
+    public void algebraicConnectivityOfDisconnectedGraphIsZero() {
+        // Two disjoint edges: {a-b, c-d}.
+        Graph<String, Edge> g = emptyGraph();
+        addEdge(g, "a", "b");
+        addEdge(g, "c", "d");
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(g);
+        a.compute();
+        // Fiedler value (second-smallest Laplacian eigenvalue) is 0 for a
+        // disconnected graph.
+        assertEquals(0.0, a.getAlgebraicConnectivity(), TOL);
+        assertEquals(2, a.getComponentCount());
+    }
+
+    @Test
+    public void algebraicConnectivityOfTriangleIsThree() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(completeGraph(3));
+        a.compute();
+        // L(K_3) spectrum = {0, 3, 3}. Fiedler value = 3.
+        assertEquals(3.0, a.getAlgebraicConnectivity(), TOL);
+        assertEquals(1, a.getComponentCount());
+    }
+
+    @Test
+    public void graphEnergyOfTriangle() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(completeGraph(3));
+        a.compute();
+        // Energy(K_3) = |2| + |-1| + |-1| = 4.
+        assertEquals(4.0, a.getGraphEnergy(), TOL);
+    }
+
+    // ─── Edge cases ─────────────────────────────────────────────────
+
+    @Test
+    public void emptyGraphIsHandled() {
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(emptyGraph());
+        a.compute();
+        assertEquals(0, a.getAdjacencyEigenvalues().length);
+        assertEquals(0, a.getLaplacianEigenvalues().length);
+        assertEquals(0.0, a.getSpectralRadius(), TOL);
+        assertEquals(0, a.getComponentCount());
+    }
+
+    @Test
+    public void singleVertexGraphIsHandled() {
+        Graph<String, Edge> g = emptyGraph();
+        g.addVertex("only");
+        GraphSpectrumAnalyzer a = new GraphSpectrumAnalyzer(g);
+        a.compute();
+        // Single vertex: adjacency = [0], Laplacian = [0].
+        assertArrayEquals(new double[]{0.0}, a.getAdjacencyEigenvalues(), TOL);
+        assertArrayEquals(new double[]{0.0}, a.getLaplacianEigenvalues(), TOL);
+        assertEquals(1, a.getComponentCount());
+    }
+
+    @Test(expected = IllegalArgumentException.class)
+    public void nullGraphRejected() {
+        new GraphSpectrumAnalyzer(null);
+    }
+
+    @Test(expected = IllegalStateException.class)
+    public void queryBeforeComputeRejected() {
+        new GraphSpectrumAnalyzer(completeGraph(3)).getSpectralRadius();
+    }
+}