snn-opt 0.4.0

A spiking neural network solver for constrained convex optimization (QP/LP).

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 — playing 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 publication series (PCA, Ridge, TDSVM, Norm, SVM, CF, KRR, Procrustes — see Applications).

This repository is intended both as a research artifact — 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 and projection dynamics.

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 box-clipping shortcut for problems like SVM.

Convergence and projection dynamics

Three diagnostic figures, regenerated from benchmarks/, give a quick visual sense of what the solver actually does:

Convergence on a random 50-D QP with 30 inequalities. The objective gap drops geometrically over a few thousand iterations; the iterate stability ‖x_{t+1}−x_t‖ mirrors it; constraint violation stays at machine-precision floor throughout.
Projection-spike raster on a 4-D box-constrained problem whose unconstrained optimum lies outside the box. Each row is one inequality, each marker a spike (sized by displacement); only the four "active" faces of the box fire. The bottom panel is the corresponding objective gap. This is the literal sense in which the solver is spiking.
Warm-start speedup on a sequence of 30 drifting QPs (a stylized MPC workload). Cold-started solves take ~260 iterations each; warm-started solves drop to ~140 from the second problem onward — an essentially-free 1.8× speedup, the property that makes this dynamic well-suited to receding-horizon problems.

Reproduce these from a checkout with python benchmarks/run_all.py.

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 orange spike markers at projection events. Right: the trajectory in state space superimposed on objective contours; the trajectory glides down the gradient, hits the active facet, and slides along it to the constrained optimum.

Output of examples/example_raw_mode.py. Left two panels: the trajectory under "raw" (no auto step size, no adaptive projection) and "optimized" defaults — both reach the optimum, but the optimized run does so in 201 iterations with a single projection event versus 300 / 5. Right: convergence on a log-objective scale; the optimized run gains roughly 30 orders of magnitude per 100 iterations once it leaves the boundary.

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:

backend	matvec threading
'c'	auto — OpenMP multicore when the wheel was built with it, else single-thread
'c_serial'	forced single-thread (SIMD only)
'c_openmp'	forced OpenMP multicore (raises if the wheel was built without OpenMP)

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. Box constraints are not supported with this transform (per-coordinate bounds are not rotation-invariant) — it raises a clear error if lower_bound/upper_bound are set. See docs/api.md.

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.

#	Script	Problem	Highlights
1	example1_simple_2d.py	2D quadratic with two linear cuts	Smallest possible runnable demo
1b	example1_basic_2d.py	Same problem, with trajectory plot	See examples/example1_basic_2d.png
1c	example1_advanced_2d.py	Shifted feasible region, infeasible start	Spike raster + violation plot
2	example2_3d_polytope.py	3D QP with 4 hyperplanes	Multiple active constraints, vertex solution
3	example3_linear_program.py	Box-constrained LP ($A=0$)	LP via the same machinery
4	example4_warm_start.py	Sequence of related QPs	Receding-horizon / MPC pattern; spikes drop $30 \to 0$
5	example5_infeasible_recovery.py	Infeasible initializations	Automatic projection to feasibility
6	example6_equality_constraint.py	Equality via a sandwiched band	$x_1 = a$ as a tight $\pm \varepsilon$ inequality pair
7	example7_svm_dual.py	SVM dual with kernel	Box clipping + auto step size on a real ML task
—	example_raw_mode.py	Bypass auto-config	Compares raw vs. optimized solver settings

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 — one-page summary of each SNN-X paper with links to PDFs.
• docs/api.md — hand-curated API reference.
• 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.4.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.