Saddle-Point Problems#

Stationarity Problems#

Two-block stationarity problems: barrier-smoothed residual and coupled Newton step.

Cvxium’s single-block Problem minimizes one objective over one decision variable. This module states a different problem: find a joint stationary point of a smooth objective ζ(x, y) in two variable blocks, a large structured x-block and a small y-block, each confined to the strict interior of its own inequality domain. With ζ concave in x and convex in y, the stationary point is a saddle point,

V = max_x min_y ζ(x, y).

TwoBlockObjective is the abstract objective (value, block gradients, and curvature blocks), with QuadraticSaddleObjective as a concrete quadratic example. StationarityProblem bundles an objective with its block domains and an optional equality on the x-block, and assembles the residuals and coupled Newton systems for both the barrier method and the primal-dual method, without running any loop. TwoBlockStationaritySolver is the protocol the engines (StationaritySolver and PrimalDualStationaritySolver) satisfy, StationarityConvergenceError is what they raise on non-convergence, and two_block_result_message formats their result messages.

class cvxium.saddle.stationarity.TwoBlockObjective[source]#

Bases: ABC

Abstract smooth objective ζ(x, y) in two variable blocks.

The decision variables split into a large, typically structured x-block of dimension n and a small y-block of dimension m, ordered [x; y]. Subclasses supply the objective’s value, its two block gradients, and its three second-order blocks. The interface carries no minimax vocabulary: a joint stationary point of ζ is a point where both block gradients vanish, and whether it is a minimum, maximum or saddle depends on the signs of the curvature the subclass returns, not on this interface.

The one curvature sign the interface fixes is the x-block’s: neg_hessian_xx returns the negated x-curvature -∇^2_xx ζ, which a consumer assembles as a positive definite block. Everything else is returned raw (un-negated), and the coupled Newton system applies the remaining signs.

Notes

StationarityProblem builds the barrier-smoothed joint stationarity residual and the coupled Newton step from these pieces. Supplying the x-curvature already negated lets the large positive definite block flow straight into a CoupledSaddleSystem with no sign flip at the call site.

abstract property dimension_x#

Return the dimension n of the x-block, the large structured block.

Returns:: n – Dimension of the x-block.
Return type:: int

abstract property dimension_y#

Return the dimension m of the y-block, the small block.

Returns:: m – Dimension of the y-block.
Return type:: int

abstractmethod evaluate(x, y)[source]#

Evaluate the objective ζ(x, y).

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

value – The objective value ζ(x, y).

Return type:

float

abstractmethod gradient_x(x, y)[source]#

Compute the raw x-gradient ∇_x ζ(x, y).

Implementations return the plain partial gradient, un-negated; the consumer applies any sign it needs.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

grad – The gradient ∇_x ζ(x, y), of length n.

Return type:

npt.NDArray[np.float64]

abstractmethod gradient_y(x, y)[source]#

Compute the raw y-gradient ∇_y ζ(x, y).

Implementations return the plain partial gradient, un-negated; the consumer applies any sign it needs.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

grad – The gradient ∇_y ζ(x, y), of length m.

Return type:

npt.NDArray[np.float64]

abstractmethod neg_hessian_xx(x, y)[source]#

Return the negated x-curvature -∇^2_xx ζ(x, y) as a MatrixSystem.

The consumer relies on this block being positive definite: it is the large M_xx block of a CoupledSaddleSystem. A subclass must return -∇^2_xx ζ, not +∇^2_xx ζ.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

system – The negated x-curvature -∇^2_xx ζ(x, y), of order n.

Return type:

MatrixSystem

abstractmethod hessian_yy(x, y)[source]#

Return the raw y-curvature ∇^2_yy ζ(x, y) as a MatrixSystem.

Returned un-negated. The consumer relies on this block being positive definite: it is the small H_yy block of a CoupledSaddleSystem.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

system – The y-curvature ∇^2_yy ζ(x, y), of order m.

Return type:

MatrixSystem

abstractmethod mixed_hessian(x, y)[source]#

Return the mixed curvature ∇^2_xy ζ(x, y) as a dense matrix.

The matrix is tall and skinny (n ≫ m). It is dense, but cheap despite its n rows: the coupled system’s Schur step applies it as m batched right-hand-side columns, in one M_xx solve.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

mixed – The mixed curvature ∇^2_xy ζ(x, y), of shape (n, m).

Return type:

npt.NDArray[np.float64]

class cvxium.saddle.stationarity.QuadraticSaddleObjective(a, b, P, C, R)[source]#

Bases: TwoBlockObjective

Quadratic two-block objective, concave in x and convex in y.

Models

ζ(x, y) = a^T x + b^T y - (1/2) x^T P x + x^T C y + (1/2) y^T R y,

with P (n × n) symmetric positive definite, R (m × m) symmetric positive definite, and C (n × m) the coupling. The -(1/2) x^T P x term makes ζ concave in x and the +(1/2) y^T R y term makes it convex in y, so a joint stationary point is a saddle point. The block gradients and curvatures are

∇_x ζ = a - P x + C y,      ∇_y ζ = b + C^T x + R y,
-∇^2_xx ζ = P,              ∇^2_yy ζ = R,
∇^2_xy ζ = C.

Parameters:

a (npt.NDArray[np.float64]) – Linear x-coefficient, of length n. Must be one-dimensional, or the constructor raises ValueError.
b (npt.NDArray[np.float64]) – Linear y-coefficient, of length m. Must be one-dimensional, or the constructor raises ValueError.
P (npt.NDArray[np.float64]) – Symmetric positive definite x-curvature, of shape (n, n). ζ carries -P as its x-curvature; P itself is what neg_hessian_xx returns. A wrong shape raises ValueError.
C (npt.NDArray[np.float64]) – Coupling, of shape (n, m). A wrong shape raises ValueError.
R (npt.NDArray[np.float64]) – Symmetric positive definite y-curvature, of shape (m, m). A wrong shape raises ValueError.

Notes

P and R are symmetrized, as (P + P^T) / 2, so a slightly asymmetric input still gives a symmetric curvature. Definiteness is not checked. All arrays are copied at construction, so later changes to the caller’s arrays do not affect the objective.

property dimension_x#

Return the dimension n of the x-block, the length of a.

Returns:: n – Dimension of the x-block.
Return type:: int

property dimension_y#

Return the dimension m of the y-block, the length of b.

Returns:: m – Dimension of the y-block.
Return type:: int

evaluate(x, y)[source]#

Evaluate the quadratic objective ζ(x, y).

Computes

ζ(x, y) = a^T x + b^T y - (1/2) x^T P x + x^T C y + (1/2) y^T R y.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

value – The objective value.

Return type:

float

gradient_x(x, y)[source]#

Compute the raw x-gradient a - P x + C y.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

grad – The gradient ∇_x ζ, of length n.

Return type:

npt.NDArray[np.float64]

gradient_y(x, y)[source]#

Compute the raw y-gradient b + C^T x + R y.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.

Returns:

grad – The gradient ∇_y ζ, of length m.

Return type:

npt.NDArray[np.float64]

neg_hessian_xx(x, y)[source]#

Return the negated x-curvature -∇^2_xx ζ = P, which is constant.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n; unused, since the curvature is constant.
y (npt.NDArray[np.float64]) – The y-block, of length m; unused.

Returns:

system – A DenseSystem holding the symmetrized P.

Return type:

MatrixSystem

hessian_yy(x, y)[source]#

Return the y-curvature ∇^2_yy ζ = R, which is constant.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n; unused, since the curvature is constant.
y (npt.NDArray[np.float64]) – The y-block, of length m; unused.

Returns:

system – A DenseSystem holding the symmetrized R.

Return type:

MatrixSystem

mixed_hessian(x, y)[source]#

Return the mixed curvature ∇^2_xy ζ = C, which is constant.

Returns a copy, so a caller that mutates the result cannot corrupt the stored C and later evaluations.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n; unused, since the curvature is constant.
y (npt.NDArray[np.float64]) – The y-block, of length m; unused.

Returns:

mixed – A copy of C, of shape (n, m).

Return type:

npt.NDArray[np.float64]

class cvxium.saddle.stationarity.StationarityProblem(objective, x_domain=<factory>, y_domain=<factory>, x_equality=None)[source]#

Bases: object

Two-block stationarity problem with log-barrier domains.

Bundles a smooth TwoBlockObjective ζ(x, y) with the inequality ConstraintLayer lists whose strict interiors are the domains of the two blocks: x_domain constrains x, and y_domain constrains y. Each domain enters through its log barrier

φ(w) = -∑ log(-f_i(w)),

weighted by 1 / t, where the caller supplies t on each call (the outer loop that drives t up lives in the solvers).

Like Problem, this is a passive container: it states the barrier-smoothed joint stationarity residual and assembles the coupled Newton system, but runs no outer loop, line search, tolerance test or Phase I. The barrier-smoothed residual is

r = [r_x; r_y],

with

r_x = ∇_x ζ - (1/t) g_X(x),
r_y = ∇_y ζ + (1/t) g_Y(y),

where g_X, the sum of the x-domain layers’ barrier_gradient, is the gradient of the x-domain barrier, and g_Y likewise for the y-domain. The concave x-block subtracts its barrier gradient; the convex y-block adds it. The Jacobian of r is

 _              _
|  -M_xx   B     |
|_ B^T     H_yy _|

with

M_xx = -∇^2_xx ζ + (1/t) H_X,   H_yy = ∇^2_yy ζ + (1/t) H_Y,   B = ∇^2_xy ζ,

where H_X and H_Y are the barrier Hessians of the two domains, and M_xx and H_yy are positive definite. A Newton step solves this coupled saddle against (-r_x, -r_y).

The primal-dual counterparts (primal_dual_residual and build_primal_dual_coupled_system) replace the barrier terms with explicit inequality duals and reduce to the same saddle structure.

Fields:

objective (TwoBlockObjective) – The smooth two-block objective ζ(x, y).
x_domain (Sequence[ConstraintLayer]) – Inequality layers whose strict interior is the domain of x. May be empty (an unconstrained x-block). Stored as a tuple after construction, so a caller cannot mutate it and bypass the dimension validation.
y_domain (Sequence[ConstraintLayer]) – Inequality layers for the y-block. May be empty. Also stored as a tuple.
x_equality (EqualityConstraintLayer | None) – An optional linear equality A_eq x = b_eq on the x-block; the y-block has no equality. When present, the residual gains a multiplier term and a primal-feasibility row, and the Newton system is bordered (see Notes). When None, the problem is the plain two-block saddle. The multiplier λ changes from one Newton step to the next and belongs to the solver, not to this passive container, so the methods take it as an argument.

Notes

With an equality A_eq x = b_eq and multiplier λ, the residual extends to

r = [r_x; r_y; r_eq],

with

r_x = ∇_x ζ - (1/t) g_X + A_eq^T λ,
r_y = ∇_y ζ + (1/t) g_Y,
r_eq = A_eq x - b_eq,

and the Newton Jacobian becomes the equality-bordered coupled saddle

 _                       _
|  -M_xx   B      A_eq^T  |
|  B^T     H_yy   0       |
|_ A_eq    0      0      _|

an EqualityCoupledSaddleSystem rather than a CoupledSaddleSystem. An infeasible-start Newton method handles the equality naturally: x need not satisfy A_eq x = b_eq initially (the r_eq row drives the residual to zero), and λ starts at zero.

solve(solver, x0=None, y0=None)[source]#

Solve the problem with a two-block engine; shorthand for solver.solve.

Like Problem.solve: a caller states the two-block problem once and picks the algorithm by passing a TwoBlockStationaritySolver, the barrier StationaritySolver or the Mehrotra PrimalDualStationaritySolver.

Parameters:

solver (TwoBlockStationaritySolver) – The two-block engine to run.
x0 (npt.NDArray[np.float64] | None) – x-block warm start, forwarded to the engine.
y0 (npt.NDArray[np.float64] | None) – y-block warm start, forwarded to the engine.

Returns:

result – The engine’s result.

Return type:

TwoBlockSolveResult

Raises:

StationarityConvergenceError – Propagated from solver.solve when the solve does not converge.
ProblemInfeasibleError – Propagated from solver.solve when a block’s domain has no strictly interior point.

property num_equalities#

Count the equality rows p on the x-block.

Returns:: p – Number of rows of x_equality; zero when it is None.
Return type:: int

property dimension_x#

Return the dimension n of the x-block, taken from the objective.

Returns:: n – The objective’s dimension_x.
Return type:: int

property dimension_y#

Return the dimension m of the y-block, taken from the objective.

Returns:: m – The objective’s dimension_y.
Return type:: int

residual(x, y, t, lam=None)[source]#

Compute the barrier-smoothed joint stationarity residual.

Without an x_equality this is [r_x; r_y], with

r_x = ∇_x ζ - (1/t) g_X,   r_y = ∇_y ζ + (1/t) g_Y.

With one, it extends to [r_x; r_y; r_eq]: r_x gains the multiplier term A_eq^T λ, and the equality row is r_eq = A_eq x - b_eq.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n, strictly inside x_domain (not checked; see _barrier_gradient).
y (npt.NDArray[np.float64]) – The y-block, of length m, strictly inside y_domain (not checked).
t (float) – The inverse barrier weight, positive; the barrier enters with weight 1 / t.
lam (npt.NDArray[np.float64] | None) – The equality multiplier λ, of length p. Required when x_equality is set, and not allowed otherwise.

Returns:

residual – [r_x; r_y], of length n + m, or [r_x; r_y; r_eq], of length n + m + p, when x_equality is set.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When t ≤ 0, when lam is missing although x_equality is set (or supplied although it is not), or when lam has the wrong length.

build_coupled_system(x, y, t)[source]#

Assemble the barrier coupled Newton system for one step.

Without an x_equality this is the plain saddle

[[-M_xx, B], [B^T, H_yy]],

a CoupledSaddleSystem. With one it is the equality-bordered saddle

[[-M_xx, B, A_eq^T], [B^T, H_yy, 0], [A_eq, 0, 0]],

an EqualityCoupledSaddleSystem.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n, strictly inside x_domain.
y (npt.NDArray[np.float64]) – The y-block, of length m, strictly inside y_domain.
t (float) – The inverse barrier weight, positive; the barrier Hessians enter with weight 1 / t.

Returns:

system –

The coupled saddle with

M_xx = -∇^2_xx ζ + (1/t) H_X,
H_yy = ∇^2_yy ζ + (1/t) H_Y,
B = ∇^2_xy ζ,

bordered by A_eq when x_equality is set.

Return type:

CoupledSaddleSystem | EqualityCoupledSaddleSystem

Raises:

ValueError – When t ≤ 0, or when x or y is not strictly inside its domain (propagated from _barrier_hessian).

newton_step(x, y, t)[source]#

Compute one Newton step (Δx, Δy) for the barrier-smoothed residual.

Solves the coupled saddle against (-r_x, -r_y), the update that drives the residual to zero to first order.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n, strictly inside x_domain.
y (npt.NDArray[np.float64]) – The y-block, of length m, strictly inside y_domain.
t (float) – The inverse barrier weight, positive.

Returns:

delta_x (npt.NDArray[np.float64]) – The x-block step, of length n.
delta_y (npt.NDArray[np.float64]) – The y-block step, of length m.

Raises:

ValueError – When the problem has an x_equality: the equality multiplier carries over from step to step, so a stateless single step is ill-defined. Drive such a problem through StationaritySolver, which owns the multiplier, or assemble the step from residual and build_coupled_system. Also propagated from those two methods when t ≤ 0 or a block is not strictly inside its domain.
NewtonStepError – Propagated from CoupledSaddleSystem.solve when M_xx or the Schur complement cannot be solved.

primal_dual_residual(x, y, lam_x, lam_y, lam_eq=None)[source]#

Compute the primal-dual joint stationarity residual, with explicit duals.

The primal-dual counterpart of residual: the barrier’s smoothed gradient terms are replaced by the explicit constraint forces G^T λ built from the per-block inequality duals. Without an x_equality this is [r_x; r_y]; with one it extends to [r_x; r_y; r_eq], where r_x gains A_eq^T λ_eq and the equality row is r_eq = A_eq x - b_eq:

r_x = ∇_x ζ - G_X^T λ_x  (+ A_eq^T λ_eq),
r_y = ∇_y ζ + G_Y^T λ_y,

where G_X and G_Y are the x-domain and y-domain constraint Jacobians. The concave x-block subtracts its constraint force and the convex y-block adds it, matching the sign convention of residual. These are only the dual-feasibility (stationarity) rows; the complementarity rows S Λ 1 = μ 1 belong to the solver, which carries the slacks, duals and centering.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n.
y (npt.NDArray[np.float64]) – The y-block, of length m.
lam_x (npt.NDArray[np.float64]) – The x-domain inequality duals, one per x-domain constraint, stacked in layer order; empty when x_domain is empty.
lam_y (npt.NDArray[np.float64]) – The y-domain inequality duals, in the same convention.
lam_eq (npt.NDArray[np.float64] | None) – The equality multiplier λ_eq, of length p. Required when x_equality is set, and not allowed otherwise.

Returns:

residual – [r_x; r_y], of length n + m, or [r_x; r_y; r_eq], of length n + m + p, when x_equality is set.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When lam_x or lam_y has the wrong length, or when lam_eq is missing although x_equality is set (or supplied although it is not), or has the wrong length.

build_primal_dual_coupled_system(x, y, lam_x, lam_y)[source]#

Assemble the primal-dual coupled Newton system, with explicit duals.

The primal-dual counterpart of build_coupled_system: the weights the barrier method derives from t are replaced by the explicit duals in lam_x and lam_y, so the reduced Newton system has the same two-block saddle structure

[[-M_xx, B], [B^T, H_yy]]

(or its equality-bordered form), with

M_xx = -∇^2_xx ζ + ∑ λ_x,i C_x,i,
H_yy = ∇^2_yy ζ + ∑ λ_y,j C_y,j,
B = ∇^2_xy ζ,

where

C_i = ∇^2 f_i - (1 / f_i) ∇f_i ∇f_i^T

is the shared curvature atom. On the central path (λ_i = -1 / (t f_i)) this reproduces build_coupled_system.

The equality multiplier enters only the residual, not the Newton matrix (the border A_eq is structural), so this method, like build_coupled_system, takes no equality multiplier.

Parameters:

x (npt.NDArray[np.float64]) – The x-block, of length n, strictly inside x_domain (the atom’s 1 / f_i term requires it).
y (npt.NDArray[np.float64]) – The y-block, of length m, strictly inside y_domain.
lam_x (npt.NDArray[np.float64]) – The x-domain inequality duals, stacked in layer order as in primal_dual_residual; strictly positive for a genuine primal-dual iterate.
lam_y (npt.NDArray[np.float64]) – The y-domain inequality duals, in the same convention.

Returns:

system – The coupled saddle with M_xx, H_yy and B as above, bordered by A_eq when x_equality is set.

Return type:

CoupledSaddleSystem | EqualityCoupledSaddleSystem

Raises:

ValueError – When lam_x or lam_y has the wrong length, or when a block point is not strictly inside its domain.

exception cvxium.saddle.stationarity.StationarityConvergenceError(message, suboptimality, result)[source]#

Bases: OptimizationError

Error raised when a two-block stationarity solve fails to converge.

The stationarity engines (StationaritySolver and PrimalDualStationaritySolver) return a TwoBlockSolveResult only when the joint stationarity problem is solved to the user’s tolerance: an OPTIMAL result, or an ACCEPTABLE one at the soft tolerance. Every other termination (an exhausted outer or inner iteration budget, or a centering or Newton step that can make no further progress) is a non-convergence failure, and raises this exception rather than returning a status a caller might overlook.

This mirrors the barrier method’s InteriorPointMethodError: running out of iterations is never proof that the problem is infeasible, so it is never reported as an INFEASIBLE status. Infeasibility needs a genuine certificate, raised as ProblemInfeasibleError from Phase I.

Parameters:

message (str) – Description of the failure.
suboptimality (float) – The failure metric at the point of giving up: the barrier gap on outer exhaustion, or the stationarity residual or KKT score on a stalled step.
result (TwoBlockSolveResult) – The best-effort iterate reached before the solve gave up, with status SolveStatus.MAX_ITERATIONS. It did not reach the user’s tolerance but is often close, so a caller that repackages a terminal status, or continues a sweep past a hard subproblem, can catch the error and inspect result instead of re-running with a looser tolerance. Stored as result; its x-block also becomes the base class’s last_iterate.

class cvxium.saddle.stationarity.TwoBlockStationaritySolver(*args, **kwargs)[source]#

Bases: Protocol

Structural interface for a two-block stationarity solver.

A two-block StationarityProblem with a start (x0, y0) does not fit the single-Problem, single-x0 signature of Cvxium’s Solver, so the two-block engines, the barrier StationaritySolver and the primal-dual PrimalDualStationaritySolver, sit outside that hierarchy. This protocol names the shape they share, so a caller can hold either without knowing which engine it drives.

An implementer solves the joint stationarity problem from optional per-block warm starts and returns a TwoBlockSolveResult. With structural typing, an engine conforms by shape alone, with no base class and no registration. runtime_checkable also lets a caller test conformance with isinstance; that checks only that the method exists, not its signature, which is left to the static type checker.

solve(problem, x0=None, y0=None)[source]#

Solve a two-block stationarity problem from optional warm starts.

Parameters:

problem (StationarityProblem) – The two-block stationarity problem to solve.
x0 (npt.NDArray[np.float64] | None) – x-block warm start; an implementer may feasibilize it, or choose a default when it is None.
y0 (npt.NDArray[np.float64] | None) – y-block warm start, with the same semantics as x0.

Returns:

result – The converged iterate, work counts, per-block duals, convergence diagnostics and status, returned only when the solve reaches the user’s (soft) tolerance. An engine that computes the saddle value reports it, with its bracket, in the result’s objective.

Return type:

TwoBlockSolveResult

Raises:

StationarityConvergenceError – When the solve fails to converge (an exhausted budget or a stalled step); the best-effort iterate is on the exception’s result.
ProblemInfeasibleError – When a block’s domain has no strictly interior point.

Barrier Solver#

Infeasible-start Newton barrier solver for a two-block StationarityProblem.

StationaritySolver finds a joint stationary point (a saddle point, when the objective is concave in x and convex in y) of a StationarityProblem. The outer loop raises the barrier parameter t geometrically. At each t, an inner centering loop takes infeasible-start Newton steps on the barrier-smoothed residual, solving the coupled saddle system the problem assembles, and accepts each step by a residual-norm backtracking line search. A start that is not strictly interior is first made so by a per-block Phase I.

At each outer step the solver also brackets the saddle value V in two ways: a closed-form bracket that is exact on the central path, and a certified bracket from two per-block dual solves (_BlockDual) that holds off it too. _BarrierTrace accumulates both, with the per-step convergence trajectories, for the result’s diagnostics.

class cvxium.saddle.stationarity_solver.StationaritySolver(settings=None, check_inertia=True)[source]#

Bases: object

Infeasible-start Newton barrier solver for a two-block stationarity problem.

Drives a StationarityProblem (a smooth two-block objective ζ(x, y) with log-barrier domains) to a joint stationary point. The outer loop raises the barrier parameter t geometrically. Each inner centering solves the coupled Newton system

[[-M_xx, B], [B^T, H_yy]]

and takes an infeasible-start step that reduces ‖r‖: a cap that keeps both blocks inside their domains, then a residual-norm Armijo backtrack. The solve converges when the barrier gap n_barriers / t falls within

max(``settings.outer.tolerance``, ``settings.outer.tolerance_rel`` |ζ|),

an absolute floor, optionally relaxed to a fraction of the objective’s scale (see _gap_tolerance).

The class stands apart from Solver: a two-block StationarityProblem with a start (x0, y0) does not fit Solver.solve’s single-Problem, single-x0 signature, and the ‖r‖ centering has no single scalar objective to backtrack on (the same reasoning that keeps CoupledSaddleSystem off the SaddleSystem base).

Parameters:

settings (OptimizationSettings | None) – Algorithm settings; None uses OptimizationSettings().
check_inertia (bool) – Whether to check the coupled saddle’s block curvature signs once, at the first Newton step of the first outer iteration (see _check_inertia).

solve(problem, x0=None, y0=None)[source]#

Solve a two-block problem from a start (x0, y0), feasibilizing as needed.

A block’s start need not be strictly interior. A start already strictly inside its domain is used as given (the infeasible-start Newton loop drives any x_equality residual to zero), while a None or non-interior start triggers a per-block Phase I that produces a strictly interior point, respecting x_equality for the x-block (see _resolve_block). With no inequality constraints at all, the solver centers once and reports the exact saddle value.

Parameters:

problem (StationarityProblem) – The two-block stationarity problem.
x0 (npt.NDArray[np.float64] | None) – x-block start. None, or a non-interior guess, is feasibilized by Phase I, which a supplied guess seeds; when x_domain is empty, None becomes zeros.
y0 (npt.NDArray[np.float64] | None) – y-block start, with the same semantics as x0 (the y-block has no equality).

Returns:

result – The converged iterate, work counts, per-block duals, saddle value with its bracket, diagnostics and status (OPTIMAL or ACCEPTABLE).

Return type:

TwoBlockSolveResult

Raises:

StationarityConvergenceError – When the barrier loop fails to converge: centering stalls, or the outer budget runs out with the gap above tolerance. The best-effort iterate is on the exception’s result, with status MAX_ITERATIONS.
ProblemInfeasibleError – When Phase I finds that a block’s domain has no strictly interior point.
ValueError – When the inertia check fails, or when x0 or y0 does not match its block’s dimension.

Primal-Dual Solver#

Mehrotra predictor-corrector primal-dual solver for two-block stationarity problems.

Defines PrimalDualStationaritySolver, the primal-dual sibling of the barrier engine StationaritySolver. Both conform to the TwoBlockStationaritySolver protocol, solve(problem, x0, y0) -> TwoBlockSolveResult, so a downstream driver or frontier wrapper can hold either without knowing which engine it drives. Where the barrier engine nests an inner centering inside an outer loop that raises t, this engine follows the central path directly: it carries explicit per-block inequality duals λ^x and λ^y, and takes one Mehrotra predictor-corrector step per iteration on the coupled saddle, driving the surrogate gap

η = s^x · λ^x + s^y · λ^y

and the stationarity residual to zero together.

One iteration. The slacks s^x = -f^x(x) and s^y = -f^y(y) are implicit (Boyd and Vandenberghe, 2004, section 11.7). The stationarity residual (r_x, r_y) comes from StationarityProblem.primal_dual_residual, and the coupled saddle

[[-M_xx, B], [B^T, H_yy]]

from StationarityProblem.build_primal_dual_coupled_system. Eliminating the complementarity rows S Λ 1 = μ 1 gives the reduced right-hand side

b_x = -r_x - J_x^T (r_cx / s^x),   b_y = -r_y + J_y^T (r_cy / s^y),

which the saddle solves; the dual and slack steps follow as

Δλ^x = S_x^{-1} (-r_cx + Λ^x J_x Δx),   Δs^x = -J_x Δx

(mirrored for y). The concave x-block subtracts its constraint force and the convex y-block adds it, matching the sign convention of primal_dual_residual. An affine predictor (r_c = s ∘ λ) sets a per-block Mehrotra centering

σ_x = (μ^x_aff / μ_x)^3,   σ_y = (μ^y_aff / μ_y)^3,

where each block uses its own surrogate gap and σ is floored (_SIGMA_FLOOR), so a block-imbalanced saddle still recenters both blocks. The corrector recenters each block to its own σ μ and carries the second-order term Δs_aff ∘ Δλ_aff. One factorization backs both solves.

Equality path. A StationarityProblem with an x_equality A_eq x = b_eq takes the bordered path: the residual gains the primal-feasibility row r_eq = A_eq x - b_eq and an A_eq^T ν term in r_x, and the coupled saddle borders to EqualityCoupledSaddleSystem, whose solve(b_x, b_y, b_eq) returns the equality-multiplier step Δν alongside the primal steps. The reduced x and y right-hand sides are unchanged, and b_eq = -r_eq. The multiplier ν is free (unsigned), so it takes the x-block step length with no fraction-to-boundary cap or clamp.

There is no Levenberg damping (unlike the single-block PrimalDualSolver): QuadraticSaddleObjective and barrier-domain problems are well conditioned by construction, and a singular reduced saddle surfaces as a NewtonStepError from the coupled solve. The start-resolution and inertia helpers are module-level mirrors of StationaritySolver’s methods.

References

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

Mehrotra, Sanjay, On the Implementation of a Primal-Dual Interior Point Method, SIAM Journal on Optimization 2, 1992.

class cvxium.saddle.primal_dual_stationarity_solver.PrimalDualStationaritySolver(settings=None, check_inertia=True)[source]#

Bases: object

Mehrotra predictor-corrector primal-dual method for a two-block saddle.

Drives a StationarityProblem (a smooth two-block objective ζ(x, y) with log-barrier domains, and optionally an equality on x) to a joint stationary point by following the central path with explicit per-block inequality duals. Conforms to the TwoBlockStationaritySolver protocol; the module docstring describes one iteration.

Standalone rather than a Solver subclass, for the same reason as StationaritySolver: the two-block problem and its (x0, y0) start do not fit the single-Problem, single-x0 signature of Solver.solve.

Parameters:

settings (OptimizationSettings | None) – Algorithm settings; None uses OptimizationSettings(). The solver reads inner.max_iterations (iteration cap), outer.tolerance (strict tolerance on the surrogate gap and the residuals), outer.tolerance_soft (soft-tolerance fallback), pd_fraction_to_boundary_min (floor of the Mehrotra boundary factor), armijo.beta and armijo.min_step (endgame backtracking and step floor), trace, and observability.label (trace name). The settings also configure the Phase I BarrierSolver.
check_inertia (bool) – Whether to check the coupled saddle’s inertia at the first iteration (see _check_inertia).

solve(problem, x0=None, y0=None)[source]#

Solve a two-block stationarity problem from optional per-block warm starts.

A start already strictly interior to its domain is used as given; otherwise a per-block Phase I finds one (the x-block Phase I also lands on A_eq x = b_eq when an x_equality is set). The bordered _run_equality loop runs when an x_equality is present, and the plain _run loop otherwise.

Strict interiority does not imply a usable start. A warm start can satisfy every f < 0 and still hug one face so closely that the first Mehrotra step cannot fit between it and the boundary, so the solve fails without improving on its start. When a run fails that way (_made_no_progress) and some block’s start was the caller’s own, solve replaces those blocks’ starts with Phase I points (_centered_start) and runs once more. The retry is best-effort: if it fails in any way, the original failure is raised.

Parameters:

problem (StationarityProblem) – The two-block stationarity problem, optionally with an x_equality.
x0 (npt.NDArray[np.float64] | None) – Warm start for the x-block, of length n; Phase I runs when it is None or not strictly interior.
y0 (npt.NDArray[np.float64] | None) – Warm start for the y-block, of length m, with the same semantics.

Returns:

result – The iterate, iteration count, per-block duals, saddle value, surrogate-gap suboptimality and primal-dual diagnostics. Returned only with status OPTIMAL, or ACCEPTABLE when the best iterate met the soft tolerance.

Return type:

TwoBlockSolveResult

Raises:

StationarityConvergenceError – When the iteration fails to converge to the soft tolerance: an exhausted budget, a slack reaching its boundary, or a coupled system admitting no usable Newton step. The best iterate rides on the exception’s result (status MAX_ITERATIONS), and a NewtonStepError that triggered the stop is chained as the cause.
ProblemInfeasibleError – When a block’s domain admits no strictly interior point.
ValueError – When x0 or y0 has the wrong shape, or when the inertia check fails on the first run.