API reference

Modeling

moment_sos.variables(n)[source]

Return a list of n coordinate polynomials x_0 … x_{n-1}.

class moment_sos.Poly(n, terms=None)[source]

A polynomial in n variables: {exponent tuple: coefficient}.

static constant(n, c)[source]

Constant polynomial c in n variables.

static variable(n, i)[source]

Coordinate polynomial x_i in n variables.

degree()[source]

Total degree (0 for the zero polynomial).

substitute_affine(c, d)[source]

Return p(c + d * t) as a polynomial in t.

class moment_sos.Constraint(poly, kind)[source]

poly >= 0 (kind=’ineq’) or poly == 0 (kind=’eq’).

class moment_sos.Problem(objective, constraints=(), sense='max', box=None, boolean=False, add_box_constraints=True, name='')[source]

max (or min) f(x) s.t. g_i(x) >= 0, h_j(x) = 0, x in a box (or x in {-1,+1}^n).

Parameters:

  • objective (Poly)

  • constraints (iterable of Constraint (build them with >=, <=, == on Poly))

  • sense ("max" or "min")

  • box ((lo, hi) scalars or length-n sequences. Required unless boolean.) – The feasible set MUST lie in the box (the certified bound relies on it).

  • boolean (True (all variables) or a length-n boolean sequence marking the variables that are) – spins x_i in {-1,+1}; x_i^2 = 1 is then used as a rewriting rule (mixed problems allowed).

  • add_box_constraints (add (hi-x)(x-lo) >= 0 for every variable (recommended).)

to_original(t)[source]

Map normalised coordinates t in [-1,1]^n back to the original variables.

to_internal(x)[source]

Map original coordinates x to the normalised box [-1,1]^n used internally.

value(x)[source]

Original objective at an original-coordinate point.

violation(x)[source]

Largest violation of the original constraints / box at x.

Solving

moment_sos.solve(problem, order=2, cliques='dense', etas=(0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.0001, 3e-05, 1e-05), iters_per_stage=500, history_size=30, time_limit=None, device='cpu', theta0=None, verbose=False, warm_start=True, warm_iters=None, stop_at=None, tight_weights=False, term_sparsity=None, extra_basis=None)[source]

Certified bound through the order-order moment-SOS relaxation (entropic smoothing + L-BFGS).

Parameters:

  • problem (Problem)

  • order (int) – Relaxation order; must satisfy 2 * order >= degree of every polynomial (see min_order()).

  • cliques ("dense", "auto" or list of index lists) – One clique, the correlative-sparsity cliques, or explicit cliques.

  • etas (sequence of float) – Smoothing parameters of the continuation (SHORT_ETAS starts coarser).

  • iters_per_stage (int) – L-BFGS iterations per smoothing value.

  • time_limit (float, optional) – Seconds; the best certified bound so far is returned.

  • warm_start (bool) – Solve order 1, 2, … first and embed; the start is already a valid certificate equal to the lower-order bound.

  • stop_at (float, optional) – Stop as soon as the certified bound is at least this good (<= stop_at for a max problem, >= for a min problem), e.g. best_cut + 1 - 1e-9 to certify an integer optimum.

  • term_sparsity (int or "max", optional) – Split each clique’s moment block into overlapping sub-blocks by term-sparsity closure (steps iterations, or until stable). The bound stays valid for any block structure but the relaxation is weaker than the dense one; warm start is disabled. See the sparsity guide for when it helps.

  • device (str) – "cpu"; "cuda" is accepted but untested.

  • history_size (int)

  • theta0 (Tensor | None)

  • verbose (bool)

  • warm_iters (int | None)

  • tight_weights (bool)

Returns:

Result.bound is a valid bound (up to floating-point eigensolver error) at the best iterate seen; every iterate is a valid certificate, so a time_limit simply stops early.

Return type:

Result

moment_sos.solve_alm(problem, order=2, cliques='dense', rho0=1000.0, outer=300, inner=600, time_limit=None, device='cpu', theta0=None, warm_start=False, verbose=False, tol=1e-07, adapt_rho=False, shift='uniform', tight_weights=False, inner_solver='lbfgs', term_sparsity=None, extra_basis=None)[source]

Dual augmented-Lagrangian (boundary-point type) solver for min sigma + sum_b w_b t_b s.t. Q_b(theta) + t_b I >= 0.

  • inner problems: smooth (C^1), bound-constrained, solved by L-BFGS-B (inner_solver=”lbfgs”) or by a semismooth Newton-CG method with exact spectral Hessian-vector products (inner_solver=”newton”);

  • penalty parameter doubled whenever the primal infeasibility does not drop by 30% (residual balancing, adapt_rho, is available but misbehaves with inexact inner solves);

  • warm start: order r is started from the full ALM state (theta, t, X) of order r-1, embedded by zero padding;

  • every outer iterate gives a certified bound U(theta); the best one is returned.

shift“uniform” (default) – one scalar shift t_b per block, G_b = Q_b + t_b I, value inf U = rho_r^tr;
“diag” – one vector d_b in R^{s_b} per block (Proposition prop:diag), G_b = Q_b + diag(d_b),

objective sigma + sum_b (w_b/s_b) 1^T d_b with w_b/s_b = 1 (moment) or B_i (loc), d_b >= 0 for ‘plus’ blocks and free for constant-trace (Boolean) blocks; value inf U^d = rho_r^d (diagonal cuts). The certified bound at an iterate is U^d(theta, d) = sigma + sum_b (w_b/s_b) 1^T d_b + sum_b w_b Phi_b(Q_b + diag d_b) (the uniform-shift certificate applied to the shifted blocks; valid by Theorem thm:certificate), and Result.U is the best such bound over the outer iterates. Result.extra[“d”] holds the final shift vector (concatenated over blocks, offsets Result.extra[“d_off”]) and Result.extra[“U_uniform”] the best uniform certificate U(theta) seen along the way.

Parameters:
Return type:

Result

moment_sos.min_order(problem)[source]

Smallest relaxation order that can represent every polynomial of the problem.

class moment_sos.Result(problem: 'Problem', relaxation: 'Relaxation', oracle: 'Oracle', theta: 'torch.Tensor', U: 'float', history: 'list' = <factory>, time: 'float' = 0.0, n_eval: 'int' = 0, extra: 'dict' = <factory>)[source]
Parameters:
property bound

an upper bound if sense=’max’, a lower bound if sense=’min’.

Type:

Certified bound on the optimum of the ORIGINAL problem

round(**kw)[source]

Heuristic primal recovery: best feasible point found from the dual iterate (see moment_sos.rounding).

rigorous()[source]

Bound that holds in exact arithmetic for the stored normalised problem (see moment_sos.verify).

summary()[source]

One-line description of the bound and of the relaxation.

moment_sos.DEFAULT_ETAS = (0.1, 0.03, 0.01, 0.003, 0.001, 0.0003, 0.0001, 3e-05, 1e-05)

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

moment_sos.SHORT_ETAS = (0.001, 0.0003, 0.0001, 3e-05, 1e-05)

Built-in immutable sequence.

If no argument is given, the constructor returns an empty tuple. If iterable is specified the tuple is initialized from iterable’s items.

If the argument is a tuple, the return value is the same object.

Adaptive lifting (triangle walk)

moment_sos.walk(problem, rounds=10, per_round=None, rule='nms', mode='cuts', max_candidates=4000, drop_inactive=True, seed=0, time_limit=None, verbose=False, **solver_kw)[source]

Walk from order 1 towards order 2 by adding triples chosen by rule; every round gives a certified bound.

Parameters:

  • problem (Problem) – Fully spin (boolean=True) problem with an objective of degree at most 2 (MAX-CUT, Ising, QUBO in spins).

  • rounds (int) – Number of rounds after the initial order-1 solve.

  • per_round (int, optional) – Triples added per round K (default: the number of variables n).

  • rule ({"nms", "topk", "random"})

  • mode ({"cuts", "basis"}) – Impose the triples as triangle inequalities (order-1 relaxation with scalar localizing blocks) or as pair monomials in the moment basis (partial order-2 lifting).

  • drop_inactive (bool) – (mode="cuts") discard triples whose four inequalities have zero multiplier after a round.

  • time_limit (float, optional) – Total seconds; the walk stops after the round that exceeds it.

  • **solver_kw – Passed to moment_sos.solve_alm() (the primal moment matrix needed for the separation is read from its multipliers).

Return type:

WalkResult

class moment_sos.WalkResult(bounds=<factory>, times=<factory>, sizes=<factory>, triples=None, result=None, rule='', mode='', stopped='')[source]

Outcome of walk(). bounds[r] is the certified bound after round r (r = 0 is order 1).

Parameters:
moment_sos.select_triples(Y, K, rule='nms', max_candidates=4000, rng=None)[source]

Pick up to K triples of vertices for the next round (see the module docstring for the rules).

Returns an (m, 3) int array, m <= K; empty when no triangle inequality is violated (rule "random" ignores violations and always returns K distinct triples).

moment_sos.triangle_violations(Y, max_candidates=4000, chunk=8, tol=1e-06)[source]

Exact separation of the triangle inequalities over all triples i < j < k.

Parameters:

  • Y ((n, n) array) – Moment matrix entries Y_ij = E[x_i x_j] (diagonal ignored).

  • max_candidates (int) – Keep this many most violated triples.

  • chunk (int) – Number of i rows processed at once (memory/speed trade-off).

  • tol (float) – Only triples violated by more than tol are returned.

Returns:

tris – The m <= max_candidates most violated triples, sorted by decreasing violation max_s [-1 - s . (Y_ij, Y_ik, Y_jk)] over the four valid sign patterns.

Return type:

(m, 3) int array, viol : (m,) float array

Relaxation internals

class moment_sos.Relaxation(prob, order, cliques='dense', tight_weights=False, term_sparsity=None, extra_basis=None)[source]

Pure-numpy index structure of the order-order relaxation of a Problem.

Parameters:
moment_sos.correlative_sparsity_cliques(prob)[source]

Maximal cliques of a chordal extension (min-fill-in) of the correlative sparsity graph.

Parameters:

prob (Problem)

class moment_sos.Oracle(relax, device='cpu', dtype=torch.float64)[source]
Parameters:

relax (Relaxation)

moment_sos.lift_theta(res, relax_high, oracle_high)[source]

Embed a dual point of order r into order r+1 (zero-padded Gram matrices). The embedded point is a dual point of the larger relaxation; its certified bound equals the old one when no diagonal shift was active and is slightly larger otherwise (the shift weight grows with the block size).

Parameters:
Return type:

Tensor

Verification and references

moment_sos.rigorous_bound(result)[source]

Rigorous (exact-arithmetic) upper bound U_rig for the normalised problem and the corresponding bound for the user’s problem (upper bound for max, lower bound for min). Returns a dict.

moment_sos.verified_lmin_lower_bound(Q, margin=None, max_doublings=40)[source]

ell <= lambda_min(Q) for a symmetric float64 matrix Q, or -inf if the factorisation never succeeds.

moment_sos.solve_dual_ipm(R, solver=None, verbose=False)[source]

Solve the exact dual min sigma + sum_b w_b t_b s.t. Q_b(theta) + t_b I >= 0 with an interior-point solver (cvxpy). Returns (value in the internal normalisation, seconds).

Parameters:

R (Relaxation)

Serialization

moment_sos.problem_to_dict(problem)[source]

Serialise a Problem (the user’s original, un-normalised data).

moment_sos.problem_from_dict(d)[source]

Build a Problem from the dict format described in the module docstring.

moment_sos.load_problem(path)[source]
moment_sos.save_problem(problem, path)[source]
moment_sos.result_to_dict(res, rigorous=False, round_point=True)[source]

Machine-readable summary of a Result (bound, kind, point, history).

Neural-network encodings

Polynomial-optimization encodings of ReLU-network questions.

  • relu_verification – max of a linear function of the logits over an l_inf ball around x0. Exact ReLU constraints z>=0, z>=pre, z*(z-pre)=0 (+ triangle) for unstable neurons only; stable neurons are substituted away. Order 1 of the relaxation is the SDP of Raghunathan et al. (a Shor relaxation); orders >= 2 are strictly stronger.

  • lipschitz_problem – squared l2 Lipschitz bound max over |u|<=1, s in {0,1}^w of |W1 diag(s) W0 u|^2 of a one-hidden-layer ReLU network (the activation-pattern bound that LipSDP relaxes).

moment_sos.nn.ibp(Ws, bs, lo, hi)[source]

Interval bound propagation: list of pre-activation (lo, hi) per layer (last layer = logits).

moment_sos.nn.forward(Ws, bs, x)[source]

Forward pass of a ReLU network given weight matrices Ws and biases bs (last layer linear).

moment_sos.nn.relu_verification(Ws, bs, x0, eps, c, input_lo=None, input_hi=None, triangle=True, name='verify')[source]

Return (problem, info). c is the vector whose product with the logits is maximised.

moment_sos.nn.verification_exact_milp(Ws, bs, x0, eps, c, input_lo=None, input_hi=None, time_limit=120)[source]

Exact max of c^T logits over the l_inf ball by big-M MILP (scipy/HiGHS). Returns (value, status).

moment_sos.nn.verification_pgd(Ws, bs, x0, eps, c, input_lo=None, input_hi=None, steps=300, restarts=20, seed=0)[source]

Projected-gradient lower bound on the same maximisation (feasible points).

moment_sos.nn.lipschitz_problem(W0, W1, name='lipschitz', mixed=True)[source]

Squared l2 Lipschitz bound as a POP: max |W1 diag(s) W0 u|^2 s.t. |u|<=1, s in {0,1}^w (s = (1+t)/2, t spin). Variables [u, t].

moment_sos.nn.lipschitz_exact_pattern(W0, W1)[source]

max over activation patterns of ||W1 diag(s) W0||_2^2 by enumeration (w <= ~20).

moment_sos.nn.lipsdp_neuron(W0, W1)[source]

LipSDP-Neuron bound (Fazlyab et al. 2019) on the squared Lipschitz constant, via cvxpy (MOSEK if installed, else Clarabel).