1"""Small primal-dual active-set solver for strictly convex inequality QPs.
 2
 3Solves min .5 v'Hv + c'v subject to A v >= b using the KKT sign convention
 4H v + c - A' mu = 0, mu >= 0.
 5"""
 6import numpy as np
 7
 8
 9def solve_active_set(H, c, A, b, active_init=None, tol=1e-10, max_iter=100):
10    H = np.asarray(H, float); c = np.asarray(c, float)
11    A = np.asarray(A, float); b = np.asarray(b, float)
12    n, m = H.shape[0], A.shape[0]
13    active = set([] if active_init is None else active_init)
14    history = []
15    v = np.linalg.solve(H, -c)
16    mu = np.zeros(m)
17    for it in range(1, max_iter + 1):
18        inds = sorted(active)
19        if inds:
20            K = np.block([[H, -A[inds].T], [A[inds], np.zeros((len(inds), len(inds)))]])
21            rhs = np.r_[-c, b[inds]]
22            try:
23                sol = np.linalg.solve(K, rhs)
24            except np.linalg.LinAlgError:
25                # A dependent active set is not useful; drop its last member.
26                active.remove(inds[-1]); history.append(("remove_singular", inds[-1])); continue
27            v = sol[:n]; mu[:] = 0.; mu[inds] = sol[n:]
28        else:
29            v = np.linalg.solve(H, -c); mu[:] = 0.
30        residual = A @ v - b
31        bad_inactive = [(residual[j], j) for j in range(m) if j not in active and residual[j] < -tol]
32        bad_active = [(mu[j], j) for j in active if mu[j] < -tol]
33        history.append(("solve", tuple(sorted(active)), float(residual.min()) if m else 0.,
34                        float(mu.min()) if m else 0.))
35        if not bad_inactive and not bad_active:
36            return v, mu, sorted(active), it, history
37        if bad_active:
38            # Most negative multiplier must leave the active set.
39            _, j = min(bad_active)
40            active.remove(j); history.append(("remove", j))
41        else:
42            # Most violated inactive constraint enters.
43            _, j = min(bad_inactive)
44            active.add(j); history.append(("add", j))
45    raise RuntimeError("PDAS did not converge")
46
47
48def kkt_residual(H, c, A, b, v, mu):
49    stationarity = H @ v + c - A.T @ mu
50    primal = np.maximum(0., b - A @ v)
51    dual = np.maximum(0., -mu)
52    complementarity = np.abs(mu * (A @ v - b))
53    return max(np.max(np.abs(stationarity)), np.max(primal), np.max(dual), np.max(complementarity))