snn_opt

A spiking neural network solver for constrained convex optimization.

Abstract

snn_opt is a Python implementation of the spiking neural network (SNN) → convex optimization equivalence, developed as part of an ongoing research program on neuromorphic computation for classical machine-learning problems. It solves quadratic and linear programs of the form

$$ \min_{x \in \mathbb{R}^n}\ \tfrac{1}{2}, x^\top A x + b^\top x \quad\text{subject to}\quad C x + d \le 0, $$

by alternating gradient descent, which plays the role of leaky-integrate membrane drift, with discrete projection events that clamp the trajectory to the constraint boundary, the optimization analogue of an integrate-and-fire spike. The construction follows Mancoo, Boerlin and Machens (NeurIPS 2020) and is the canonical solver underlying the SNN-X research program, a series of classical machine-learning problems recast as constrained convex programs and solved by these dynamics (see Applications).

This repository is intended both as a research artifact, since every published SNN-X paper can be reproduced from the code here, and as a teaching resource for students entering the area: it ships with annotated examples, a self-contained mathematical writeup, and a benchmark suite that visualizes convergence, projection dynamics, and the solver's accuracy limits.

The problem

Given a positive semi-definite Hessian $A \in \mathbb{R}^{n\times n}$, a linear cost $b \in \mathbb{R}^n$, and $m$ linear inequality constraints stacked into $C \in \mathbb{R}^{m \times n}$ and $d \in \mathbb{R}^m$, we seek

$$ x^\star ;=; \arg\min_{x}\ \tfrac{1}{2}, x^\top A x + b^\top x \quad\text{s.t.}\quad c_i^\top x + d_i \le 0,\ i = 1,\dots,m. $$

The class subsumes box-constrained QPs (set $C = [I; -I]$), linear programs ($A = 0$), kernel-ridge subproblems, support-vector machine duals, projected-gradient flows on polytopes, and the bulk of the inner solves that arise in receding-horizon control.

The spiking idea, in one picture

The continuous-time dynamics

$$ \dot x ;=; -\nabla f(x) ;-; C^\top s(t) $$

models a population of $n$ leaky integrators driven by the gradient $\nabla f(x) = Ax + b$, with a corrective spike train $s(t)$ that fires whenever an inequality $c_i^\top x + d_i$ would otherwise become positive. Each spike applies a minimal projection that re-enters the feasible set; spike inter-arrival times encode constraint traffic. Discretized with forward Euler and an adaptive step that reaches the boundary exactly, this becomes a fast projected-gradient solver with diagnostics that double as a neural raster plot.

See docs/theory.md for the full derivation, including the eigenvalue-based step-size choice that eliminates k0 as a hyperparameter and the treatment of bound constraints as implicit facets of the same projection sweep.

Convergence and projection dynamics

Four diagnostic figures, regenerated from benchmarks/ with python benchmarks/run_all.py, give a quick visual sense of what the solver actually does. Every objective gap below is measured against an exact optimum computed by the active-set KKT solve in benchmarks/qpref.py, never against a long run of snn_opt itself; scoring the solver against its own fixed point cannot reveal a standing offset between that fixed point and the true minimiser, and on these problems there is one.

Convergence on a random 50-D QP with 30 inequalities (7 active at the optimum). The gap descends geometrically for about 1800 iterations, then the iterate settles into a period-2 limit cycle: it alternates between two points whose gaps are 3.0e-4 and 6.7e-4, which is why the two branches in panel (a) are drawn separately and why panel (b) flatlines rather than decaying. The run is jointly feasible throughout (max row distance 8.2e-7) and reports converged=False at the 4000-iteration cap: the scale-invariant KKT certificate (v0.6.0) measures the limit cycle's relative optimality defect at 8.9e-4, above the default kkt_rel_tol = 1e-4, and the value is a true floor (unchanged at 40k iterations). Loosen kkt_rel_tol, or shrink the floor itself with k0_scale, to certify on this problem. See Accuracy and tuning for where the floor comes from.

Projection-spike raster on an 8-D QP over a 16-facet polytope whose unconstrained minimiser sits well outside the feasible set. Each marker is one projection event. The network visibly searches for the active set: row 11 fires a burst around t = 4..10 and then falls silent, while rows 5, 10 and 3 are recruited at t ~= 21, 23 and 41 and fire on every step thereafter. Those three persistent rows are exactly the active set of the true optimum. This is the practical payoff of the spiking view, and the literal sense in which the solver is spiking.

Warm-start speedup on a sequence of 30 drifting QPs, a stylized MPC workload, measured under the v0.6.0 KKT stopping criterion. From the second problem onward, warm starting cuts a 221-iteration cold solve to 101 iterations, an essentially free 2.19x, and wall time falls in step (2.03x). Iterations are the headline because they are deterministic; the wall-time panel is the median of five timed runs per problem, since a single pass picks up scheduler noise indistinguishable from signal.

Accuracy against step size, at three iteration budgets. The spiking dynamics converge to a fixed point of the discretised flow, which is offset from the exact minimiser by an amount that shrinks with the gradient step k0. Smaller k0 also means more iterations are needed to arrive, so each budget has a knee, and the knee moves left as the budget grows. The shipped default (k0_scale = 0.5) sits to the right of every knee: on this problem it leaves a gap of 6.7e-4, while k0_scale = 0.02 reaches 1.2e-5 given 80k iterations.

A closer look: trajectory and raw-vs-optimized modes

Two figures generated by the example scripts give a more concrete sense of what the dynamics look like in 2-D, where everything is easy to visualize:

State evolution and 2-D trajectory for examples/example1_basic_2d.py, a constrained QP whose unconstrained minimum lies outside the feasible polytope. Left: per-component value over time, with the projection events shown as a rug along the bottom. The sawtooth in each trace is the integrate-and-fire dynamic itself: drift away from the boundary, spike back onto it, repeat. Right: the trajectory in state space over objective contours, gliding down the gradient, meeting the active facet and sliding along it to the constrained optimum. This example deliberately runs projection_method='fixed', which produces many small corrections rather than one exact jump, because that is what makes the spiking behaviour visible.

Output of examples/example_raw_mode.py, comparing a fixed gradient and projection step against the defaults (auto k0 from the Lipschitz constant, exact projection to the boundary). Both reach the optimum on a problem whose constraint is active at the solution: the optimized run gets within 1e-16 of the exact objective in roughly 25 iterations, while the raw run needs the full 300 to reach 1e-15.

Installation

snn_opt requires Python 3.9+, NumPy, and SciPy. The fastest path is PyPI:

pip install snn-opt                # core, prebuilt wheel (no compiler needed)
pip install "snn-opt[examples]"    # also installs matplotlib for examples
pip install "snn-opt[dev]"         # examples + cvxpy + pytest + ruff

The PyPI distribution name is snn-opt (hyphenated, lowercase, per PEP 503); the Python import name is snn_opt. So you pip install snn-opt and then import snn_opt.

For an editable install from a checkout (development workflow):

git clone https://github.com/ahkhan03/SNN_opt.git
cd SNN_opt
pip install -e .                   # core
pip install -e ".[examples]"       # also installs matplotlib for examples
pip install -e ".[dev]"            # examples + cvxpy + pytest + ruff

For a specific commit (reproducibility for papers/collaborators):

pip install "git+https://github.com/ahkhan03/SNN_opt.git@<commit-sha>"

The package can also be run without installation: every example and test sits next to a small sys.path bootstrap that points at src/. Smoke test:

python tests/test_installation.py

Compiled C++ backend

The PyPI wheels ship a precompiled C++ kernel (snn_opt._kernel) that accelerates the inner adaptive-projection loop by roughly an order of magnitude over the pure-Python path. Opt in via the backend keyword:

result = solve_qp(A, b, C, d, x0, backend='c')        # compiled kernel (auto)
result = solve_qp(A, b, C, d, x0, backend='python')   # reference (default)

The compiled kernel comes in three numerically identical variants that differ only in how the inner matrix–vector products are threaded:

Only the matvec is data-parallel; the Euler recurrence and the greedy projection are inherently serial (an Amdahl ceiling of roughly 2–3× on a few cores). Because per-call thread fork/join only pays off on large systems, the multicore path is automatically skipped below a work threshold, so 'c' matches the serial path on small/medium problems and only spins up threads on large ones. Multithreading honours OMP_NUM_THREADS; snn_opt._kernel.HAS_OPENMP and snn_opt._kernel.max_threads() report the build's capability.

All backends are kept in lockstep by the parity test suite (tests/test_c_backend_parity.py). The C kernel supports dense problems with projection_method='adaptive'; sparse and non-adaptive paths transparently use Python. The same kernel source is HLS-compatible and is the basis for the FPGA deployment track. When the precompiled kernel is unavailable on your platform (rare), the 'c*' backends raise a clear error and the Python backend continues to work.

Problem transforms (eigenbasis)

Orthogonal to the backend, the transform axis rewrites the problem into an equivalent one that is cheaper to solve and maps the solution back. Transforms are an explicit opt-in (SolverConfig.transform); the canonical solver stays the default, and a transform composes with any backend.

from snn_opt import solve_qp
result = solve_qp(A, b, C, d, x0, ...)  # canonical (default)

from snn_opt import SNNSolver, SolverConfig, OptimizationProblem
cfg = SolverConfig(transform='eigenbasis', backend='c')
result = SNNSolver(OptimizationProblem(A, b, C, d), cfg).solve(x0)

EigenbasisTransform (transform='eigenbasis') rotates a symmetric-PSD Hessian into its eigenbasis (A = VΛVᵀ), collapsing the dominant O(n²) A @ x gradient step into an O(n) elementwise product; the projection is unchanged because the constraint Gram is rotation-invariant. The win grows with problem size.

Since v0.5.0 the transform does accept box bounds. Per-coordinate bounds are not rotation-invariant, so they cannot stay implicit: they are materialized as explicit rotated unit-norm rows, growing m by up to 2n. The O(1) implicit-facet advantage is deliberately surrendered under a transform, which is the trade to be aware of when combining the two. See docs/api.md.

What v0.5.0 changed

v0.5.0 is a structural correctness release, and the behaviour it fixes is worth understanding before relying on results from an earlier version.

Before v0.5.0, bound constraints were enforced by a terminal clip applied after the halfspace sweep, with nothing re-projecting behind it. Composing the two is not a projection onto their intersection (the classical POCS failure), so on a problem where a bound and an interacting row are simultaneously active, the solver could stall at a point feasible for neither and report an objective that undercuts the true optimum. Bounds are now implicit unit-normal facets inside one unified projection sweep, and the terminal clip is gone.

Three result fields expose feasibility, optimality diagnostics, and projection termination on any nontrivial problem:

result.joint_feasible            # feasibility of rows AND bounds together
result.kkt_residual              # scale-invariant KKT certificate (v0.6.0)
result.projection_budget_exhausted

• joint_feasible is the honest feasibility flag. Pre-0.5 the convergence gate looked at rows only, so a box violation could not fail it.
• kkt_residual is the scale-invariant KKT certificate at the final point (see the next section); with the default settings, converged=True means exactly that this certificate passed, together with feasibility and the plateau criterion, at three consecutive checkpoints. The older stationarity_residual diagnostic is retained for one compatibility release but mixes units and can depend on constraint row order; prefer kkt_residual.
• projection_budget_exhausted reports that the sweep hit its watchdog. max_projection_iters is now a safety cap (default None, auto-sized), and hitting it aborts the solve rather than being reported as convergence.

projection_method='fixed' combined with bounds now raises, because the legacy fixed-step path cannot enforce bounds correctly without the clip that was removed.

Scale-invariant convergence certification (v0.6.0)

Before v0.6.0, converged required an absolute projected-gradient norm below 1e-6. That test had two structural defects, found when an MPC user ran QPs whose gradient scale is ~1e10: (a) rescaling the objective rescales every gradient, so on large-scale problems the flag could never fire at any solution quality; and (b) the projected-gradient heuristic removes each active facet's gradient component independently, so at a constrained optimum with correlated active normals it stalls at a cross-term residue and is structurally nonzero even at the exact optimum. converged=False therefore said nothing about solution quality; solves on perfectly solvable problems ran to their iteration cap by construction.

The v0.6.0 criterion is a KKT-cone certificate: one nonnegative least-squares fit of -∇f(x) onto the cone of all unit-normalized facet normals, augmented with a complementarity row so slack facets cannot absorb the gradient, accepted when

r_kkt  <=  kkt_abs_tol + kkt_rel_tol * max(‖A x‖, ‖b‖, ‖Nᵀμ‖)

Both sides carry gradient units, so while the relative term dominates the threshold the decision is invariant under positive objective rescaling, constraint row order, row duplication, and per-row scaling: the same problem certifies identically at natural scale and at 1e10x. (The kkt_abs_tol floor deliberately takes over at near-zero gradient scales, the intentional fallback that lets a genuinely-zero problem terminate.) The fit runs host-side on every backend: the compiled kernel advances the dynamics in checkpoint-sized chunks and the same Python policy evaluates each checkpoint, so converged means one thing everywhere (the FPGA reference is unchanged and reports fixed-horizon results, which the host can certify with the same function). The cheap plateau/feasibility gates are evaluated first, so the certificate's NNLS cost is confined to near-termination checkpoints, and end-to-end overhead is negligible (see benchmarks/).

Migration: results from v0.5 remain reproducible with ConvergenceConfig(optimality_test="legacy_projected_gradient"), which preserves the old test verbatim. converged=False runs from v0.5 can legitimately become converged=True (or stop ~100x earlier) under the new criterion; nothing about the dynamics changed, only the stopping decision. The default kkt_rel_tol = 1e-4 is calibrated to the O(k0) fixed-point floor of the default step size: it certifies the quality the dynamics genuinely reach, roughly 1e-3 relative solution error on well-conditioned problems. It is a residual tolerance, not an error bound: on a nearly singular Hessian a small residual can coexist with a larger solution displacement.

KV260 reference implementation

fpga/kv260_v05/ contains the restricted fixed-horizon Kria K26 reference that was physically qualified for the SNN-MSRP study. It preserves the v0.5 unified-projection semantics, binary64 input boundary, fixed-point contract, 200 MHz build recipe, and exact compatibility ABI used by that experiment. Its README records the measured qualification surface and the unsupported cases. It is not a general FPGA backend for solve_qp.

Accuracy and tuning

The spiking dynamics converge to a fixed point of the discretised flow, not to the exact minimiser of the QP. On well-conditioned problems the two are close; they are not identical, and the difference is set by the gradient step size k0 = k0_scale / L.

Concretely, on the 50-D benchmark of Figure 1 with the shipped defaults, the solver reaches a period-2 limit cycle whose objective gap against the exact optimum alternates between 3.0e-4 and 6.7e-4, and stays there: the value is identical to ten significant figures at 20k and at 100k iterations. It is jointly feasible the whole time. So the limitation is accuracy of the fixed point, not feasibility.

What to do about it, in order of usefulness:

1. Read kkt_residual / kkt_scale as the optimality verdict. It is the scale-invariant KKT defect of the final point, comparable across problems and objective scalings; converged=True certifies it below kkt_rel_tol. On this benchmark it reports 8.9e-4, an honest measurement of the limit cycle, which is why the run does not certify at the default 1e-4.
2. Lower k0_scale, and raise the iteration budget with it. Figure 4 maps the trade. On that problem, 0.5 gives 6.7e-4 and 0.02 gives 1.2e-5, but only if the budget is large enough to arrive; at 5k iterations the same 0.02 setting is far worse than the default. Tune the pair, never k0_scale alone.
3. Polish externally if you need machine precision. Once the active set is correct (and it usually is, see Figure 2), the exact optimum follows from one equality-constrained KKT solve on those rows. That is exactly what benchmarks/qpref.py does, in well under a millisecond on these sizes.

Two known limitations are worth stating plainly. Ill-conditioned or stiff QPs are the harder case: the native adaptive stepping can fail to reach tolerance and return an infeasible point, and naive Jacobi/diagonal preconditioning conflicts with the adaptive step-size rule rather than fixing it. When a run reports converged=False, read the diagnosis in order: check joint_feasible and projection_budget_exhausted first (feasibility failures and watchdog aborts are their own categories), then read kkt_residual / kkt_scale; since v0.6.0 that number is scale-invariant, so "how far from optimal" is finally a well-posed question at any problem scaling.

Quick start

import numpy as np
from snn_opt import solve_qp

# Minimize ||x||^2 subject to  x_1 + 2 x_2 <= 1  (and that's it).
A  = np.eye(2)
b  = np.zeros(2)
C  = np.array([[1.0, 2.0]])
d  = np.array([-1.0])
x0 = np.array([1.0, 1.0])

result = solve_qp(A, b, C, d, x0, max_iterations=1000)

print(result.summary())             # converged?  iterations?  spikes?
print("x* =", result.final_x)
print("f* =", result.final_objective)

For repeated solves (warm-started receding-horizon problems), construct an SNNSolver once and call .solve(x0) per problem instance; see examples/example4_warm_start.py.

Examples

All scripts live under examples/ and are runnable as plain python examples/example_name.py.

Run them all in sequence:

python examples/run_all_examples.py

Documentation

• docs/theory.md: derivation of the SNN/convex-optimization equivalence, step-size analysis, projection geometry, convergence criteria.
• docs/applications.md: catalogue of the published work that uses this solver.
• docs/api.md: hand-curated API reference.
• benchmarks/README.md: what each figure shows and how to regenerate it.
• https://snn.ahkhan.me: companion site, designed for a broader audience (students, curious researchers).

Applications

The framework is currently demonstrated in:

• Khan, Mohammed & Li (2025), Portfolio Optimization: A Neurodynamic Approach Based on Spiking Neural Networks, Biomimetics, 10(12):808. doi:10.3390/biomimetics10120808. Portfolio selection cast as a constrained QP and solved by the spiking dynamics implemented here.

Additional applications are in preparation. As they reach publication they will be added to docs/applications.md.

Citing this work

If snn_opt plays a role in your research or teaching, please cite both the software and the framework paper:

@software{khan2026snnopt,
  author  = {Khan, Ameer Hamza and Li, Shuai},
  title   = {snn\_opt: A Spiking Neural Network Solver for Constrained Convex Optimization},
  year    = {2026},
  version = {0.6.0},
  url     = {https://github.com/ahkhan03/SNN_opt},
  license = {Apache-2.0},
}

@inproceedings{mancoo2020understanding,
  author    = {Mancoo, Allan and Boerlin, Martin and Machens, Christian K.},
  title     = {Understanding Spiking Networks Through Convex Optimization},
  booktitle = {Advances in Neural Information Processing Systems (NeurIPS)},
  year      = {2020},
}

The full per-paper bibliography of the SNN-X series is maintained at docs/applications.md.

License

Apache-2.0, see LICENSE. Permissive, with an explicit patent grant; suitable for both academic and commercial reuse.

Acknowledgments

Developed at the School of Artificial Intelligence, Taizhou University.

This codebase implements the SNN-QP research program led by Prof. Shuai Li (IEEE Fellow; Faculty of Information Technology and Electrical Engineering, University of Oulu, Finland), whose work on neurodynamic optimization originated this line of inquiry. The mathematical framework follows Mancoo, Boerlin and Machens (NeurIPS 2020) and the broader projection-neural-network lineage (Hopfield–Tank, Kennedy–Chua, Xia–Wang, Liu–Wang). Pull requests, bug reports, and citations of the SNN-X papers in your own work are all warmly welcomed.