Constraint Layers #

Raises:

ValueError – Concrete layers raise when lmbda is None.

abstractmethod hessian_curvature(w, lmbda)[source]#

Return the weighted constraint curvature ∑ λ_i ∇^2 f_i as a MatrixSystem.

This is the second-derivative term alone, without the Jacobian outer-product term that hessian_contribution adds. It is the constraint side of the Lagrangian Hessian

∇^2 L = ∇^2 f_0 + ∑ λ_i ∇^2 f_i,

which the barrier method minimizes over w when it evaluates the dual.

Abstract with no default, so every layer states its own curvature: affine layers (such as BoundsConstraintLayer, ImbalanceConstraintLayer and LinearInequalityConstraintLayer) return a ZeroSystem, and curved layers (such as NormConstraintLayer and SquaredNormConstraintLayer) return their real curvature. Requiring the declaration keeps a curved layer from silently inheriting zero curvature and mis-assembling the dual Hessian.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M. The curvature of a quadratic layer does not depend on w, but the signature allows curvature that does.
lmbda (npt.NDArray[np.float64]) – Per-constraint weights λ, of length m.

Returns:

system – The weighted curvature ∑ λ_i ∇^2 f_i.

Return type:

barrier_objective(w, shift=0.0)[source]#

Evaluate the log barrier of the shifted constraints.

Computes

φ(w) = -∑ log(shift - f_i(w)).

Parameters:

w (npt.NDArray[np.float64]) – Point at which the barrier is evaluated, of length M.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

phi – The barrier value; inf when any shifted constraint is not strictly satisfied.

Return type:

barrier_gradient(w, shift=0.0)[source]#

Compute the gradient of the log barrier with respect to w.

With c_i = f_i(w) - shift, negative for strictly feasible w,

∇φ(w) = -G^T (1 / c),

where 1 / c is taken componentwise.

Parameters:

w (npt.NDArray[np.float64]) – Strictly feasible point at which the gradient is evaluated, of length M.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

grad_phi – The barrier gradient, of length M.

Return type:

npt.NDArray[np.float64]

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps the layer strictly feasible.

Finds the largest α such that every constraint satisfies

f_i(w + α dw) - (shift + α ds) < 0.

The pair (shift, ds) lets the relaxation move with the step: slack-form Phase I, which moves its slack s jointly with w, passes shift = s and ds equal to the slack’s search direction. The defaults shift = 0 and ds = 0 check the layer’s own constraints.

The default implementation backtracks from 1 / β by factors of β, evaluating constraints at each trial point. It handles any nonlinearity but is slow, and the largest value it can report is 1 / β: enough for the caller’s β α to reach the full unit step, but not the true maximum. Subclasses with closed-form bounds should override it, report the true maximum (inf when no constraint binds), and may ignore line_search.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Settings supplying β and the minimum step; None uses the defaults.

Returns:

alpha – The largest feasible trial step, at most 1 / β; a value no larger than the minimum step when no trial step is feasible.

Return type:

is_satisfied(w)[source]#

Test whether every constraint in this layer holds strictly at w.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when f_i(w) < 0 for every constraint.
Return type:: bool

class cvxium.problems.constraint_layers.SlackConstraintLayer[source]#

Bases: ConstraintLayer

Constraint layer whose Phase I relaxes the constraints with a slack variable.

When a slack layer is the active layer, Problem.make_phase1 builds a slack-form Phase I problem: it augments the decision variable with a slack s, relaxes the layer’s constraints to

f_i(w) ≤ s,

and minimizes s. The layer remains a constraint in that problem, never the objective, so it needs no objective face, and slack form works for any layer.

Every slack layer defined in Cvxium is affine, so is_linear is True here; a lifted layer (LiftedConstraintLayer) instead forwards the value of the layer it wraps.

is_linear = True#

Whether every constraint in this layer is affine in w, so that ∇^2 f_i = 0.

class cvxium.problems.constraint_layers.DirectConstraintLayer[source]#

Bases: ConstraintLayer, Objective

Single-constraint layer whose Phase I minimizes the constraint value directly.

When a direct layer is the active layer, Problem.make_phase1 builds a direct-form Phase I problem: the layer becomes the objective, and minimizing its value f(w), subject to the remaining (background) layers, drives it below zero without a slack variable.

A direct layer is therefore also an Objective, and this class supplies that face: num_constraints is fixed at 1, so the single-constraint requirement is structural; evaluate and gradient return the constraint value and gradient row; and at supplies the constraint’s curvature as the objective Hessian. Subclasses provide the constraint face (dimension, constraints, grad_constraints, hessian_contribution and hessian_curvature) and may offer barrier_parameter_hint and interior_point.

num_constraints(M)[source]#

Return one, the number of constraints in a direct layer.

Parameters:: M (int) – Dimension of the decision variable; the count does not depend on it.
Returns:: m – Always 1.
Return type:: int

evaluate(w)[source]#

Evaluate the objective face, the layer’s constraint value f(w).

Parameters:: w (npt.NDArray[np.float64]) – Point at which the constraint is evaluated, of length M.
Returns:: f – The unshifted constraint value.
Return type:: float

gradient(w)[source]#

Compute the objective-face gradient ∇f(w), the constraint’s gradient row.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the gradient is evaluated, of length M.
Returns:: grad_f – The constraint gradient, of length M.
Return type:: npt.NDArray[np.float64]

at(w)[source]#

Return the per-point context for using the constraint value as an objective.

When a direct layer is the objective of a direct-form Phase I problem, build_hessian_system needs the objective’s structured Hessian. That is the constraint’s own curvature ∇^2 f(w), with no Jacobian outer-product term (the layer acts as an objective, not as a barrier-wrapped constraint), which is hessian_curvature at unit weight. The inherited Objective.at would route the point’s Hessian to the constraint-face hessian_contribution, which raises without lmbda; this override supplies the curvature instead.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the objective is evaluated, of length M.
Returns:: point – The objective’s value, gradient and structured Hessian at w.
Return type:: ObjectivePoint

barrier_parameter_hint(x0, num_ineq, outer_tolerance)[source]#

Return a lower bound on the initial barrier parameter in direct-form Phase I.

When this layer is the objective of a direct-form Phase I problem, BarrierSolver raises its initial barrier parameter t0 to max(t0, hint) under every adaptive initialization strategy. Without such a floor, t0 can be too small relative to the objective’s magnitude: the centering then under-progresses on the objective, and the solver may spuriously certify infeasibility.

A subclass with a closed-form version of the heuristic

t0 ≈ m / (f(x0) - p_star),

where p_star is a lower bound on the optimal value (Boyd and Vandenberghe, 2004, section 11.3.1), can override this hook to return a problem-aware bound.

Parameters:

x0 (npt.NDArray[np.float64]) – Strictly feasible starting point of the Phase I problem.
num_ineq (int) – Number of inequality constraints in the Phase I problem (the background layers’ constraints), the m in the heuristic.
outer_tolerance (float) – Outer duality-gap tolerance, a floor for denominators that may be tiny when x0 is nearly optimal.

Returns:

hint – Lower bound on t0, or None (the default) when the layer has none to offer.

Return type:

float | None

References

Boyd, Stephen and Vandenberghe, Lieven, Convex Optimization, Cambridge University Press, 2004.

interior_point(equality=None)[source]#

Return a strictly interior point of this layer, or None if none is known.

A direct layer that knows a point w with f(w) < 0 in closed form can override this hook to supply it. When the layer is the only inequality of a feasibility problem, Solver._solve_feasibility takes that point as the answer, certifying feasibility with no Phase I solve. This also sidesteps the direct-form Phase I in that case, whose objective-only shape (the sole layer hoisted to the objective, leaving no inequality layer) the barrier solver cannot handle.

When the feasibility problem also carries equality constraints A w = b, they are passed as equality, so a layer that knows the equality-constrained solution in closed form (for example, the minimum-residual point of an ellipsoid restricted to A w = b) can return a point satisfying both. The caller rechecks the strict inequality and the equality residual, and falls back to the iterative Phase I if either fails. A layer that cannot satisfy the equality in closed form should return None when equality is given.

Parameters:: equality (EqualityConstraintLayer | None) – The problem’s equality constraints, or None for an inequality-only problem.
Returns:: w – A point with f(w) < 0 that also satisfies equality when one is given; None by default, since most layers know no such point.
Return type:: npt.NDArray[np.float64] | None

Linear Constraints#

Affine constraint layers for Cvxium’s interior-point solvers.

Collects the layers whose feasible regions are defined by affine relations on the decision variable w:

EqualityConstraintLayer: equality constraints A w = b. It holds the data, a cached SVD of A, and the factory that builds the per-step KKT saddle; EqualitySolver does the solving. It is not a ConstraintLayer: a Problem carries it separately.
BoundsConstraintLayer: box bounds lb ≤ w ≤ ub, with a diagonal barrier Hessian.
ImbalanceConstraintLayer: imbalance constraints |(B w - c)_i| ≤ ψ_i.
LinearInequalityConstraintLayer: two-sided linear inequalities cl ≤ A w ≤ cu, with a low-rank barrier Hessian.
BandedLinearInequalityLayer: one-sided inequalities A w ≤ cu whose barrier Hessian A^T diag(d) A is banded under the given column order, assembled as a BandedSystem and solved by banded Cholesky instead of a width-m Woodbury update.
ArrowLinearInequalityLayer: like the banded layer, but with a dense K-column border on the banded core, assembled as an ArrowSystem.

Every inequality layer here subclasses SlackConstraintLayer, so when one is the active layer, Problem.make_phase1 builds a slack-form Phase I problem for it.

class cvxium.problems.linear_constraints.EqualityConstraintLayer(A, b, saddle_factory=None)[source]#

Bases: object

Linear equality constraints A w = b, with a cached SVD of A.

Holds the equality data, the lazily computed and cached singular value decomposition of A that the equality-handling algorithms consume, and the factory that builds the per-step KKT saddle. It does no solving itself: the projection and KKT solves live in EqualitySolver, which operates on a layer. Separating the data and its factorization from the algorithm lets several problems share one layer, and thus one warm SVD, without recomputation.

Parameters:

A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – Equality constraint matrix, of shape (p, M), dense or sparse. A sparse A is useful only together with a saddle_factory that exploits it (such as SparseSaddleFactory); the default factory densifies it, and svd_A is dense-only.
b (npt.NDArray[np.float64]) – Equality right-hand side; must be one-dimensional of length p (checked by assert).
saddle_factory (SaddleFactory | None) – The helper that builds the per-step KKT saddle from a Hessian block and A. Injecting it is how the caller chooses which structure of A to exploit without this layer knowing or branching on it. None uses SchurComplementSaddleFactory (a dense Schur complement).

property num_eq_constraints#

Return the number p of equality constraints, the row count of A.

Returns:: p – Number of equality constraints.
Return type:: int

property dimension#

Return the decision-variable dimension M, the column count of A.

Returns:: M – The decision-variable dimension.
Return type:: int

property saddle_factory#

Return the factory that builds this layer’s per-step KKT saddle.

Problem passes it on when it builds the augmented equality layer of a slack-form Phase I problem, so that layer exploits the same structure.

Returns:: saddle_factory – The factory given at construction, or the default SchurComplementSaddleFactory.
Return type:: SaddleFactory

saddle_system(hessian)[source]#

Build the per-step KKT saddle for this equality over a Hessian block.

Pairs the equality data A with the supplied Hessian H into the bordered Newton system

     _        _
K = |  H   A^T |
    |_ A    0 _|

This parallels ConstraintLayer.hessian_contribution: the equality object builds the per-step solvable, just as a constraint layer builds its Hessian contribution. The injected saddle_factory builds the concrete saddle, so which structure of A is exploited is the caller’s choice, not this layer’s.

Parameters:: hessian (MatrixSystem) – The Hessian block H, of order M.
Returns:: saddle – The saddle-point system K, ready to solve.
Return type:: SaddleSystem

svd_A()[source]#

Compute and cache the thin SVD of A.

Delegates to aspect_ratio_svd in numerical_helpers, which takes the direct thin SVD, or, for a wide A (at least 1,000 columns and at most a tenth as many rows as columns), a QR factorization of A^T followed by the SVD of the small triangular factor. The result is cached on the layer, so every problem sharing the layer reuses it. EqualitySolver uses it for the equality projection and the barrier-parameter initialization.

Returns:

U (npt.NDArray[np.float64]) – Left singular vectors, of shape (p, k), where k = min(p, M).
s (npt.NDArray[np.float64]) – Singular values, of length k, in descending order.
Vh (npt.NDArray[np.float64]) – Right singular vectors, transposed, of shape (k, M), with A = U diag(s) Vh.

Raises:

NotImplementedError – When A is sparse; a sparse A must use the injected saddle path instead.

class cvxium.problems.linear_constraints.BoundsConstraintLayer(lb=None, ub=None, dimension=None)[source]#

Box-bound constraints lb ≤ w ≤ ub.

Either bound may be omitted (None), but not both. The constraint vector stacks the active sides, lower block first:

lower:   f_i(w) = lb_i - w_i ≤ 0,
upper:   f_i(w) = w_i - ub_i ≤ 0.

With both sides present and every bound finite, the vector has length 2 M, laid out as [lower block, upper block]; with one side present it has length M.

An infinite entry (-inf in lb, +inf in ub) means that component has no bound on that side and contributes no constraint: each block then holds only the components whose bound is finite, in index order. Keeping a row w_i - inf ≤ 0 would not be harmless: its barrier term -log(inf) makes the barrier objective -inf everywhere, which disables the line search, and its multiplier 0 times f_i = -inf makes the Lagrangian NaN.

Every Jacobian row is ±e_i, so the barrier Hessian contribution is diagonal and is returned as a DiagonalSystem.

Parameters:

lb (float | list[float] | npt.NDArray[np.float64] | None) – Lower bound on each component of w; a scalar broadcasts to all components, and None disables the lower bound. The constructor raises ValueError when both lb and ub are None, when they are not broadcast-compatible, or when some component has lb ≥ ub (a degenerate or infeasible box).
ub (float | list[float] | npt.NDArray[np.float64] | None) – Upper bound on each component of w; a scalar broadcasts to all components, and None disables the upper bound.
dimension (int | None) – Decision-variable dimension M. Required when neither bound is an array, since a scalar bound pins no size on its own. When an array bound pins M, an explicit dimension is only a cross-check. The constructor raises ValueError when no dimension can be resolved, when two array bounds disagree on length, or when dimension disagrees with an array bound.

property dimension#

Return the decision-variable dimension M resolved at construction.

Returns:: M – The decision-variable dimension.
Return type:: int

num_constraints(M)[source]#

Count the constraints: M per active side, less the infinite entries.

Parameters:: M (int) – Dimension of the decision variable.
Returns:: m – The number of finite bounds; 2 M when both sides are present and finite.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the stacked bound constraints, shifted by shift.

Returns the lower block followed by the upper block,

lb_i - w_i - shift,   w_i - ub_i - shift,

over the active sides and the components with a finite bound.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraints are evaluated, of length M.
shift (float) – Relaxation subtracted from every constraint value.

Returns:

c – The shifted constraint values, of length m.

Return type:

npt.NDArray[np.float64]

grad_constraints(w)[source]#

Compute the dense constraint Jacobian, rows -e_i (lower) and +e_i (upper).

Forms an M × M identity, so it costs O(M^2) memory; the solvers use the structured products and hessian_contribution instead.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M; only its length is used.
Returns:: G – The constraint Jacobian, of shape (m, M).
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the constraint Jacobian, stacking -v and v.

Returns -v_i over the finite lower bounds followed by v_i over the finite upper bounds, in O(M) time.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M; only its length is used.
v (npt.NDArray[np.float64]) – Vector of length M to multiply.

Returns:

Gv – The product G v, of length m.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed constraint Jacobian.

Splits u into its lower block u_lo and upper block u_up and returns

-u_lo + u_up,

with each block scattered to the components it covers, in O(M) time.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M; only its length is used.
u (npt.NDArray[np.float64]) – Vector of length m to multiply, in the layer’s constraint order.

Returns:

GTu – The product G^T u, of length M.

Return type:

npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the weighted diagonal barrier curvature as a DiagonalSystem.

Each bound is affine with gradient ±e_i, so its atom -(1 / f_i) ∇f_i ∇f_i^T is diagonal. With the slacks

s_lo = w - lb + shift,   s_up = ub - w + shift,

which are the negated shifted constraint values, coordinate i receives

η_i = g_lo,i + g_up,i,
g_lo,i = λ_lo,i / s_lo,i,   g_up,i = λ_up,i / s_up,i.

The barrier method’s λ = 1 / s recovers 1 / s_lo,i^2 + 1 / s_up,i^2; the primal-dual method passes its duals directly.

The system also carries η^{-1} for the stable assembly path. With one side present it is s / λ, an exact reciprocal. With both sides it factors out the dominant term,

1 / η_i = (1 / g_max) / (g_lo,i / g_max + g_up,i / g_max),
g_max = max(g_lo,i, g_up,i),

which stays finite even when g_lo,i + g_up,i overflows because a slack is tiny. With infinite bound entries, a component bounded on one side gets that side’s term alone, and one bounded on neither gets η_i = 0 and an infinite reciprocal.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M, strictly inside the shifted box.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, of length m, in the layer’s constraint order (lower block, then upper block). Required; the None default exists only so the signature matches Objective.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – The diagonal η, with its stably computed reciprocal.

Return type:

DiagonalSystem

Raises:

ValueError – When lmbda is None.

hessian_curvature(w, lmbda)[source]#

Return zero curvature, since every bound is affine (∇^2 f_i = 0).

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M; only its length is used.
lmbda (npt.NDArray[np.float64]) – Per-constraint weights λ, of length m; unused.

Returns:

system – A ZeroSystem of order M.

Return type:

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps the shifted bounds strictly satisfied.

Finds the largest α such that, on the active sides,

lb - shift - α ds < w + α dw < ub + shift + α ds.

The lower side can bind only where dw_i + ds < 0 and the upper side only where dw_i - ds > 0; the step is the smallest closed-form ratio over those components, so line_search is ignored. An infinite bound entry gives an infinite ratio and never binds.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly inside the shifted box, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Unused; accepted for signature compatibility.

Returns:

alpha – The largest feasible step; inf when no bound binds.

Return type:

is_satisfied(w)[source]#

Check whether lb < w < ub holds strictly on the active sides.

An infinite bound entry always holds.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when every finite bound holds strictly.
Return type:: bool

class cvxium.problems.linear_constraints.ImbalanceConstraintLayer(B, c, psi)[source]#

Imbalance constraints |(B w - c)_i| ≤ ψ_i, the infinity-norm bound ‖B w - c‖ ≤ ψ.

Encoded as 2 q linear inequalities: a first block bounding each residual from below, then a second bounding it from above,

f_i(w)     = -(B w - c)_i - ψ_i ≤ 0,   i = 1, ..., q,
f_{q+i}(w) =  (B w - c)_i - ψ_i ≤ 0,   i = 1, ..., q.

The Jacobian is [-B; B], and the barrier Hessian contribution is a width-2q low-rank update, suited to q much smaller than M.

Parameters:

B (npt.NDArray[np.float64]) – Constraint matrix, of shape (q, M).
c (npt.NDArray[np.float64]) – Constraint target, of length q.
psi (float | list[float] | npt.NDArray[np.float64]) – Tolerance ψ on each row, of length q; a scalar broadcasts to all rows. Must be positive for the constraints to have a strict interior.

property q#

Return the number q of imbalance rows, the row count of B.

Returns:: q – Number of imbalance rows; the layer has 2 q constraints.
Return type:: int

property dimension#

Return the decision-variable dimension M, the column count of B.

Returns:: M – The decision-variable dimension.
Return type:: int

property psi#

Return the per-row tolerance ψ as an array.

Returns:: psi – The tolerance on each row, of length q.
Return type:: npt.NDArray[np.float64]

num_constraints(M)[source]#

Count the constraints, 2 q (one lower and one upper per row).

Parameters:: M (int) – Dimension of the decision variable; the count does not depend on it.
Returns:: m – Always 2 q.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the stacked imbalance constraints, shifted by shift.

Returns -(B w - c) - ψ - shift followed by (B w - c) - ψ - shift.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraints are evaluated, of length M.
shift (float) – Relaxation subtracted from every constraint value.

Returns:

c – The shifted constraint values, of length 2 q.

Return type:

npt.NDArray[np.float64]

grad_constraints(w)[source]#

Compute the constraint Jacobian G = [-B; B].

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused, since the layer is affine.
Returns:: G – The constraint Jacobian, of shape (2 q, M).
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the constraint Jacobian [-B; B].

Forms B v once and returns [-B v, B v].

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused.
v (npt.NDArray[np.float64]) – Vector of length M to multiply.

Returns:

Gv – The product G v, of length 2 q.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed constraint Jacobian.

Computes

[-B; B]^T u = B^T (u_up - u_lo),

where u_lo holds the first q entries of u and u_up the last q, with one product by B^T.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused.
u (npt.NDArray[np.float64]) – Vector of length 2 q to multiply.

Returns:

GTu – The product G^T u, of length M.

Return type:

npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the imbalance rows’ barrier curvature as a low-rank update.

The layer stands in for a handful of two-sided rows, far fewer than the dimension, so its weighted curvature

∑ (-λ_i / f_i) ∇f_i ∇f_i^T

is a genuine width-2q low-rank update on a zero base (see _low_rank_barrier_contribution).

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, of length 2 q, in the layer’s constraint order. Required; the None default mirrors ConstraintLayer.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – A ZeroSystem carrying the single low-rank update.

Return type:

Raises:

ValueError – When lmbda is None. Propagated from _low_rank_barrier_contribution.

hessian_curvature(w, lmbda)[source]#

Return zero curvature, since the imbalance rows are affine (∇^2 f_i = 0).

The low-rank structure in hessian_contribution is the barrier Jacobian outer-product term -(1 / f_i) ∇f_i ∇f_i^T, not constraint curvature.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M; only its length is used.
lmbda (npt.NDArray[np.float64]) – Per-constraint weights λ, of length 2 q; unused.

Returns:

system – A ZeroSystem of order M.

Return type:

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps every imbalance row strictly satisfied.

Finds the largest α such that, componentwise,

|B (w + α dw) - c| < ψ + shift + α ds.

With r = B w - c, the upper side requires

α (B dw - ds) < ψ + shift - r,

and the lower side requires

α (-B dw - ds) < ψ + shift + r.

Each side binds only on rows where the coefficient of α is positive; the step is the smallest closed-form ratio over those rows, so line_search is ignored.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible for the shifted constraints, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Unused; accepted for signature compatibility.

Returns:

alpha – The largest feasible step; inf when no row binds.

Return type:

is_satisfied(w)[source]#

Check whether |(B w - c)_i| < ψ_i holds strictly on every row.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when every row holds strictly.
Return type:: bool

class cvxium.problems.linear_constraints.LinearInequalityConstraintLayer(A, cl=None, cu=None)[source]#

Two-sided linear inequalities cl ≤ A w ≤ cu.

Either side of a row may be omitted by passing ±inf in the corresponding entry of cl or cu. The constraint vector stacks the active sides, upper rows first:

upper rows (finite cu_i):   f(w) = (A w)_i - cu_i ≤ 0,
lower rows (finite cl_i):   f(w) = cl_i - (A w)_i ≤ 0.

A genuinely two-sided row contributes one constraint to each block, so the layer has num_cu + num_cl constraints. The barrier Hessian contribution is a low-rank update of that width, suited to few rows relative to M.

Parameters:

A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – Constraint matrix, of shape (m, M), dense or sparse. Must be two-dimensional, or the constructor raises ValueError.
cl (npt.NDArray[np.float64] | None) – Lower bound on each row, of length m; -inf disables a row’s lower side, and None disables every lower side. The constructor raises ValueError when cl or cu does not have length m, when some row has both sides infinite (it would contribute no constraint, so passing None for both is rejected), or when some row has cl_i ≥ cu_i (an empty strict interior).
cu (npt.NDArray[np.float64] | None) – Upper bound on each row, of length m; +inf disables a row’s upper side, and None disables every upper side.

Notes

A may be dense or any scipy sparse matrix. The matrix-vector products in constraints, grad_constraints_multiply, grad_constraints_transpose_multiply, btls_max_step and is_satisfied apply A with the @ operator, so a sparse A stays sparse there. grad_constraints densifies A, and hessian_contribution calls it (through _low_rank_barrier_contribution) to form the dense low-rank factor G^T, of shape (M, num_cu + num_cl), at every Newton step. A sparse A therefore saves work in the products but not in the Hessian assembly; BandedLinearInequalityLayer and ArrowLinearInequalityLayer assemble their Hessians without densifying A.

property dimension#

Return the decision-variable dimension M, the column count of A.

Returns:: M – The decision-variable dimension.
Return type:: int

property num_cu#

Return the number of active upper-bound rows (finite entries of cu).

Returns:: num_cu – Number of rows with a finite upper bound.
Return type:: int

property num_cl#

Return the number of active lower-bound rows (finite entries of cl).

Returns:: num_cl – Number of rows with a finite lower bound.
Return type:: int

num_constraints(M)[source]#

Count the active constraints across both sides.

Parameters:: M (int) – Dimension of the decision variable; must equal the column count of A.
Returns:: m – num_cu + num_cl.
Return type:: int
Raises:: ValueError – When M differs from the layer’s dimension.

constraints(w, shift=0.0)[source]#

Evaluate the stacked row constraints, shifted by shift.

Returns the upper rows followed by the lower rows,

(A w)_i - cu_i - shift,   cl_i - (A w)_i - shift,

over the rows with a finite bound on that side, from one product A w.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraints are evaluated, of length M.
shift (float) – Relaxation subtracted from every constraint value.

Returns:

c – The shifted constraint values, of length num_cu + num_cl.

Return type:

npt.NDArray[np.float64]

grad_constraints(w)[source]#

Return the dense constraint Jacobian G = [A_cu; -A_cl].

A_cu holds the rows of A with a finite upper bound and A_cl those with a finite lower bound. The layer is affine, so G is formed (and a sparse A densified) once, on the first call, and cached; hessian_contribution asks for it at every Newton step. The cached array is read-only, so a caller that needs to modify G must copy it, and the layer treats A as immutable after construction.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused, since the layer is affine.
Returns:: G – The constraint Jacobian, of shape (num_cu + num_cl, M); read-only.
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the constraint Jacobian, with one product A v.

Returns (A v) over the upper rows followed by -(A v) over the lower rows.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused.
v (npt.NDArray[np.float64]) – Vector of length M to multiply.

Returns:

Gv – The product G v, of length num_cu + num_cl.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed constraint Jacobian.

Scatters u to the row space (adding the upper-row entries and subtracting the lower-row entries, so a two-sided row receives their difference), then applies A^T once.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated; unused.
u (npt.NDArray[np.float64]) – Vector of length num_cu + num_cl to multiply, in the layer’s constraint order.

Returns:

GTu – The product G^T u, of length M.

Return type:

npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the linear-inequality barrier curvature as a low-rank update.

The layer assumes few rows relative to the dimension, so its weighted curvature

∑ (-λ_i / f_i) ∇f_i ∇f_i^T

is a genuine low-rank update, of width num_cu + num_cl, on a zero base (see _low_rank_barrier_contribution). The low-rank factor is the dense G^T, formed by grad_constraints even when A is sparse.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, in the layer’s constraint order. Required; the None default mirrors ConstraintLayer.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – A ZeroSystem carrying the single low-rank update.

Return type:

Raises:

ValueError – When lmbda is None. Propagated from _low_rank_barrier_contribution.

hessian_curvature(w, lmbda)[source]#

Return zero curvature, since the rows are affine (∇^2 f_i = 0).

The low-rank structure in hessian_contribution is the barrier Jacobian outer-product term, not constraint curvature.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M; only its length is used.
lmbda (npt.NDArray[np.float64]) – Per-constraint weights λ; unused.

Returns:

system – A ZeroSystem of order M.

Return type:

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps every active row strictly satisfied.

Finds the largest α such that, on the active rows,

cl - shift - α ds < A (w + α dw) < cu + shift + α ds.

The upper side binds only where (A dw)_i - ds > 0 and the lower side only where -(A dw)_i - ds > 0; the step is the smallest closed-form ratio over those rows, so line_search is ignored.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible for the shifted constraints, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Unused; accepted for signature compatibility.

Returns:

alpha – The largest feasible step; inf when no row binds.

Return type:

is_satisfied(w)[source]#

Check whether cl < A w < cu holds strictly on the active rows.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when every finite row bound holds strictly.
Return type:: bool

class cvxium.problems.linear_constraints.BandedLinearInequalityLayer(A, cu, *, ridge_scale=1e-08, ridge_floor=1e-12)[source]#

Bases: LinearInequalityConstraintLayer

One-sided linear inequalities A w ≤ cu whose barrier Hessian is banded.

A drop-in LinearInequalityConstraintLayer for problems whose columns are ordered so that the barrier Hessian

A^T diag(d) A,   d_r = -λ_r / f_r,

is banded: for example, difference inequalities w_a - w_b ≤ cu_r with the variables sorted so each row’s two nonzeros sit close together. The base class assembles that term as a width-m low-rank (Woodbury) update, whose m × m capacitance matrix costs O(m^3) to factor; this layer instead assembles a lower-banded BandedSystem, which banded Cholesky factors in O(M L^2), where L is the bandwidth.

Only the upper side is supported; pass difference-structured rows so the Hessian is genuinely banded. The bandwidth is found from A at construction, as the largest column span within a row. Constraint values, gradients and the step bound are inherited unchanged.

Parameters:

A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – Constraint matrix, of shape (m, M), dense or sparse, with columns ordered so that A^T diag(d) A is banded.
cu (npt.NDArray[np.float64]) – Upper bound on each row, of length m. Every entry must be finite, or the constructor raises ValueError (a row with an infinite bound and no lower side contributes no constraint).
ridge_scale (float) – Relative ridge, multiplied by the largest row weight d_r and added to the Hessian diagonal. A difference-constraint Hessian is a graph Laplacian, singular on the constant vector; when the rest of the problem adds only negligible diagonal (a generous box, say), banded Cholesky sees a matrix that is not positive definite and fails. The ridge keeps it strictly positive definite.
ridge_floor (float) – Absolute ridge added on top of the relative one.

Notes

At construction, _row_pairs enumerates the P column pairs of the rows and maps each to a fixed slot in the (L + 1) × M band storage. Each Newton step scales the pair products by the row weights and accumulates them with one np.bincount, costing O(P + L M), before the O(M L^2) banded factorization.

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Assemble A^T diag(-λ / f) A, plus a ridge, as a lower-banded system.

Scales the precomputed pair products by the row weights d_r = -λ_r / f_r and scatters them into lower band storage, then adds

ridge = ``ridge_scale`` max_r d_r + ``ridge_floor``

to the diagonal. Since the layer is affine, this is the whole contribution apart from the ridge.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, of length m. Required; the None default mirrors ConstraintLayer.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – A BandedSystem of order M and bandwidth L, in lower storage.

Return type:

Raises:

ValueError – When lmbda is None.

class cvxium.problems.linear_constraints.ArrowLinearInequalityLayer(A, cu, *, core_dimension, ridge_scale=1e-08, ridge_floor=1e-12)[source]#

Bases: LinearInequalityConstraintLayer

One-sided linear inequalities with a banded-plus-border barrier Hessian.

Like BandedLinearInequalityLayer, but the columns of A split in two: the first M (core_dimension) form a banded core, and the trailing K couple densely to it, for n = M + K columns in all. The barrier Hessian A^T diag(d) A then has the bordered (arrow) structure

     _         _
H = |  C     B  |
    |_ B^T   D _|

with C banded of order M, B a dense M × K border, and D a dense K × K corner. This arises when K shared variables enter every row of an otherwise difference-structured system. The Hessian is returned as an ArrowSystem, whose Schur-complement solve touches the core only through its banded factorization.

A small relative ridge is added to the core diagonal and the corner diagonal: a difference-constraint core is a graph Laplacian, singular on the constant vector, and the arrow’s K × K Schur complement is factored by Cholesky with no positive-definite fallback, so the ridge keeps both blocks numerically positive definite.

Parameters:

A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – Constraint matrix, of shape (m, n), dense or sparse. Its first core_dimension columns are the banded core and its trailing K columns the border.
cu (npt.NDArray[np.float64]) – Upper bound on each row, of length m. Every entry must be finite, or the constructor raises ValueError.
core_dimension (int) – Order M of the banded core. Must satisfy 0 < M < n, or the constructor raises ValueError.
ridge_scale (float) – Relative ridge, multiplied by the largest row weight d_r and added to the core and corner diagonals.
ridge_floor (float) – Absolute ridge added on top of the relative one.

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Assemble the barrier Hessian as a banded core, dense border, and corner.

Sorts the precomputed column pairs (a ≤ b) of the barrier Hessian A^T diag(d) A, with row weights

d_r = -λ_r / f_r,

by where they land: both columns in the core go to the lower-banded core, a core column paired with a border column goes to the M × K border, and two border columns go to the K × K corner (accumulated in its upper triangle, then symmetrized). Then adds

ridge = ``ridge_scale`` max_r d_r + ``ridge_floor``

to the core and corner diagonals.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length n.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, of length m. Required; the None default mirrors ConstraintLayer.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – An ArrowSystem of order n, with a BandedSystem core of order M.

Return type:

Raises:

ValueError – When lmbda is None.

Quadratic Constraints#

Quadratic and sublevel-set constraint layers for Cvxium’s interior-point solvers.

Defines three single-constraint layers, each a DirectConstraintLayer:

NormConstraintLayer: the Euclidean-norm bound ‖B w - c‖ ≤ φ, encoded in squared form.
SquaredNormConstraintLayer: the bound ‖w‖^2 ≤ φ, a cheaper special case with B = I and c = 0.
ObjectiveConstraintLayer: the sublevel set f(w) ≤ φ of any convex Objective.

Because each is a direct layer, it is also an Objective. When one is the active layer, Problem.make_phase1 builds a direct-form Phase I problem that minimizes the constraint value itself, subject to the remaining layers, and BarrierSolver solves it. Each layer supplies structured Hessians (a low-rank update, a diagonal plus a rank-one update, or the wrapped objective’s own structure) so that Newton steps exploit them.

class cvxium.problems.quadratic_constraints.NormConstraintLayer(B, c, phi)[source]#

Bases: DirectConstraintLayer

Euclidean-norm constraint ‖B w - c‖ ≤ φ, encoded in squared form.

The user-visible constraint is ‖B w - c‖ ≤ φ. Internally it is the single constraint

f(w) = ‖B w - c‖^2 - φ^2 ≤ 0,

which has the same feasible set (for φ > 0) but is smooth where B w = c. Its derivatives are

∇f = 2 B^T v,   ∇^2 f = 2 B^T B,   where v = B w - c,

so the Hessian is constant and of rank at most q. The squared form avoids the division by ‖v‖ that the unsquared form needs, and the numerical trouble it causes near B w = c.

As a DirectConstraintLayer, the layer serves as its own Phase I objective: a direct-form Phase I minimizes f(w), and its objective Hessian is hessian_curvature at unit weight, the rank-q update 2 B^T B on a zero base, which Newton steps exploit.

The shift and ds arguments of the layer protocol must be zero here, and constraints, hessian_contribution and btls_max_step assert it. A shift does not commute with squaring: the bound

‖B w - c‖^2 ≤ φ^2 + shift

is not equivalent to ‖B w - c‖ ≤ φ + shift, so accepting a non-zero shift would silently change the feasible region. No caller needs one: Problem.make_phase1 shifts only the active layer of a slack-form Phase I, and a direct layer always takes the direct form.

Parameters:

B (npt.NDArray[np.float64]) – Constraint matrix, of shape (q, M). Must be two-dimensional, or the constructor raises ValueError.
c (npt.NDArray[np.float64]) – Constraint target, of length q. Any other shape makes the constructor raise ValueError.
phi (float) – Norm bound φ. Must be positive, or the constructor raises ValueError.

property q#

Return the number of rows q of B.

Returns:: q – The row count of B.
Return type:: int

property dimension#

Return the decision-variable dimension M, the column count of B.

Returns:: M – The column count of B.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the squared-form constraint value ‖B w - c‖^2 - φ^2.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraint is evaluated, of length M.
shift (float) – Must be zero; see the class description.

Returns:

c – The constraint value, of length 1; negative when w strictly satisfies the bound.

Return type:

npt.NDArray[np.float64]

Raises:

AssertionError – When shift is not zero.

grad_constraints(w)[source]#

Compute the constraint Jacobian G = (2 B^T (B w - c))^T.

The squared form is smooth everywhere; G is zero where B w = c.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
Returns:: G – The constraint Jacobian, of shape (1, M).
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the constraint Jacobian.

Computes

G v = 2 (B w - c)^T (B v)

without forming G.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
v (npt.NDArray[np.float64]) – Vector of length M to multiply.

Returns:

Gv – The product G v, of length 1.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed constraint Jacobian.

Computes

G^T u = 2 u_0 B^T (B w - c)

without forming G. An empty u gives the zero vector.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
u (npt.NDArray[np.float64]) – Vector of length 1 (or empty) to multiply.

Returns:

GTu – The product G^T u, of length M.

Return type:

npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the weighted barrier curvature as a rank-(q + 1) update on zero.

The constraint f = ‖v‖^2 - φ^2, with v = B w - c, has

∇f = 2 B^T v,   ∇^2 f = 2 B^T B,

so its curvature atom is

C = 2 B^T B - (4 / f) (B^T v) (B^T v)^T.

Weighted by λ, this is a ZeroSystem carrying one low-rank update with factor

κ = [B^T | B^T v],   of shape (M, q + 1),

and weights

d = (2 λ, ..., 2 λ, -4 λ / f).

The barrier method passes λ = -1 / f, giving

(-2 / f) B^T B + (4 / f^2) (B^T v) (B^T v)^T,

while the primal-dual method passes its dual variable. Every weight is positive on the strict interior, where f < 0 and λ > 0, so the update has no downdate and no mixed inertia.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible, of length M.
lmbda (npt.NDArray[np.float64] | None) – The weight λ, of length 1. Required; the None default exists only so the signature matches Objective.hessian_contribution.
shift (float) – Must be zero; see the class description.

Returns:

system – The weighted curvature λ C.

Return type:

Raises:

ValueError – When lmbda is None, or when w is not strictly feasible (f ≥ 0).
AssertionError – When shift is not zero.

hessian_curvature(w, lmbda)[source]#

Return the weighted constraint curvature 2 λ B^T B as a rank-q update.

Since ∇^2 f = 2 B^T B, the curvature is a ZeroSystem carrying one low-rank update with factor B^T and every weight 2 λ, positive semidefinite for λ ≥ 0 and independent of w. At unit weight it is also the objective Hessian of a direct-form Phase I on this layer.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M; only its length is used.
lmbda (npt.NDArray[np.float64]) – The weight λ, of length 1.

Returns:

system – The weighted curvature 2 λ B^T B.

Return type:

interior_point(equality=None)[source]#

Return the minimum-residual point of the ellipsoid if strictly interior.

The constraint value f(w) is smallest where ‖B w - c‖ is smallest, so the deepest interior point, by constraint value, is the minimum-residual point on whatever set the feasibility problem imposes:

Without equality constraints, the least-squares solution w_star = B^+ c, computed from a thin SVD of B with a Tikhonov-damped pseudo-inverse (tikhonov_filter) so that a near-degenerate B stays numerically stable. It is the ellipsoid’s center when B is square and nonsingular.
With equality constraints A w = b, the equality-constrained least-squares solution, the w block of the KKT system shown below. It is solved through the equality layer’s saddle_system, so a sparse A stays sparse.

The KKT system is

[ 2 B^T B   A^T ] [ w ]   [ 2 B^T c ]
[    A       0  ] [ ν ] = [    b    ].

The point is returned only when it is strictly interior (‖B w_star - c‖ < φ) and, with equality constraints, satisfies |A w_star - b| < 1e-10 componentwise, the same gate Solver._solve_feasibility applies before accepting it. Otherwise, including when the saddle solve raises NewtonStepError (as it can for a rank-deficient B), the method returns None and the caller falls back to the iterative Phase I.

Parameters:: equality (EqualityConstraintLayer | None) – The problem’s equality constraints A w = b, or None for an inequality-only feasibility problem.
Returns:: w – The strictly interior minimum-residual point, satisfying A w = b when equality is given, or None when that point is not strictly interior or not found.
Return type:: npt.NDArray[np.float64] | None

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps ‖B (w + α dw) - c‖^2 below φ^2.

With v = B w - c, the constraint at step α is the quadratic inequality

a α^2 + 2 b α + e < 0,   a = ‖B dw‖^2,   b = v^T B dw,   e = ‖v‖^2 - φ^2.

Strict feasibility gives e < 0, so when a > 0 the roots have opposite signs and the positive root is the step to the boundary. When a = 0 (B dw = 0) the constraint value does not change along dw and the step is unbounded. The bound is closed form, so line_search is ignored.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Must be zero; see the class description.
ds (float) – Must be zero; see the class description.
line_search (LineSearchSettings | None) – Ignored.

Returns:

alpha – The step at which the constraint reaches the boundary; inf when B dw = 0.

Return type:

Raises:

AssertionError – When shift or ds is not zero.
ValueError – Propagated from math.sqrt when w is not strictly feasible and the discriminant is negative.

Notes

The positive root is

(-b + sqrt(b^2 - a e)) / a.

For b > 0 the numerator cancels catastrophically when |a e| is much smaller than b^2, so in that case the code uses the algebraically equal form

-e / (b + sqrt(b^2 - a e)).

is_satisfied(w)[source]#

Test whether ‖B w - c‖ < φ strictly.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when ‖B w - c‖^2 < φ^2.
Return type:: bool

class cvxium.problems.quadratic_constraints.SquaredNormConstraintLayer(phi, dimension)[source]#

Bases: DirectConstraintLayer

Squared-Euclidean-norm constraint ‖w‖^2 ≤ φ.

Encoded as the single constraint

f(w) = w^T w - φ ≤ 0.

It has the same feasible set as NormConstraintLayer(B=I, c=0, phi=sqrt(φ)). The dedicated class exists because with B = I and c = 0 the math collapses to the gradient 2 w and the Hessian 2 I, and the general class would spend matrix products on what are scalar operations.

As a DirectConstraintLayer, the layer serves as its own Phase I objective, with the diagonal Hessian 2 I, and it offers a barrier_parameter_hint for that problem. Unlike NormConstraintLayer, it accepts a shift, which relaxes the bound to w^T w ≤ φ + shift.

Parameters:

phi (float) – The bound φ on the squared norm. Must be positive, or the constructor raises ValueError.
dimension (int) – Decision-variable dimension M. Required because the scalar phi does not fix the length of w.

property dimension#

Return the decision-variable dimension M given at construction.

Returns:: M – The dimension passed to the constructor.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the shifted constraint value w^T w - φ - shift.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraint is evaluated, of length M.
shift (float) – Relaxation subtracted from the constraint value.

Returns:

c – The shifted constraint value, of length 1.

Return type:

npt.NDArray[np.float64]

grad_constraints(w)[source]#

Compute the constraint Jacobian G = (2 w)^T.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
Returns:: G – The constraint Jacobian, of shape (1, M).
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the constraint Jacobian.

Computes G v = 2 w^T v without forming G.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
v (npt.NDArray[np.float64]) – Vector of length M to multiply.

Returns:

Gv – The product G v, of length 1.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed constraint Jacobian.

Computes G^T u = 2 u_0 w without forming G. An empty u gives the zero vector.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
u (npt.NDArray[np.float64]) – Vector of length 1 (or empty) to multiply.

Returns:

GTu – The product G^T u, of length M.

Return type:

npt.NDArray[np.float64]

barrier_parameter_hint(x0, num_ineq, outer_tolerance)[source]#

Return a lower bound on t0 for a direct-form Phase I on this layer.

When this layer is the objective of a direct-form Phase I problem, BarrierSolver raises its initial barrier parameter t0 to max(t0, hint) under every adaptive initialization strategy (all but UNIT). Without that floor, t0 can be too small for the objective’s scale, and the centering under-progresses on the objective.

The hint is a closed-form lower bound on the heuristic

t0 ≈ m / (f(x0) - p_star)

(Boyd and Vandenberghe, 2004, section 11.3.1). The Phase I objective satisfies

f(w) = w^T w - φ ≥ -φ,

so p_star ≥ -φ and f(x0) - p_star ≤ x0^T x0, giving

hint = m / max(outer_tolerance, x0^T x0),

where m is num_ineq. The outer_tolerance floor avoids dividing by a tiny x0^T x0 when x0 is near the origin.

Parameters:

x0 (npt.NDArray[np.float64]) – Strictly feasible starting point of the Phase I problem.
num_ineq (int) – Number of inequality constraints in the Phase I problem (the background layers’ constraints), the m in the heuristic.
outer_tolerance (float) – Outer duality-gap tolerance, a floor for the denominator.

Returns:

hint – The lower bound m / max(outer_tolerance, x0^T x0); never None for this layer.

Return type:

float | None

References

Boyd, Stephen and Vandenberghe, Lieven, Convex Optimization, Cambridge University Press, 2004.

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the weighted barrier curvature as a diagonal plus a rank-one update.

With c = w^T w - φ - shift, negative on the strict interior, the constraint has ∇f = 2 w and ∇^2 f = 2 I, so its curvature atom is

C = 2 I - (4 / c) w w^T.

Weighted by λ, this is a DiagonalSystem with every entry 2 λ, carrying one rank-one update with factor w and weight -4 λ / c. With the barrier method’s weight

λ = -1 / c

this becomes

-(2 / c) I + (4 / c^2) w w^T,

while the primal-dual method passes its dual variable as λ. The diagonal’s reciprocal is passed explicitly as 1 / (2 λ), the same value DiagonalSystem computes by default; at the barrier weight it equals -c / 2.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible for the shifted bound, of length M.
lmbda (npt.NDArray[np.float64] | None) – The weight λ, of length 1, positive. Required; the None default exists only so the signature matches Objective.hessian_contribution.
shift (float) – Relaxation applied to the constraint value, as in constraints.

Returns:

system – The weighted curvature λ C.

Return type:

Raises:

ValueError – When lmbda is None, or when w is not strictly feasible for the shifted bound (c ≥ 0).
ZeroDivisionError – When λ = 0, from the reciprocal 1 / (2 λ).

hessian_curvature(w, lmbda)[source]#

Return the weighted constraint curvature 2 λ I as a diagonal.

Since ∇^2 f = 2 I, the curvature is a DiagonalSystem with every entry 2 λ, independent of w.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M; only its length is used.
lmbda (npt.NDArray[np.float64]) – The weight λ, of length 1.

Returns:

system – The weighted curvature 2 λ I.

Return type:

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest step that keeps the shifted squared-norm bound strict.

Finds the largest α with

‖w + α dw‖^2 < φ + shift + α ds,

that is,

a α^2 + 2 b α + e < 0,

with

a = dw^T dw,   b = w^T dw - ds / 2,   e = w^T w - φ - shift.

Strict feasibility gives e < 0. When a > 0 the roots have opposite signs, and the positive root

(-b + sqrt(b^2 - a e)) / a

is the bound. When a = 0 (dw = 0) the left side is e - α ds, which reaches zero at α = e / ds when ds < 0 and never otherwise. The bound is closed form, so line_search is ignored.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible for the shifted bound, of length M.
dw (npt.NDArray[np.float64]) – Step direction, of length M.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Ignored.

Returns:

alpha – The step at which the constraint reaches the boundary; inf when it never does.

Return type:

Raises:

ValueError – Propagated from math.sqrt when w is not strictly feasible and the discriminant is negative.

is_satisfied(w)[source]#

Test whether w^T w < φ strictly.

The test uses the unshifted bound.

Parameters:: w (npt.NDArray[np.float64]) – Point to test, of length M.
Returns:: satisfied – True when w^T w < φ.
Return type:: bool

class cvxium.problems.quadratic_constraints.ObjectiveConstraintLayer(objective, upper_bound)[source]#

Bases: DirectConstraintLayer

Sublevel-set constraint f(w) ≤ φ built from a convex Objective.

Wraps any convex Objective f as the single constraint

f(w) - φ ≤ 0,

reusing the objective’s value, gradient and structured Hessian. It mirrors the dual role of DirectConstraintLayer: a direct layer exposes a constraint value as an objective, and this layer exposes an objective as a constraint. A caller can therefore promote an objective to a constraint at the Problem level (as in an ε-constraint method) without a dedicated constraint class per objective form.

The constraint set is convex when f is convex; the caller is responsible for that. To encode a superlevel constraint g(w) ≥ φ on a concave g, wrap f = -g, which is convex, and pass -φ as the bound.

As a DirectConstraintLayer, the layer serves as its own Phase I objective, the scalar f(w) - φ, whose objective Hessian is ∇^2 f. The Jacobian products, btls_max_step and is_satisfied are the ConstraintLayer defaults; in particular, the step bound backtracks rather than solving in closed form.

Parameters:

objective (Objective) – The convex objective f whose sublevel set defines the constraint. Its dimension fixes M.
upper_bound (float) – The bound φ.

Notes

The barrier curvature atom at weight λ is

λ ∇^2 f - (λ / f_c) ∇f ∇f^T,

with f_c = f(w) - φ - shift < 0 on the strict interior. Since λ > 0 and f_c < 0, the rank-one term is an update, not a downdate.

property dimension#

Return the decision-variable dimension M of the wrapped objective.

Returns:: M – The wrapped objective’s dimension.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the shifted constraint value f(w) - φ - shift.

Parameters:

w (npt.NDArray[np.float64]) – Point at which the constraint is evaluated, of length M.
shift (float) – Relaxation subtracted from the constraint value.

Returns:

c – The shifted constraint value, of length 1.

Return type:

npt.NDArray[np.float64]

grad_constraints(w)[source]#

Compute the constraint Jacobian G = ∇f(w)^T.

Parameters:: w (npt.NDArray[np.float64]) – Point at which the Jacobian is evaluated, of length M.
Returns:: G – The wrapped objective’s gradient as a row, of shape (1, M).
Return type:: npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Return the weighted barrier curvature of the sublevel-set constraint.

Computes

λ ∇^2 f - (λ / f_c) ∇f ∇f^T,   f_c = f(w) - φ - shift,

as the wrapped objective’s hessian_contribution scaled by λ, carrying one rank-one update with factor ∇f(w) and weight -λ / f_c, so the structure of the objective’s Hessian is preserved. See the class Notes.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, strictly feasible for the shifted bound, of length M.
lmbda (npt.NDArray[np.float64] | None) – The weight λ, of length 1. Required; the None default exists only so the signature matches Objective.hessian_contribution.
shift (float) – Relaxation applied to the constraint value, as in constraints.

Returns:

system – The weighted barrier curvature.

Return type:

Raises:

ValueError – When lmbda is None, or when w is not strictly feasible for the shifted bound (f_c ≥ 0).

hessian_curvature(w, lmbda)[source]#

Return the weighted constraint curvature λ ∇^2 f(w).

This is the wrapped objective’s hessian_contribution at w scaled by λ, without the Jacobian outer-product term.

Parameters:

w (npt.NDArray[np.float64]) – Current iterate, of length M.
lmbda (npt.NDArray[np.float64]) – The weight λ, of length 1.

Returns:

system – The weighted curvature λ ∇^2 f(w).

Return type:

Lifted Constraints#

Dimension-lifting adapters that express slack-form Phase I as a Problem.

When the active layer of a Problem (its last inequality layer, the one whose feasibility Phase I seeks) is a SlackConstraintLayer, Problem.make_phase1 builds a slack-form Phase I problem: an ordinary Problem over the augmented variable

w_aug = [w; s]

of dimension M + 1, which minimizes the slack s. This module supplies the layers of that problem:

LiftedConstraintLayer lifts each background layer (already satisfied at the starting point) from dimension M to the augmented dimension; the layer sees only w.
_SlackActiveLayer wraps the active layer, relaxing its constraints to f_i(w) ≤ s.
_SlackUpperBoundLayer bounds the slack, s ≤ s_max, so the barrier subproblem has a unique solution.

LiftedConstraintLayer is public: a caller formulating its own problem over an augmented variable (an epigraph [w; t], say) lifts its existing layers with it, and its objective with LiftedObjective. The other two classes are construction details of Problem.make_phase1 and not part of the public API.

class cvxium.solvers.phase1.LiftedConstraintLayer(base, base_dimension, extra_dimension=1)[source]#

Constraint layer lifted to act on a larger decision variable.

The wrapped layer sees only the first base_dimension components of the augmented vector and ignores the trailing ones. The augmented Jacobian gets zero columns for the extra components, and the augmented Hessian contributions get a zero block.

Problem.make_phase1 uses it to embed the background layers, which act on w of dimension M, into a slack-form Phase I problem over w_aug = [w; s]. Callers that augment w with auxiliary variables of their own use it the same way, with extra_dimension set to the number of auxiliaries.

The class subclasses SlackConstraintLayer, so when a lifted layer is the last inequality layer of a caller-built problem, Problem.make_phase1 treats it like any other slack layer and builds a slack-form Phase I problem around it.

Parameters:

base (ConstraintLayer) – The layer being lifted. It sees w_aug[:base_dimension] only.
base_dimension (int) – Dimension M the base layer was constructed for.
extra_dimension (int) – Number of trailing components the base layer ignores. The default, 1, is the single slack of a slack-form Phase I. Must be at least 1, or the constructor raises ValueError.

Notes

is_linear is forwarded from the base layer, so a solver choosing its dual path does not mistake a lifted squared-norm layer for an affine one.

property dimension#

Return the augmented dimension, M plus the trailing components.

Returns:: n – base_dimension + extra_dimension.
Return type:: int

num_constraints(M)[source]#

Count the constraints, which are the base layer’s.

Parameters:: M (int) – Augmented dimension of the decision variable; ignored, since the base layer is queried at base_dimension.
Returns:: m – Number of constraints in the base layer.
Return type:: int

constraints(w, shift=0.0)[source]#

Evaluate the base layer’s shifted constraints on the leading components.

Parameters:

w (npt.NDArray[np.float64]) – Augmented point; only w[:base_dimension] is read.
shift (float) – Relaxation subtracted from every constraint value, as in ConstraintLayer.constraints.

Returns:

c – The base layer’s shifted constraint values, of length m.

Return type:

npt.NDArray[np.float64]

btls_max_step(w, dw, shift=0.0, ds=0.0, line_search=None)[source]#

Find the largest feasible step by delegating to the base layer.

The trailing components do not enter this layer’s constraints, so the augmented step is feasible exactly when its truncation is. Delegating, rather than inheriting the default backtracking search, lets a base layer with a closed-form bound keep it under lifting.

Parameters:

w (npt.NDArray[np.float64]) – Current augmented iterate, strictly feasible; only w[:base_dimension] is read.
dw (npt.NDArray[np.float64]) – Augmented step direction; only dw[:base_dimension] is read.
shift (float) – Relaxation at α = 0.
ds (float) – Rate at which the relaxation grows with α.
line_search (LineSearchSettings | None) – Settings passed through to the base layer; None uses the defaults.

Returns:

alpha – The base layer’s btls_max_step on the truncated point and direction.

Return type:

grad_constraints(w)[source]#

Compute the augmented Jacobian: base columns, then zero columns.

Parameters:: w (npt.NDArray[np.float64]) – Augmented point; only w[:base_dimension] is read.
Returns:: G – The base Jacobian padded with len(w) - base_dimension zero columns, of shape (m, len(w)).
Return type:: npt.NDArray[np.float64]

grad_constraints_multiply(w, v)[source]#

Multiply a vector by the augmented Jacobian.

The extra columns are zero, so

G v = G_base v[:base_dimension].

Parameters:

w (npt.NDArray[np.float64]) – Augmented point; only w[:base_dimension] is read.
v (npt.NDArray[np.float64]) – Augmented vector to multiply; only v[:base_dimension] is read.

Returns:

Gv – The product G v, of length m.

Return type:

npt.NDArray[np.float64]

grad_constraints_transpose_multiply(w, u)[source]#

Multiply a vector by the transposed augmented Jacobian.

Computes G_base^T u and pads it with zeros in the extra components.

Parameters:

w (npt.NDArray[np.float64]) – Augmented point; only w[:base_dimension] is read.
u (npt.NDArray[np.float64]) – Vector of length m to multiply.

Returns:

GTu – The product G^T u, of length len(w).

Return type:

npt.NDArray[np.float64]

hessian_contribution(w, lmbda=None, shift=0.0)[source]#

Lift the base layer’s Hessian contribution to the augmented dimension.

The base layer has curvature only in its own base_dimension coordinates, so its MatrixSystem becomes the first block of a BlockDiagonalSystem, with a ZeroSystem over the trailing coordinates. Problem.build_hessian_system sums this with the other layers’ full-dimension contributions, so the embedded base composes correctly whatever its structure (diagonal, low rank, and so on).

Parameters:

w (npt.NDArray[np.float64]) – Augmented iterate; only w[:base_dimension] is read.
lmbda (npt.NDArray[np.float64] | None) – Per-constraint weights λ, of length m, passed to the base layer. Required; the None default mirrors ConstraintLayer.hessian_contribution.
shift (float) – Relaxation applied to the constraint values, as in constraints.

Returns:

system – The block-diagonal lifted contribution, or the base contribution itself when len(w) equals base_dimension.

Return type:

Raises:

ValueError – Propagated from the base layer’s hessian_contribution when lmbda is None.

hessian_curvature(w, lmbda)[source]#

Lift the base layer’s weighted curvature to the augmented dimension.

The lifted constraint f_i(w_aug[:M]) has curvature only in the base coordinates, so the base hessian_curvature becomes the first block of a BlockDiagonalSystem, with a ZeroSystem over the trailing coordinates, as in hessian_contribution.

Parameters:

w (npt.NDArray[np.float64]) – Augmented iterate; only w[:base_dimension] is read.
lmbda (npt.NDArray[np.float64]) – Per-constraint weights λ, of length m.

Returns:

system – The block-diagonal lifted curvature, or the base curvature itself when len(w) equals base_dimension.

Return type: