@euriklis/quadprog

quadprog in TypeScript — a typed, dependency‑free quadratic‑programming (Goldfarb–Idnani) solver for Node and Bun.

A small, dependency‑free convex quadratic‑programming solver for JavaScript and TypeScript (Node and Bun). Written in TypeScript, it ships type declarations (.d.ts) — so solveQP, solveQPFast, and the QPResult shape are fully typed out of the box. It implements the Goldfarb–Idnani dual active‑set method and ships two entry points:

• solveQP — the pure‑scalar solver. Zero dependencies, deterministic, fast.
• solveQPFast — the same algorithm, but the cubic matrix factorization in the start‑up phase runs on a SharedArrayBuffer worker pool. Pays off for large dense problems (n ≳ 512); below that it transparently calls solveQP.

Both return bit‑for‑bit identical results.

It solves

minimize    ½ xᵀ D x − dᵀ x          x ∈ ℝⁿ
subject to  A₁ᵀ x  =  b₁             (the first  meq  constraints)
            A₂ᵀ x  ≥  b₂             (the remaining ones)

with D symmetric positive‑definite. This covers portfolio optimization with constraints, constrained least squares, the dual problem of support‑vector machines, and RBF networks, among many others.

Installation

npm install @euriklis/quadprog@latest --save

Usage

A is an n × q matrix whose column i is the normal of constraint i (so A[variable][constraint]). b has length q, d length n, D is n × n. meq (default 0) is the number of leading equality constraints.

import { solveQP } from "@euriklis/quadprog";

// minimize ½‖x‖² − dᵀx  s.t.  the three inequalities Aᵀx ≥ b
const D = [[1, 0, 0], [0, 1, 0], [0, 0, 1]];
const d = [0, 5, 0];
const A = [[-4, 2, 0], [-3, 1, -2], [0, 0, 1]];
const b = [-8, 2, 0];

const r = solveQP(D, d, A, b);
// r.solution                  → [0.476190, 1.047619, 2.095238]
// r.value                     → −2.380952   (½xᵀDx − dᵀx at the optimum)
// r["Lagrangian multipliers"] → [0, 0.238095, 2.095238]
// r["active constraints"]     → [2, 1, 0]   (first two are active)
// r["count of active constraints"] → 2

For large dense problems, solveQPFast is a drop‑in async replacement:

import { solveQPFast } from "@euriklis/quadprog";

const r = await solveQPFast(D, d, A, b);   // n < 512 → delegates to solveQP

Return value

| field | meaning | |---|---| | solution | the minimizer x* (length n) | | value | the objective ½x*ᵀDx* − dᵀx* | | unconstrained_solution | D⁻¹d, the minimizer ignoring all constraints | | Lagrangian multipliers | one per constraint, 0 if inactive | | active constraints / count of active constraints | which constraints bind at x* | | iterations | [main iterations, constraints dropped] | | ierr | exit code: 0 success, 1 infeasible (inconsistent constraints), 2 D not positive‑definite | | message | human‑readable form of ierr |

Check ierr before trusting solution. On an infeasible problem (ierr === 1) or a non‑SPD D (ierr === 2) the returned solution does not solve the program. Branch on the numeric ierr rather than matching message.

const r = solveQP(D, d, A, b);
if (r.ierr !== 0) throw new Error(r.message);

Examples

Runnable scripts live in examples/ (run with node or bun):

| file | what it shows | |---|---| | 01-basic.js | a plain inequality‑constrained QP | | 02-equality.js | an equality constraint via meq (projection onto a line) | | 03-portfolio.js | long‑only minimum‑variance portfolio (budget equality + no‑short bounds) | | 04-fast-large.js | solveQPFast on a large dense problem, vs solveQP |

node examples/03-portfolio.js
# weights  : [ '0.4530', '0.1693', '0.2673', '0.1105' ]
# Σ weights: 1.000000 (= 1, budget)
# all ≥ 0  : true
# variance : 0.020424

How it compares

The two established options on npm for QP in JavaScript are Alberto Santini's quadprog (the long‑standing 1‑indexed port of the same Goldfarb–Idnani Fortran) and osqp (the OSQP first‑order solver compiled to WebAssembly). They solve different regimes — pick by problem shape:

| | @euriklis/quadprog | quadprog (Santini) | osqp (WASM) | |---|---|---|---| | Method | dual active‑set (Goldfarb–Idnani) | dual active‑set (same) | first‑order (ADMM) | | Best for | small–medium dense QP, exact | small dense QP | large sparse QP | | Accuracy | machine‑precision (direct) | machine‑precision (direct) | tolerance‑based, needs tuning | | Language / types | TypeScript, ships .d.ts | JS, no types | has types (WASM bindings) | | Module | ESM (Node + Bun) | CommonJS | ESM + WASM | | API indexing | 0‑indexed (idiomatic JS) | 1‑indexed (arrays padded) | matrix triplets | | Runtime deps | none | none | WASM binary + async init | | Parallel start‑up | solveQPFast (worker pool) | — | — (single WASM thread) | | Infeasibility | flagged via ierr | flagged via ierr | flagged via status | | License | MIT | MIT | Apache‑2.0 |

Which should you use?

• Dense D, exact answer, no build step → this package. Same algorithm as the reference, but modern (ESM + types + 0‑indexed), validated value‑for‑value against it across 864 randomised problems, with an optional parallel factorization for large n.
• Very large, sparse constraint matrices → osqp. A first‑order method with sparse linear algebra scales past what a dense active‑set solver targets (at the cost of an approximate, tolerance‑controlled solution and a WASM dependency).
• Already on Santini's quadprog and happy → this is a drop‑in with the same results; you gain types, ESM, 0‑indexing, and solveQPFast.

The algorithm

This section explains why the method works. It assumes you are comfortable with positive‑definite matrices, the Cholesky factorization, QR / Givens rotations, and Lagrange multipliers.

New to the active set or the dual problem? Both get a plain‑language, picture‑first primer in docs/concepts.md — read that first if either term is unfamiliar.

1. The problem is a single point

Write the objective f(x) = ½ xᵀ D x − dᵀ x. Because D is symmetric positive‑definite, f is strictly convex, the feasible region

Ω = { x : aᵢᵀx ≥ bᵢ  (i ≥ meq),  aᵢᵀx = bᵢ  (i < meq) }

is a convex polyhedron, and the constrained minimum — if Ω ≠ ∅ — is the unique global one. There is no question of local minima.

Its gradient is ∇f(x) = D x − d. Setting it to zero gives the unconstrained minimum

x⁰ = D⁻¹ d.

If x⁰ ∈ Ω we are already done. The interesting case is when x⁰ violates some constraints and the optimum sits on the boundary of Ω.

2. What the optimum looks like — the KKT conditions

A feasible x* is the optimum iff there exist multipliers λᵢ ≥ 0 such that

   D x* − d = Σ λᵢ aᵢ          (stationarity: the gradient is a combination of
                                active constraint normals)
   λᵢ ( aᵢᵀx* − bᵢ ) = 0       (complementary slackness: λᵢ = 0 unless
                                constraint i is tight)
   λᵢ ≥ 0  for inequalities    (dual feasibility)

Because the problem is convex, these conditions are not just necessary but sufficient. So the whole job is: find the set of constraints that are tight at the optimum — the active set 𝒜 — together with multipliers λ ≥ 0 that balance the gradient. Once 𝒜 is known the solution is pure linear algebra: minimize f subject to aᵢᵀx = bᵢ for i ∈ 𝒜.

3. Two ways to search, and why the dual one is used

A primal active‑set method keeps x feasible and walks along the boundary, swapping constraints in and out until the multipliers come out non‑negative. It needs a feasible starting point, which itself costs a solve.

Goldfarb and Idnani turn this around. The dual method keeps the multiplier side healthy (λ ≥ 0) and works towards feasibility:

Start at the unconstrained minimum x⁰ = D⁻¹d with an empty active set (λ = 0, trivially dual‑feasible) and, one constraint at a time, repair the worst violation while never letting any λᵢ go negative.

Its great advantage: no feasible starting point is needed — x⁰ is free — and every iterate is already the exact optimum of the QP restricted to the current active set.

4. The right coordinate system: J = R⁻¹

Factor the (constant) matrix once,

D = Rᵀ R         (Cholesky, R upper‑triangular),     J := R⁻¹.

Then D⁻¹ = J Jᵀ and, crucially, Jᵀ D J = I: the columns of J form an orthonormal basis in the inner product defined by D. Working in this basis turns the awkward D‑weighted geometry into ordinary Euclidean geometry.

While the algorithm runs it maintains, for the current active set of size k, a QR‑type factorization of the active constraint normals N = [ a_{𝒜(1)} , … , a_{𝒜(k)} ] expressed in the J‑basis. Split the transformed basis as J = (Q₁ | Q₂), where the first k columns Q₁ span the active normals and the remaining Q₂ span the directions still free to move. Adding or dropping a constraint is then one sweep of Givens rotations that re‑triangularises this factor in O(n²), instead of refactorizing from scratch in O(n³).

5. One iteration

Let x be the current point with active set 𝒜 and multipliers u ≥ 0.

1. Pick the most‑violated constraint. Among all constraints compute the (norm‑scaled) violation aₚᵀx − bₚ; let p be the most negative. If none is negative, every KKT condition holds — stop, x is the optimum.

2. Step direction. With n⁺ = aₚ the violated normal, split it in the J‑basis into the part lying in the free subspace and the part in the active subspace:

   z = Q₂ Q₂ᵀ J n⁺      (primal direction: how x should move)
   r = R⁻¹ Q₁ᵀ n⁺       (dual direction:  how the active λ react)

   z is the direction in which moving x reduces the violation of p without disturbing the already‑active constraints; r says how the current multipliers change as we move.

3. Step length t = min(t₁, t₂):

   • t₂ (primal) — the distance along z until constraint p becomes exactly tight, t₂ = −(aₚᵀx − bₚ) / (zᵀ aₚ).
   • t₁ (dual) — the largest step before some active multiplier would turn negative, t₁ = min_{ i : rᵢ > 0 } uᵢ / rᵢ. The minimizer it₁ is the constraint that would be driven infeasible in the dual.

4. Take the step.

   • Full step (t₂ ≤ t₁, and finite): move x ← x + t₂ z, update the multipliers, and add p to the active set (a Givens update of the factorization). The active set grew by one; go to 1.
   • Partial step (t₁ < t₂): we cannot reach tightness yet — the blocking constraint it₁ must leave first. Update u ← u − t₁ r, drop it₁ (a Givens down‑date), and recompute the direction for the same p.
   • If z = 0 and no dual step is possible, the constraints are inconsistent — Ω = ∅, report infeasible.

5. Repeat. Each main iteration either makes a constraint active for good or removes one that should not have been; the dual objective increases monotonically, so the process is finite.

When the loop ends, the active multipliers are read off from u, every inactive constraint gets λ = 0, and x is returned as x*.

6. Cost, and what solveQPFast parallelizes

Per iteration the work is O(n²) (the Givens sweeps and matrix–vector products); there are typically O(n + q) iterations. The one‑off start‑up, however, is O(n³): the Cholesky D = RᵀR and the triangular inverse J = R⁻¹.

The active‑set loop is inherently sequential, but that start‑up is not. solveQPFast keeps the identical loop and only accelerates the factorization:

• a blocked Cholesky whose trailing‑matrix update A₂₂ ← A₂₂ − L₂₁L₂₁ᵀ is a matrix multiply, and
• a blocked triangular inverse, R = [[A,B],[0,C]] ⇒ R⁻¹ = [[A⁻¹, −A⁻¹BC⁻¹],[0, C⁻¹]], whose only cubic term −A⁻¹BC⁻¹ is two matrix multiplies,

with every matrix multiply dispatched across a SharedArrayBuffer worker pool. Below n = 512 the dispatch overhead is not worth it and solveQPFast simply calls solveQP.

The pool's workers are unref'd, so they never keep the process alive — a script that awaits solveQPFast exits on its own once it finishes; no process.exit is needed. A long‑running host that wants to release the threads eagerly can import { shutdown } from "@euriklis/quadprog" and call shutdown().

Correctness

Every result is validated, value for value, against Alberto Santini's quadprog — the long‑standing, heavily‑used 1‑indexed translation of Berwin Turlach's reference Fortran qpgen2 — across 864 randomised problems (square and rectangular, n = 1…12, q = 1…24, with and without equality constraints), plus the analytic edge cases (a single active bound, an equality‑constrained projection, a fully‑determined active set). Solutions and Lagrange multipliers agree to machine precision (≈ 10⁻¹⁵). solveQPFast is checked the same way for n ∈ {512, 768, 1024}, including feasible constraint‑heavy problems that run the active‑set loop hundreds of times through the parallel factorization path. Run the suite with npm test.

Feasibility caveat — solveQPFast. On a feasible problem solveQPFast equals solveQP to machine precision. On an infeasible one (no x satisfies the constraints) there is no meaningful answer: every solver returns a constraint‑violating point, and because the parallel and scalar factorizations differ in their last bits, a long chaotic active‑set path can drive those two garbage outputs apart. That is garbage‑in/garbage‑out, not a disagreement on a real solution — check min(Aᵀx − b) ≥ 0 if in doubt.

Note for users of earlier versions: the previous 0‑indexed port mistranslated several packed‑storage offsets from the 1‑based Fortran (a sentinel collision that prevented constraint 0 from ever activating, two meq off‑by‑ones, the triangular‑inverse strides, a Givens drop stride and loop bound, and the constraint count for non‑square A). These are fixed; the solver now matches the reference everywhere.

Performance

All numbers are on an AMD Ryzen 9 4900HS (16 threads), Node v20.19 / Bun 1.2.23, median ms; reproduce with the scripts in benchmark/.

How solveQP compares to Santini depends on the problem. Both share the same O(n³) factorization but differ in the O(n²)‑per‑iteration active‑set loop, so the picture splits in two:

(a) Constraint‑heavy problems — many active constraints ⇒ many iterations ⇒ the loop dominates, and the dense 0‑indexed, goto‑free layout pulls ahead (node benchmark/constraints.js, feasible box −1 ≤ x ≤ 1, q = 2n):

| n | iterations | Santini | solveQP | solveQPFast | solveQP vs Santini | |---:|---:|---:|---:|---:|---:| | 100 | 54 | 17.6 | 6.4 | 6.1 | 2.8× faster | | 200 | 112 | 179.5 | 61.2 | 58.7 | 2.9× faster | | 400 | 241 | 2082 | 776 | 746 | 2.7× faster |

Here solveQPFast ≈ solveQP: it parallelizes only the factorization, which is a minor part of a constraint‑heavy solve — the (sequential) iteration loop is the bottleneck. solveQPFast earns its keep in regime (b), below.

(b) Factorization‑dominated problems — few active constraints ⇒ ≈ 1 iteration ⇒ the time is almost entirely the Cholesky and triangular inverse, which is the same work in both, so the two are on par (node benchmark/benchmark.js, n box constraints):

| n | Santini | solveQP | |---:|---:|---:| | 200 | 5.5 | 5.5 | | 400 | 50 | 49 | | 800 | 507 | 412 | | 1024 | 1178 | 1116 |

Regime (b) is exactly where solveQPFast helps: it runs that factorization on the worker pool. Same box problems, solveQPFast vs solveQP (identical results; below n = 512 it just calls solveQP, hence ≈ 1× there):

| n | Node | Bun | |---:|---:|---:| | 600 | 1.4× | 2.9× | | 800 | 1.6× | 5.2× | | 1024 | 2.2× | 4.8× |

Run node benchmark/compare.js for the full Node‑vs‑Bun, three‑solver table.

Dependencies

None. solveQPFast uses only the built‑in worker_threads and SharedArrayBuffer, available in both Node and Bun.

License

MIT. Provided free of charge; the author is not liable for any damages arising from its use.

Contact

Questions, bugs, suggestions: [email protected].