diff --git a/src/Core/Graph.fs b/src/Core/Graph.fs index 52a16f6..cdb23a7 100644 --- a/src/Core/Graph.fs +++ b/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 diff --git a/tests/Tests.FSharp/Algebra/Graph.Tests.fs b/tests/Tests.FSharp/Algebra/Graph.Tests.fs index 8168a23..ca170ab 100644 --- a/tests/Tests.FSharp/Algebra/Graph.Tests.fs +++ b/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 ───────── + +[] +let ``largestEigenvalue returns None for empty graph`` () = + let g : Graph = Graph.empty + Graph.largestEigenvalue 1e-9 100 g |> should equal (None: double option) + +[] +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" + +[] +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" + +[] +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)