HOPSO — Harmonic Oscillator based Particle Swarm Optimization

Reference implementation of the optimizer introduced in Harmonic Oscillator based Particle Swarm Optimization by Yury Chernyak, Ijaz Ahamed Mohammad, Nikolas Masnicak, Matej Pivoluska, and Martin Plesch.

---

The idea

Standard PSO moves particles by a velocity update rule with no intrinsic notion of energy, which lets velocities grow without bound and makes convergence hard to control. HOPSO instead treats each particle as a damped harmonic oscillator swinging about an attractor placed between its personal best and the swarm's global best:

x(t) = A · e^(−λt) · cos(ωt + θ) + attractor

The amplitude A decays exponentially, so the swarm contracts on its own — exploration early, exploitation late — without a hand-tuned inertia schedule. A floor on A, proportional to the personal-best/global-best separation, keeps particles from collapsing onto the attractor prematurely. When a particle improves, its clock resets to t = 0 and it swings wide again.

Install and run

git clone https://github.com/hbar-systems/hbar.science-hopso.git
cd hbar.science-hopso
pip install numpy tqdm
pytest tests/ -q      # 8 tests, ~1 second

Minimal use — results come back through output lists, one entry per run:

import numpy as np
from hopso import hopso

def rastrigin(x):
    z = np.asarray(x) - 1.5
    return 10*len(z) + np.sum(z**2 - 10*np.cos(2*np.pi*z))

# hp = [w1, w2, tm, lambda]
e_min, vectors, velocities, vel_mag, gbest, amps, pos = [], [], [], [], [], [], []
hopso(rastrigin, [1.0, 1.0, 2*np.pi, 0.05],
      num_particles=12, runs=10, dimension=6, max_cut=2.0,
      e_min=e_min, vectors=vectors, velocities=velocities, vel_mag=vel_mag,
      gbest=gbest, amps=amps, pos=pos, max_iterations=150)

print(min(e_min))          # best energy found across runs
print(vectors[np.argmin(e_min)])   # the corresponding parameter vector

Search bounds are hard-coded to [0, π] in both implementations — initial positions are sampled there and positions are clipped there each iteration. The source marks these lines # CHANGE BOUNDS ACCORDINGLY. Rescale your objective, or edit those lines, before using HOPSO on a different domain.

What is here

path	contents
hopso.py	Reference implementation, as published. Standard Euclidean search space.
hopso_periodicity.py	Variant for periodic parameters (quantum circuit angles): circular distance and wrapping, immediate per-particle attractor updates.
examples/vqe_ansatz_demo.py	HOPSO driving a VQE ansatz — circuit construction, hopping terms, energy objective.
hpc/	Distributed MPI version with SLURM submission scripts. See hpc/README.md.
tests/	Behavioural tests pinning the published behaviour.

Parameters

name	meaning
hp[0], hp[1]	weights w1, w2 on personal best and global best when placing the attractor
hp[2]	tm — scale of the random time increment per iteration
hp[3]	λ — damping rate of the amplitude envelope
max_cut	multiplier on the amplitude floor derived from personal/global best separation
num_particles	swarm size
runs	independent restarts

---

A note on np.maximum

Two lines in hopso.py (in the amplitude-recalculation branches) call:

A[i] = np.maximum(A[i], A1[i], amp_dist[i])

NumPy interprets a third positional argument to np.maximum as the output buffer, not as a third operand. This therefore computes max(A[i], A1[i]) and writes the result into amp_dist[i]; the amp_dist floor does not enter the maximum. Every later variant in this repository — hopso_periodicity.py, examples/vqe_ansatz_demo.py, and both files in hpc/ — uses the nested form np.maximum(np.maximum(A, A1), a_dist) instead.

hopso.py is kept exactly as published. We measured the difference: rewriting those two lines to a genuine three-way maximum changes results materially and does not improve them (on a 6-dimensional Rastrigin benchmark the published form produced a better mean best-energy on two of three seeds). Since it is not an improvement and the published results were generated with this code, changing it would cost reproducibility and buy nothing.

Two consequences for users:

• To reproduce the paper, use hopso.py unmodified.
• NumPy has deprecated passing more than two positional arguments to np.maximum, so these lines will eventually raise. When that happens the correct minimal fix is np.maximum(A[i], A1[i], out=amp_dist[i]), which preserves the current behaviour exactly. Do not "fix" it to a three-way maximum without re-validating the paper's numbers.

---

Related

The VQE regularization study that uses HOPSO as one of its optimizers: hbar.science-vqe-regularization.

Authors and citation

HOPSO is joint work — see CITATION.cff for the full author list and arXiv:2410.08043 for the paper. The MPI and SLURM code in hpc/ was contributed by Ijaz Ahamed Mohammad.

License

MIT — see LICENSE.