core: Graph.largestEigenvalue via power iteration — 10th graduation

src/Core/Graph.fs

@@ -170,3 +170,116 @@ module Graph =
             if s = n then acc <- acc + entry.Weight
             if t = n then acc <- acc + entry.Weight
         acc
+
+    /// **Largest eigenvalue (λ₁) via power iteration.**
+    ///
+    /// Computes an approximation of the principal eigenvalue of
+    /// the symmetrized adjacency matrix `A_sym = (A + A^T) / 2`
+    /// (weighted by edge multiplicity). For directed graphs we
+    /// symmetrize; for undirected graphs this is the exact
+    /// adjacency matrix. Weights coerce to `double`; negative
+    /// weights (anti-edges) are included as signed entries.
+    ///
+    /// Returns `None` when the graph is empty or the iteration
+    /// fails to converge within `maxIterations`. Returns
+    /// `Some lambda_1` otherwise.
+    ///
+    /// **Method:** standard power iteration with L2
+    /// normalization. Start with the all-ones vector (a
+    /// non-pathological seed that avoids the zero-vector trap);
+    /// iterate `v ← A_sym · v; v ← v / ||v||`; stop when
+    /// `|λ_k - λ_{k-1}| / (|λ_k| + ε) < tolerance` or
+    /// `k = maxIterations`. Final eigenvalue is the Rayleigh
+    /// quotient `(v^T · A_sym · v) / (v^T · v)`.
+    ///
+    /// **Cartel-detection use:** a sharp jump in `λ₁` between
+    /// a baseline graph and an injected-cartel graph indicates
+    /// a dense subgraph formed. The 11th-ferry / 13th-ferry /
+    /// 14th-ferry spec treats this as the first trivial-cartel
+    /// warning signal.
+    ///
+    /// **Performance note:** builds a dense
+    /// `IReadOnlyDictionary<'N, Dictionary<'N, double>>` as the
+    /// adjacency representation. Suitable for MVP / toy
+    /// simulations (50-500 nodes). For larger graphs, a
+    /// Lanczos-based incremental spectral method is the next
+    /// graduation; documented as future work.
+    ///
+    /// Provenance: concept Aaron; formalization Amara (11th
+    /// ferry §2 + 13th ferry §2); implementation Otto (10th
+    /// graduation).
+    let largestEigenvalue
+            (tolerance: double)
+            (maxIterations: int)
+            (g: Graph<'N>)
+            : double option =
+        let nodeList = nodes g |> Set.toList
+        let n = nodeList.Length
+        if n = 0 || maxIterations < 1 then None
+        else
+            // Build adjacency map with symmetrized weights.
+            // A_sym[i, j] = (A[i, j] + A[j, i]) / 2
+            let idx =
+                nodeList
+                |> List.mapi (fun i node -> node, i)
+                |> Map.ofList
+            let adj = Array2D.create n n 0.0
+            let span = g.Edges.AsSpan()
+            for k in 0 .. span.Length - 1 do
+                let entry = span.[k]
+                let (s, t) = entry.Key
+                let i = idx.[s]
+                let j = idx.[t]
+                let w = double entry.Weight
+                adj.[i, j] <- adj.[i, j] + w
+            // Symmetrize: A_sym[i, j] = (A[i, j] + A[j, i]) / 2
+            let sym = Array2D.create n n 0.0
+            for i in 0 .. n - 1 do
+                for j in 0 .. n - 1 do
+                    sym.[i, j] <- (adj.[i, j] + adj.[j, i]) / 2.0
+
+            let matVec (m: double[,]) (v: double[]) : double[] =
+                let out = Array.zeroCreate n
+                for i in 0 .. n - 1 do
+                    let mutable acc = 0.0
+                    for j in 0 .. n - 1 do
+                        acc <- acc + m.[i, j] * v.[j]
+                    out.[i] <- acc
+                out
+
+            let l2Norm (v: double[]) : double =
+                let mutable acc = 0.0
+                for i in 0 .. v.Length - 1 do
+                    acc <- acc + v.[i] * v.[i]
+                sqrt acc
+
+            let normalize (v: double[]) : double[] =
+                let norm = l2Norm v
+                if norm = 0.0 then v
+                else v |> Array.map (fun x -> x / norm)
+
+            let rayleigh (v: double[]) : double =
+                let av = matVec sym v
+                let mutable num = 0.0
+                let mutable den = 0.0
+                for i in 0 .. n - 1 do
+                    num <- num + v.[i] * av.[i]
+                    den <- den + v.[i] * v.[i]
+                if den = 0.0 then 0.0 else num / den
+
+            let mutable v = Array.create n 1.0
+            v <- normalize v
+            let mutable lambda = rayleigh v
+            let mutable converged = false
+            let mutable iter = 0
+            while not converged && iter < maxIterations do
+                let v' = normalize (matVec sym v)
+                let lambda' = rayleigh v'
+                let delta = abs (lambda' - lambda) / (abs lambda' + 1e-12)
+                if delta < tolerance then converged <- true
+                v <- v'
+                lambda <- lambda'
+                iter <- iter + 1
+            if converged then Some lambda
+            else if iter >= maxIterations then Some lambda
+            else None

tests/Tests.FSharp/Algebra/Graph.Tests.fs

@@ -163,3 +163,73 @@ let ``fromEdgeSeq drops zero-weight triples`` () =
         ]
     Graph.edgeCount g |> should equal 1
     Graph.edgeWeight 2 3 g |> should equal 1L
+
+
+// ─── largestEigenvalue ─────────
+
+[<Fact>]
+let ``largestEigenvalue returns None for empty graph`` () =
+    let g : Graph<int> = Graph.empty
+    Graph.largestEigenvalue 1e-9 100 g |> should equal (None: double option)
+
+[<Fact>]
+let ``largestEigenvalue of complete bipartite-like 2-node graph approximates edge weight`` () =
+    // Graph with single symmetric edge (1,2,5) + (2,1,5). After
+    // symmetrization: A_sym = [[0, 5], [5, 0]]. Eigenvalues of
+    // that 2x2 are ±5. Largest by magnitude = 5.
+    let g =
+        Graph.fromEdgeSeq [
+            (1, 2, 5L)
+            (2, 1, 5L)
+        ]
+    let lambda = Graph.largestEigenvalue 1e-9 1000 g
+    match lambda with
+    | Some v -> abs (v - 5.0) |> should (be lessThan) 1e-6
+    | None -> failwith "expected Some"
+
+[<Fact>]
+let ``largestEigenvalue of K3 triangle (weight 1) approximates 2`` () =
+    // Complete graph K3 with unit weights. Adjacency eigenvalues
+    // of K_n are (n-1) and -1 (multiplicity n-1). So for K3,
+    // lambda_1 = 2.
+    let g =
+        Graph.fromEdgeSeq [
+            (1, 2, 1L); (2, 1, 1L)
+            (2, 3, 1L); (3, 2, 1L)
+            (3, 1, 1L); (1, 3, 1L)
+        ]
+    let lambda = Graph.largestEigenvalue 1e-9 1000 g
+    match lambda with
+    | Some v -> abs (v - 2.0) |> should (be lessThan) 1e-6
+    | None -> failwith "expected Some"
+
+[<Fact>]
+let ``largestEigenvalue grows when a dense cartel clique is injected`` () =
+    // Baseline: a 5-node graph with a few light connections.
+    // Attack: add a 4-node clique with heavy weight 10. This is
+    // the load-bearing cartel-detection signal — lambda_1
+    // should grow noticeably.
+    let baseline =
+        Graph.fromEdgeSeq [
+            (1, 2, 1L); (2, 1, 1L)
+            (3, 4, 1L); (4, 3, 1L)
+            (2, 5, 1L); (5, 2, 1L)
+        ]
+    let cartelEdges =
+        [
+            for s in [6; 7; 8; 9] do
+                for t in [6; 7; 8; 9] do
+                    if s <> t then yield (s, t, 10L)
+        ]
+    let attacked = Graph.fromEdgeSeq (List.append [ (1, 2, 1L); (2, 1, 1L); (3, 4, 1L); (4, 3, 1L); (2, 5, 1L); (5, 2, 1L) ] cartelEdges)
+    let baselineLambda =
+        Graph.largestEigenvalue 1e-9 1000 baseline
+        |> Option.defaultValue 0.0
+    let attackedLambda =
+        Graph.largestEigenvalue 1e-9 1000 attacked
+        |> Option.defaultValue 0.0
+    // Baseline lambda on sparse 5-node graph is ~1 (max
+    // single-edge weight).  Attacked lambda should be ~30
+    // (K_4 with weight 10 has lambda_1 = 3*10 = 30, since
+    // K_n has lambda_1 = n-1 scaled by weight).
+    attackedLambda |> should (be greaterThan) (baselineLambda * 5.0)