Differentiable Spectral Normalization for Large-Scale Ising Optimization

1 Introduction

Ising optimization arises in diverse applications including correlation clustering [Bansal et al., 2004], MAP inference in graphical models [Wainwright and Jordan, 2008], and energy-based learning [Hopfield, 1982; Ackley et al., 1985]. Modern instances involve sparse, signed graphs with 105 -107 variables [Lucas, 2014; Kunegis et al., 2010], where computational efficiency and solution quality are essential.

Existing approaches exhibit a tradeoff between these objectives. Spectral relaxation reduces the discrete problem to a sparse eigenvalue computation, achieving near-linear time complexity [Delorme and Poljak, 1993; Trevisan, 2012]. However, the quality of spectral solutions depends on the choice of graph normalization, typically the adjacency matrix or degree-normalized Laplacian, which is selected a priori without reference to the problem instance. Empirical studies indicate that no normalization performs well across diverse graph families [Cucuringu et al., 2019; Cucuringu et al., 2021]. Combinatorial metaheuristics such as simulated annealing [Kirkpatrick et al., 1983] and parallel tempering [Swendsen and Wang, 1986] achieve better solution quality by directly exploring the discrete energy landscape, but their reliance on sequential sampling limits scalability: runtimes grow prohibitively for graphs beyond 105 nodes.

∗Corresponding author

This paper investigates whether the normalization that defines a spectral relaxation can be optimized. We introduce Differentiable Spectral Normalization (DSN), which treats normalization as a continuous optimization variable. The method parameterizes the relaxation through a positive diagonal matrix W that defines a weighted inner product, and maximizes the spectral lower bound via projected gradient ascent. The key computational insight is that gradients of eigenvalues with respect to W admit a closed-form expression via the Hellmann-Feynman theorem [Zhang and George, 2002; Milgrom and Segal, 2002], requiring only the principal eigenpair and preserving O(|E|) complexity per iteration.

The main contributions are as follows:

• We formulate normalization selection as bound maximization over a family of diagonal metrics, providing a principled alternative to heuristic operator design.
• We derive an efficient gradient expression that avoids backpropagation through the eigensolver, enabling optimization on graphs with millions of nodes.
• We demonstrate empirically that the learned normalization yields tighter bounds that correlate strongly with improved discrete solutions (Figure 1).

On standard benchmarks (Gset, signed social networks) and synthetic graphs up to 8.4 × 106 nodes, DSN improves solution quality by 3-15% over static spectral baselines and achieves energies within 1-3% of simulated annealing at up to 190× lower runtime. These results indicate that optimizing relaxation geometry can narrow the gap between spectral efficiency and metaheuristic quality, without sacrificing the scalability advantages of eigenvalue-based methods.

The remainder of this paper is organized as follows. Section 2 defines the Ising problem and spectral relaxation framework. Section 3 derives the DSN update rule. Section 4

analyzes smoothness and convergence properties. Section 5 presents experimental results.

Related works. Approaches for large-scale Ising optimization can be broadly categorized into three paradigms: spectral relaxations, semidefinite programming, and combinatorial metaheuristics.

Spectral relaxation serves as a cornerstone for scalable graph clustering and Max-Cut due to its efficiency in leveraging sparse linear algebra [Delorme and Poljak, 1993; Trevisan, 2012]. These methods relax discrete constraints to a spherical manifold and solve for the principal eigenvectors of a graph operator. To handle signed interactions, foundational works have proposed specific operator constructions, such as the signed Laplacian and generalized eigenproblems [Kunegis et al., 2010; Cucuringu et al., 2019]. More recently, regularization techniques have been introduced to stabilize embeddings for sparse, heterogeneous graphs [Cucuringu et al., 2021]. However, a critical limitation remains: these solvers rely on fixed, axiomatic normalizations (e.g., degreebased scaling) that are selected heuristically rather than optimized. Such static choices often fail to capture the complex geometry of ground-truth solutions in massive networks. DSN generalizes this paradigm by treating the normalization matrix not as a hyperparameter, but as a continuous variable learned directly from the Ising objective.

Pioneered by Goemans and Williamson [Goemans and Williamson, 1995], Semidefinite programming (SDP) relaxations provide rigorous approximation guarantees by optimizing the embedding geometry over the cone of positive semidefinite matrices. Unlike spectral methods that project data onto a low-dimensional subspace, SDP allows for a high-dimensional representation that tightens the relaxation gap. However, standard interior-point methods for SDP scale poorly (typically O(N 3.5 ) or worse), rendering them intractable for graphs with N > 104 . Our proposed DSN strikes a balance: it can be viewed as a scalable "diagonal SDP" that restricts the metric structure to a diagonal form. This restriction allows DSN to retain the near-linear O(|E|) complexity of spectral methods while capturing the essential geometric flexibility required for high-quality solutions.

Metaheuristics such as Simulated Annealing (SA) [Kirkpatrick et al., 1983], Parallel Tempering [Swendsen and Wang, 1986], and Tabu Search directly explore the discrete energy landscape. These methods are capable of escaping local minima and finding near-optimal states on moderatesized instances. However, their performance relies heavily on extensive sequential sampling, long burn-in phases, and delicate cooling schedules [Lucas, 2014]. For massive graphs (N > 106 ), the convergence time becomes prohibitive. DSN complements these approaches by providing a fast, deterministic method to generate high-quality embeddings, which can serve as superior warm-start configurations for fine-tuning heuristics.

Optimizing objectives that depend on spectral quantities requires computing gradients of eigenvalues and eigenvectors. While sensitivity analysis for matrix decompositions is well-established [Magnus, 1985; Kato, 1995], applying it to large-scale optimization is challenging. Recent machine-learning approaches, spectral normalization [Miyato et al., 2018] and differentiable eigendecomposition [Ionescu et al., 2015; Wang et al., 2019], instead rely on automatic differentiation through iterative eigensolvers, which is memory-intensive and prone to numerical instability on large graphs. Our work circumvents these issues by leveraging the Hellmann-Feynman theorem [Zhang and George, 2002; Milgrom and Segal, 2002]. This allows us to derive exact analytical gradients using only the principal eigenpair, enabling efficient geometry updates without the need for expensive backpropagation. Complementary work applies learning to combinatorial optimization, GNN solvers [Schuetz et al., 2022], RL [Dai et al., 2017], unsupervised objectives [Karalias and Loukas, 2020], but typically requires training data; DSN is orthogonal, optimising relaxation geometry perinstance without any.

2 Preliminaries

2.1 Ising Optimization Problem

We consider the problem of finding a binary spin configuration s ∈ {-1, +1}n that minimizes the Ising energy (Hamiltonian):

In this convention, a positive interaction Jij > 0 (ferromagnetic) encourages alignment (si = sj) to lower the energy contribution to -Jij, while Jij < 0 (antiferromagnetic) encourages disagreement. Finding the exact global minimizer E ⋆ is generally NP-hard.

2.2 Signed Laplacian and Spectral Duality

To extend spectral graph theory to signed interactions, we utilize the Signed Laplacian ¯L. Let D be the diagonal signed degree matrix with entries Dii ≜ ∑j |Jij |. The Signed Laplacian is defined as ¯L ≜ D - J. Leveraging the boolean constraint s2i = 1, we establish the duality between the Laplacian quadratic form and the Ising energy:

where C = ∑i,j |Jij| is a constant total weight. Consequently, minimizing the Ising energy E(s) is mathematically equivalent to minimizing the quadratic form of the Signed Laplacian.

2.3 Spectral Relaxations and Lower Bounds

The combinatorial hardness of Eq. (1) stems from the discrete hypercube constraint s ∈ {-1,+1}n . Spectral methods address this by relaxing the domain to a continuous manifold.

A standard approach relaxes the discrete condition to a spherical constraint S n-1 = {u ∈ Rn : u⊤u = n}. Substituting the continuous variable u into the objective yields the eigenvalue problem: minu∈Rn -u⊤Ju s.t. u⊤u = n.

Using the algebraic identity min(-x) = -max(x), the optimal solution is given by the eigenvector corresponding to the largest eigenvalue λmax(J). This provides a loose spectral lower bound: E ⋆ ≥ -n · λmax(J). More generally, the relaxation geometry can be parameterized by a positive diagonal matrix W = diag(w) ≻ 0. Since any discrete spin s satisfies s⊤Ws = ∑ wis2i = Tr(W), we can relax the domain to a weighted hypersphere:

This generalization leads to a fundamental lower bound that depends on the choice of the metric W.

Lemma 1. For any positive diagonal matrix W, the optimal Ising energy is lower-bounded by the principal eigenvalue of the normalized coupling matrix:

Existing spectral methods typically commit to a fixed, heuristic geometry (W = I or W = D), often resulting in loose bounds. In contrast, our methodology (Section 3) treats B(W) as an objective function to be maximized.

3 Methodology

3.1 Learning The Relaxation Geometry via A Diagonal Metric

Classical spectral solvers fix the relaxation geometry by committing to a static operator (the unweighted matrix J or the normalized Laplacian). We instead relax the discrete hypercube onto a learnable weighted hypersphere parameterized by a positive diagonal matrix W = diag(w1, . . . , wn) ≻ 0. Consider the continuous relaxation of the Ising energy minimization problem defined in Eq. (1):

Geometrically, W defines a diagonal metric that deforms the constraint surface, effectively stretching or compressing the embedding space to align the spectral solution with the discrete problem structure. Solving this constrained optimization via Lagrange multipliers involves finding the stationary points of the Rayleigh quotient. Using the identity min(-x) = -max(x), the minimum relaxed energy is determined by the principal eigenpair (λmax, u) of the normalized operator: N(W) ≜W -1/2JW -1/2 , with ∥u∥2 = 1.

This relaxation induces a rigorous spectral lower bound on the optimal ground-state energy: E ⋆ ≥ B(W) ≜ -λmax(N(W)) · Tr(W).

Our goal is to learn the optimal normalization by maximizing this lower bound (tightening the relaxation) over a compact domain:

The domain Dε prevents degenerate scalings and ensures the numerical stability of the eigensolver.

3.2 Gradient Derivation via Eigenvalue Perturbation

A central technical challenge is that the objective B(W) depends on λmax(N(W)). Differentiating through iterative eigensolvers (Lanczos) via unrolling is memory-intensive and often numerically unstable. Instead, we leverage first-order eigenvalue perturbation theory (Hellmann-Feynman theorem). Let (λ, u) be the dominant eigenpair of N(W) with ∥u∥2 = 1. The gradient of the eigenvalue with respect to a diagonal parameter wi is: ∂λ ∂wi = u⊤∂N(W) ∂wi u. Using the identity ∂W -1/2 ∂wi = - 12 w-3/2 i Eii, where Eii ∈ Rn×n is the elementary matrix with a 1 at entry (i, i) and 0 elsewhere, the derivative of the operator is:

Substituting this into the bilinear form and utilizing N(W)u = λu, we obtain:

Applying the product rule to the full objective B(W) = -λTr(W) (where Tr(W) = ∑ wk), we derive the closedform gradient: ∂B(W) ∂wi = - ( ∂Tr(W) ∂wi λ+Tr(W) ∂λ ∂wi ) = - (1 · λ + Tr(W) (- λu2i wi )) = λ ( Tr(W ) wi u2i -1) .

This gradient expression relies solely on the leading eigenpair (λ, u), enabling efficient computation without backpropagation. Intuitively, ∇wi B acts as a localization penalty on node i: if a node i exhibits localization (large u2i ), the gradient is positive (assuming λ > 0), prompting an increase in wi. This penalizes the node in the weighted norm constraint, effectively "flattening" the eigenvector to enforce a more distributed embedding.

Algorithm 1 Differentiable Spectral Normalization (DSN)

Input: Coupling matrix J; step size η; iterations T .

Output: Discrete spin s ∈ {±1}n .

• 1: Initialize: sbest ← ∅, Ebest ←∞.
• 2: Initialize vinit randomly, W (0) ←D.
• 3: for t = 0 to T do
• 4: 1. Spectral solving: Compute (λt, ut) of:

• 5: 2. Discrete sampling: s(t) ← sign((W (t) )-1/2ut)

• 6: 3. Update best solution:
• 7: if E(s(t) ) < Ebest then
• 8: sbest ← s(t) , Ebest ← E(s(t) )
• 9: end if
• 10: 4. Gradient evaluation: Compute ∇B:

11: 5. Geometry update:

• 12: Project W (t+1) to Dε to ensure positive definiteness.
• 13: 6. Warm-start: vinit ← ut
• 14: end for
• 15: return sbest

3.3 Algorithm

We optimize B(W) via projected gradient ascent on the manifold Dε. The update rule at iteration t is:

where η > 0 is the step size and ΠDε denotes element-wise projection onto [ε, ε-1 ]. The complete procedure is summarized in Algorithm 1. The computational cost of each iteration is dominated by the sparse matrix-vector products required to compute the leading eigenpair of N(W (t) ). By using efficient iterative solvers (LOBPCG or Lanczos), the complexity remains near-linear in the number of edges |E|, preserving the scalability of classical spectral methods.

Although DSN optimizes a continuous spectral bound, the final objective is a discrete spin configuration. We employ standard sign rounding to discretize the optimized embedding. Regarding convergence, DSN performs projected gradient ascent over a compact domain. Under mild assumptions (Lipschitz continuity of the eigenvalue function, which holds given a non-zero spectral gap), the procedure is guaranteed to converge to a stationary point of the spectral bound B(W).

4 Theoretical Foundations

Definition 2. Let Dn++ = {W = diag(w) ∈ Rn×n | wi > 0, ∀i} be the cone of positive definite diagonal matrices. For any W ∈ Dn++, the weighted spectral operator N(W) ∈ Rn×n is defined as the similarity transformation of the coupling matrix J induced by the metric W:

The eigenvalues of N(W), denoted by λk(N(W)), correspond to the generalized eigenvalues of the pencil (J,W), satisfying the characteristic equation det(J - λW) = 0.

This operator generalizes standard spectral graph matrices. The choice W = I yields the coupling matrix J itself, while W = D (where Dii = ∑j |Jij |) yields the normalized signed operator associated with the random walk Laplacian. The mapping W ↦→ N(W) is smooth on Dn++, ensuring that the spectral properties evolve continuously with respect to the weights w.

Theorem 3. Let D ⊆ Dn++ be the domain where the principal eigenvalue λmax(N(W)) is simple. For any W ∈ D, the spectral lower bound function B(W) is differentiable. Let (λ, u) denote the principal eigenpair of N(W) with ∥u∥2 = 1. The partial derivative of B(W) with respect to the weight wi is given by:

Proof. The objective function is defined as B(W) = -λ(W)Tr(W), where we denote λ(W) = λmax(N(W)) for brevity. By the product rule of differentiation:

The derivative of the trace term is trivial: ∂Tr(W) ∂wi = ∂ ∂wi (∑j wj) = 1. To evaluate the eigenvalue derivative, we invoke the Hellmann-Feynman theorem extended to matrix operators. Since λ is simple, it varies smoothly with W, and its first derivative is given by the expectation of the operator derivative with respect to the eigenvector:

Recall that N(W) = W -1/2JW -1/2 . Let Eii denote the elementary matrix with 1 at entry (i, i) and 0 elsewhere. Differentiating the diagonal scaling matrix yields ∂ ∂wi W -1/2 =

- 12 w-3/2 i Eii = - 1 2wi EiiW -1/2 . Applying the product rule to the operator N(W):

Substituting this into the eigenvalue derivative expression:

Using the eigenvalue equation N(W)u = λu and the symmetry u⊤N(W) = λu⊤:

Since u⊤Eiiu = u2i , we obtain the scalar derivative ∂λ(W) ∂wi = - λu2i

wi . Finally, substituting terms back into the objective gradient equation:

Standard matrix perturbation theory establishes that the principal eigenvalue λmax(N(W)) is a real-analytic function of the matrix entries provided the eigenvalue is simple [Kato, 1995]. Since our feasible set Dε is compact, the existence of a uniform spectral gap γ > 0 (discussed below) ensures that both the objective function B(W) and its gradient are Lipschitz continuous, satisfying the regularity conditions required for optimization.

Lemma 4. Let Dε := {W = diag(w) ∈ Rn×n : ε ≤ wi ≤ ε-1 } be the feasible set. Assume there exists a uniform lower bound γ > 0 such that trajectory {W (t) } remains within the region of non-degenerate spectral gaps:

Then B(W) is L-smooth on Dε,γ. For a step size 0 < η ≤ 1/L, the projected gradient updates satisfy:

1. The objective strictly increases unless a stationary point is reached. With the standard mapping is defined as Gη(W (t) ) = 1 η (W (t) -W (t+1) ), we have:

1. {W (t) } admits accumulation points, and every W ∗ is a first-order stationary point satisfying: ⟨∇B(W ∗ ),W - W ∗ ⟩ ≤ 0, ∀W ∈ Dε.
2. Let B∗ := maxW∈Dε B(W). The minimum squared norm of gradient mapping decays at a rate of O(1/T ):

Proof. (Sketch) Dε,γ is a compact subset. Within this region, the principal eigenvalue is simple, implying that B(W) is C 2 differentiable [Kato, 1995]. Since the Hessian is continuous on a compact set, it is bounded, implying that ∇B is L-Lipschitz continuous along the trajectory.

1. By the descent lemma for L-smooth functions and the variational property of the projection ΠDε , we derive the sufficient ascent inequality: B(W (t+1) ) - B(W (t) ) ≥ 1 2η ∥W (t+1)-W (t) ∥2F . Noting that ∥Gη∥2F = 1 η2 ∥W (t)- W (t+1) ∥2F , the result follows.
2. Since B is continuous and bounded above on Dε and the sequence is strictly increasing (non-stationary), the values converge. Since {W (t) } ⊂ Dε, accumulation points exist. As ∥Gη(W (t) )∥F → 0, the fixed-point condition W ∗ = ΠDε (W ∗ +η∇B(W ∗ )) holds for any limit point W ∗ , verifying stationarity.
3. Summing the sufficient ascent inequality over T iterations yields B∗ - B(W (0) ) ≥ ηT 2 mink ∥Gη(W (k) )∥2F . Rearranging terms provides the O(1/T) rate for the squared norm.

The assumption of a non-zero spectral gap holds generically. Results in random matrix theory show that for random graphs and generic perturbations, the Laplacian spectrum is simple with high probability [Luh et al., 2025; Kato, 1995]. Eigenvalue degeneracy requires exact structural symmetry, which is rare in real-world instances. Nonetheless, for robustness, our implementation monitors the spectral gap γt = λ1(N(W)) - λ2(N(W)) at each iteration. If γt falls below a safeguard threshold 10-2 ∥J∥2, we reduce the stepsize η proportionally. Across all reported experiments, this safeguard was never triggered, confirming the validity of the smoothness assumption in practice; quantitative profiling is reported in §5.

SDP solves a convex relaxation and can produce very tight bounds; its computational cost typically scales cubically or worse, which limits scalability on large instances. In contrast, DSN optimizes a non-convex spectral objective with periteration cost dominated by sparse matrix-vector products, yielding near-linear scaling in |E|. As shown in Figure 2, DSN reached its converged bound within seconds (t ≈ 5s), whereas the SDP took longer (t ≈ 200s) to crossover. While SDP may eventually surpass DSN given sufficient runtime, the final bound gap remains small in our experiments, indicating that DSN provides a highly efficient approximation in the large-scale regime.

5 Experiments

5.1 Experimental Setup

Datasets We evaluate DSN on three instance families:

• Gset benchmark: 66 widely-used Max-Cut instances (N ∈ [800, 40000]) exhibiting diverse topologies, including random and toroidal graphs [Ye, 2003].

(b) Erd˝os-R´enyi dataset.

Figure 3: Scalability and efficiency on large-scale benchmarks (Barab´asi-Albert and Erd˝os-R´enyi). DSN scales up to 223 , uniquely bridging the gap between the speed of spectral heuristics and the quality of an exact solver. In contrast, combinatorial (Gurobi, SA) and SDP solvers exhibit prohibitive complexity, failing to scale beyond medium-sized instances. Standard spectral clustering (SpecSlap), though fast, yields significantly degraded solutions. "Gap to the best" is calculated relative to the best solution found by any solver in our evaluation.

• Real-world signed networks: Three large social graphs (Epinions, Slashdot, Wikielec) with up to 131K nodes and mixed-sign interactions, representing correlation clustering tasks [Leskovec et al., 2010].
• Synthetic stress tests: Massive scale-free (Barab´asiAlbert) and random (Erd˝os-R´enyi) graphs with mixed weights (up to N ≈ 106 ) generated to validate scalability and robustness to degree heterogeneity.

Baselines We compare DSN against methods spanning three categories: (1) Static spectral (normalized signed Laplacian); (2) Combinatorial metaheuristics: Simulated Annealing (SA), providing high-quality references for solution optimality; (3) Exact solvers (Gurobi [Gurobi Optimization, LLC, 2024], time-limited) for small instances. To ensure strict fairness, all spectral variants utilize identical postrounding refinement, and all heuristic baselines operate under strictly matched wall-clock time budgets.

Metrics The primary metric is the Discrete energy E(s) = -s⊤Js (or Cut value), allowing direct quality comparison across all solver types. To validate the proposed theoretical mechanism, we report the Spectral bound value B(W), quantifying the tightness of the learned relaxation relative to static baselines. We also measure the Optimality gap [E(s)-B(W)]/|E(s)| and Runtime to assess scalability. The bound B(W) ≤ E ⋆ ≤ E(s) further certifies that this gap upper-bounds the per-instance rounding error.

1 SDP and Gurobi averages include zeros from OOM/timeout failures. On the subset where SDP successfully converges, the en-

Implementation details. We implement DSN in PyTorch with sparse matrix support. We set ε = 10-5 for numerical stability and run T = 300 iterations of PGD with initial η = 0.1 and exponential decay (factor 0.99 every 10 steps). Results are averaged over 10 seeds. Profiling on Slashdot gives γt ≥ 2.61× safeguard, Kt ∈ [38, 56] matvecs/step (near-linear in |E|), wi ∈ [0.81, 81.5] (bulk near ¯w=1, std 1.18; hubs receive large weight per the localisation penalty, Eq. (10)), and ∥∇B∥F ∈ [3.61, 4.68], confirming Lemma 4.

Hardware Experiments were conducted on a highperformance computing node with dual Intel Xeon Gold 6246 CPUs and 1.5 TB of RAM. To accelerate spectral computations, we utilized NVIDIA A100 GPUs (40GB VRAM).

5.2 Scalability and Performance Trade-offs

We evaluate DSN against a spectrum of solvers ranging from exact combinatorial methods to fast spectral heuristics (as summarized in Table 1 and visualized in Figure 3).

Performance on standard benchmarks. On the Gset benchmark (66 instances), DSN achieves a mean energy of -5514.96, substantially outperforming both the static spectral baseline SpecSLap (-5090.54) and the SDP relaxation (-3038.54). Notably, DSN attains solution quality within 3% of Simulated Annealing (SA) (-5674.66) while requiring only 3.74 seconds compared to SA's 7.19 seconds. On the converged subset above, DSN comes within 0.55% of Gurobi's certified optimum, but Gurobi hits the 2-hour cutoff on most other Gset instances; its time-limited incumbent (-4573.78 average) consequently lags DSN by over 17%, reflecting the practical scalability ceiling of branch-and-bound at this problem size.

Efficiency on real-world signed networks. DSN demonstrates consistent advantages on large-scale real datasets. For the Epinions graph (131K nodes), DSN obtains an energy of -606858.0, matching SA's solution (-606869.0) to within 0.002% while reducing runtime by a factor of 190× (6.32s vs. 1210s). Similar patterns emerge on Slashdot and Wikielec, where DSN achieves energies of -350960.5 and

-71524.8 respectively, both within 0.1% of SA, at a fraction of the computational cost. In contrast, SDP and Gurobi fail entirely on these large networks due to memory exhaustion and complexity limits.

Scalability tomassive synthetic graphs. Figure 3 presents the central empirical finding of this work: DSN scales to graphs with N = 223 (≈ 8.4 × 106 nodes) while maintaining competitive solution quality.

• On Barab´asi-Albert graphs at scale N = 223 , DSN achieves an energy of -2134.8 × 106 in 856s. SpecSLap fails to produce results beyond 222 . At N = 216 (the limit of SA's scalability), DSN obtains an energy of -17076.4K compared to SA's -17260.0K, a gap of 1.1%, while running 107× faster (11.0s vs. 1177.7s).
• On Erd}os-R´enyi (ER) graphs, DSN shows similar robustness. At N = 223 , DSN obtains an energy of

ergies match theoretical expectations: Gurobi (-4372.15) < SDP (-4359.64) < DSN (-4348.20) < SpecSLap (-4125.70).

-2276.1 × 106 in 998 seconds. Across all comparable scales, DSN consistently captures over 99% of SA's solution quality while requiring only 3-5% of the runtime.

High-precision solvers exhibit limitations in this regime: SDP fails beyond 211 due to memory constraints, while Gurobi cannot scale past 212 within practical time budgets.

These results delineate a fundamental trade-off in largescale Ising optimization. SpecSLap occupies one extreme: sub-second runtimes, but solutions that lag 3 - 15% behind the best achievable quality. SA and Gurobi occupy the opposite extreme: near-optimal solutions but with runtimes that grow prohibitively, rendering them infeasible beyond N ≈ 216 . DSN uniquely bridges this gap, scaling to 223 while maintaining solution quality within 1-3% of metaheuristic solvers at a fraction of their computational cost.

5.3 Ablation Study

We isolate the impact of DSN's design choices via an ablation study on the Wikielec dataset (Table 2). Initializing with node degrees (W (0) = D) proves essential, improving solution quality by 2.0% over random initialization and 0.8% over identity (W = I), confirming that degree-based weighting offers a geometrically informative warm-start. Our gradient (via Hellmann-Feynman) matches Automatic Differentiation (Autodiff) in solution quality while delivering a 6.2× speedup (4.57s vs. 28.33s), by avoiding memory-intensive backpropagation through the iterative eigensolver. The exponential decay schedule yields the lowest energy (-71524.8), outperforming linear decay and fixed rates by 0.5% and 1.5%, respectively, balancing early-stage exploration with late-stage exploitation for fine-grained spectral refinement.

6 Conclusion

We introduced DSN, which treats graph normalization as a learnable diagonal metric to bridge the gap between spectral efficiency and combinatorial precision. Using the HellmannFeynman theorem, DSN achieves exact analytical gradients with near-linear complexity, enabling effective optimization on massive graphs with up to 8.4 × 106 nodes. We established that under a generic spectral-gap condition, the DSN objective is smooth and guaranteed to converge to stationary points. Empirically, DSN consistently outperforms static spectral baselines and matches state-of-the-art metaheuristics within 1-3% while reducing computational cost by up to two orders of magnitude.

Acknowledgments

This work was partially supported by NSF AMPS-2229075, Commonwealth Cyber Initiative (CCI) Awards, and the VCU Quest Award.