Cones#

Cones#

Symmetric cones and Nesterov-Todd scalings for Cvxium’s conic solvers.

Defines the cone primitives that the homogeneous self-dual solver (HsdSolver) and the second-order cone program helpers need to work over a product cone

K = K_1 × K_2 × ... × K_q.

Cone is the abstract base: the canonical interior point, membership tests, the step to the boundary, the Jordan product, the arrow-operator solve, and the Nesterov-Todd (NT) scaling W, both as an object exposing W-actions and as precomputed factors for the structured W^2 operations that assemble the Newton system. SecondOrderCone and Nonneg are the concrete cones.

NtScaling and its subclasses carry the NT scaling at one primal-dual pair (x, s): DenseNtScaling for an explicit matrix W, and two private factored forms that apply W, W^{-1}, and W^2 in time linear in the cone dimension. SocNtFactors and NonnegNtFactors hold the precomputed factors of W^2.

A new symmetric cone (for example, the positive semidefinite cone) is added by subclassing Cone, without changing the solvers.

class cvxium.problems.cone.Cone[source]#

Bases: ABC

Abstract symmetric cone supporting a conic interior-point method.

Concrete subclasses implement the cone-specific arithmetic. A conic solver works over a product cone

K = K_1 × K_2 × ... × K_q,

handling each block through this interface and otherwise staying cone-agnostic.

A symmetric cone is proper (closed, convex, pointed, with non-empty interior) and self-dual under the standard inner product, so the dual cone is K itself. Each symmetric cone carries a Jordan algebra whose identity element e_K is strictly interior and serves as the canonical interior point.

abstractmethod dimension()[source]#

Return the dimension of the cone’s embedding space.

Returns:

n – Length of an element vector.

Return type:

int

abstractmethod interior_point()[source]#

Return the canonical strictly interior point e_K.

The point is the identity element of the cone’s Jordan algebra. A conic solver uses it to build a strictly feasible starting iterate: both the primal point x and the dual slack s start at e_K, since the cone is self-dual.

Returns:

e – The identity element, of length n.

Return type:

npt.NDArray[np.float64]

abstractmethod is_in_cone(x, strict=False)[source]#

Test whether a vector lies in the cone.

Parameters:

  • x (npt.NDArray[np.float64]) – Vector of length n to test.

  • strict (bool) – When True, require x to lie in the open interior; otherwise the closed cone suffices.

Returns:

inside – True when x is in the cone (or in its interior, when strict).

Return type:

bool

Raises:

ValueError – When x does not have length n.

abstractmethod step_to_boundary(x, dx)[source]#

Find the largest step along a direction that stays in the cone.

Computes the largest α ≥ 0 such that x + α dx remains in the closed cone. A conic solver uses it to choose step lengths.

Parameters:

  • x (npt.NDArray[np.float64]) – Current point, in the cone.

  • dx (npt.NDArray[np.float64]) – Step direction, of length n.

Returns:

alpha – The largest feasible step; inf when x + α dx stays in the cone for every α ≥ 0.

Return type:

float

Raises:

ValueError – When x or dx does not have length n.

abstractmethod jordan_product(a, b)[source]#

Compute the Jordan product of two vectors in the cone’s embedding space.

Each symmetric cone has an associated Jordan algebra. Its identity element e_K (returned by interior_point) satisfies

e_K ∘ x = x

for all x. A conic solver measures complementarity through the per-block Jordan product x ∘ s.

Parameters:

  • a (npt.NDArray[np.float64]) – Left operand, of length n.

  • b (npt.NDArray[np.float64]) – Right operand, of length n.

Returns:

product – The Jordan product a ∘ b, of length n.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When a or b does not have length n.

abstractmethod nt_scaling(x, s)[source]#

Compute the Nesterov-Todd scaling at a primal-dual pair (x, s).

The NT scaling matrix W is symmetric positive definite and satisfies

W^2 s = x.

A conic solver uses W to symmetrize the Newton system and keep iterates well conditioned near the cone boundary.

Parameters:

  • x (npt.NDArray[np.float64]) – Primal point, strictly inside the cone.

  • s (npt.NDArray[np.float64]) – Dual slack, strictly inside the cone.

Returns:

scaling – An object that applies W, W^{-1}, and W^2 to vectors.

Return type:

NtScaling

Raises:

  • ValueError – When x or s does not have length n.

  • ConeBoundaryError – When x or s is not strictly inside the cone.

abstractmethod prepare_w_squared(x, s)[source]#

Precompute the Nesterov-Todd factors for the structured W^2 operations.

Callers compute the factors once per (x, s) pair and pass them to apply_w_squared and accumulate_normal_equations, avoiding a recomputation on every call. Each cone uses its own structured form of W^2, so the dense matrix is never formed.

Parameters:

  • x (npt.NDArray[np.float64]) – Primal point, strictly inside the cone.

  • s (npt.NDArray[np.float64]) – Dual slack, strictly inside the cone.

Returns:

factors – Cone-specific factors (SocNtFactors or NonnegNtFactors), opaque to the caller.

Return type:

object

Raises:

  • ValueError – When x or s does not have length n.

  • ConeBoundaryError – When x or s is not strictly inside the cone.

abstractmethod apply_w_squared(precomputed, v)[source]#

Multiply a vector by W^2 using precomputed Nesterov-Todd factors.

Uses the cone’s structured form of W^2 (rank-one plus diagonal for the second-order cone, diagonal for the nonnegative orthant), so the dense matrix is never formed. A conic solver uses it in the back-substitution step of the Newton system and to assemble per-block W^2 v products.

Parameters:

  • precomputed (object) – Factors returned by this cone’s prepare_w_squared.

  • v (npt.NDArray[np.float64]) – Vector of length n to multiply.

Returns:

W2v – The product W^2 v, of length n.

Return type:

npt.NDArray[np.float64]

abstractmethod accumulate_normal_equations(precomputed, A_block, c_block)[source]#

Compute one cone block’s contributions to the reduced Newton system.

Computes the three contractions

A_b W^2 A_b^T,   A_b W^2 c_b,   c_b^T W^2 c_b,

where A_b holds the block’s columns of the constraint matrix and c_b the block’s entries of the cost vector. The homogeneous self-dual solver reduces each Newton step to a (p + 1) × (p + 1) system whose coefficients sum these contributions over the blocks. Each cone computes them from its structured W^2 without forming it.

Parameters:

  • precomputed (object) – Factors returned by this cone’s prepare_w_squared.

  • A_block (npt.NDArray[np.float64]) – The block’s columns of the constraint matrix, of shape (p, n).

  • c_block (npt.NDArray[np.float64]) – The block’s entries of the cost vector, of length n.

Returns:

  • M_block (npt.NDArray[np.float64]) – A_b W^2 A_b^T, of shape (p, p).

  • AW2c_block (npt.NDArray[np.float64]) – A_b W^2 c_b, of length p.

  • cW2c_block (float) – c_b^T W^2 c_b.

abstractmethod solve_jordan_arrow(v, rhs)[source]#

Solve the arrow system L(v) z = rhs.

The arrow operator L(v) is the matrix of left Jordan multiplication by v,

L(v) y = v ∘ y.

A conic solver uses it to build the right-hand side of the centering step.

Parameters:

  • v (npt.NDArray[np.float64]) – Vector strictly inside the cone, so that L(v) is invertible.

  • rhs (npt.NDArray[np.float64]) – Right-hand side, of length n.

Returns:

z – The solution, of length n.

Return type:

npt.NDArray[np.float64]

class cvxium.problems.cone.SecondOrderCone(dim)[source]#

Bases: Cone

Second-order (Lorentz) cone of dimension n.

The cone is

K = { x : x_0 ≥ ‖x_{1:}‖ },

with the interior where the inequality is strict. Its Jordan product is

x ∘ y = (x^T y, x_0 y_{1:} + y_0 x_{1:}),

and the spectral decomposition of x has eigenvalues x_0 ± ‖x_{1:}‖. Throughout,

J = diag(1, -1, ..., -1),

and the Lorentz norm of x is sqrt(x^T J x).

Parameters:

dim (int) – Dimension n of the cone. Must be at least 2, or the constructor raises ValueError; the one-dimensional cone is Nonneg(1).

dimension()[source]#

Return the dimension of the cone.

Returns:

n – The dimension given at construction.

Return type:

int

interior_point()[source]#

Return the canonical interior point e_K = (1, 0, …, 0).

Returns:

e – The Jordan identity, of length n.

Return type:

npt.NDArray[np.float64]

is_in_cone(x, strict=False)[source]#

Test whether x_0 ≥ ‖x_{1:}‖, strictly when requested.

Parameters:

  • x (npt.NDArray[np.float64]) – Vector of length n to test.

  • strict (bool) – When True, require x_0 > ‖x_{1:}‖ (the open interior).

Returns:

inside – True when x is in the cone (or in its interior, when strict).

Return type:

bool

Raises:

ValueError – When x does not have length n.

step_to_boundary(x, dx)[source]#

Find the largest step along a direction that stays in the second-order cone.

Along the ray x + α dx, membership needs two conditions. The first is nonnegativity of the Lorentz quadratic

f(α) = (x_0 + α dx_0)^2 - ‖(x + α dx)_{1:}‖^2 = A α^2 + 2 B α + C,

where

A = dx_0^2 - ‖dx_{1:}‖^2,
B = x_0 dx_0 - x_{1:}^T dx_{1:},
C = x_0^2 - ‖x_{1:}‖^2 ≥ 0.

The second is the apex condition

x_0 + α dx_0 ≥ 0,

which the squared form alone misses, since f(α) ≥ 0 also holds on the reflected cone -K. The step is the smaller of the largest α satisfying each, found in closed form.

Parameters:

  • x (npt.NDArray[np.float64]) – Current point, in the cone (so that C ≥ 0).

  • dx (npt.NDArray[np.float64]) – Step direction, of length n.

Returns:

alpha – The largest feasible step; inf when x + α dx stays in the cone for every α ≥ 0.

Return type:

float

Raises:

ValueError – When x or dx does not have length n.

jordan_product(a, b)[source]#

Compute the second-order cone Jordan product.

The product is

a ∘ b = (a^T b, a_0 b_{1:} + b_0 a_{1:}).
Parameters:

  • a (npt.NDArray[np.float64]) – Left operand, of length n.

  • b (npt.NDArray[np.float64]) – Right operand, of length n.

Returns:

product – The Jordan product a ∘ b, of length n.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When a or b does not have length n.

nt_scaling(x, s)[source]#

Compute the factored Nesterov-Todd scaling for the second-order cone.

The returned scaling applies W, W^{-1}, and W^2 in O(n) time, rather than forming the dense n × n matrix and inverting and squaring it at O(n^3) cost.

Parameters:

  • x (npt.NDArray[np.float64]) – Primal point, strictly inside the cone.

  • s (npt.NDArray[np.float64]) – Dual slack, strictly inside the cone.

Returns:

scaling – The factored NT scaling at (x, s).

Return type:

NtScaling

Raises:

  • ValueError – When x or s does not have length n.

  • ConeBoundaryError – When x or s is not strictly inside the cone.

Notes

The construction follows (Andersen, Roos, and Terlaky, 2003) and (Sturm, 1999):

  1. Compute the Lorentz norms

    φ = sqrt(x^T J x),   ψ = sqrt(s^T J s),

    both positive in the strict interior.

  2. Normalize to unit Lorentz norm: x_hat = x / φ and s_hat = s / ψ.

  3. Set γ = sqrt((1 + x_hat^T s_hat) / 2).

  4. Form

    y = (x_hat + J s_hat) / (2 γ),

    the unit-Lorentz-norm Jordan square of the scaling vector w_bar.

  5. Take the Jordan square root w_bar of y (see _jordan_sqrt).

  6. Set

    W = sqrt(φ / ψ) Q(w_bar),

    where Q(w) = 2 w w^T - det(w) J is the quadratic representation, and det(w_bar) = 1 by the normalization.

The identity W^2 s = x holds because Q(w_bar) is an automorphism of the cone, so it preserves Lorentz norms: W^2 scales them by φ / ψ, taking the norm of s to that of x, while w_bar aligns the directions.

References

Andersen, Erling D., Roos, Cornelis and Terlaky, Tamás, On Implementing a Primal-Dual Interior-Point Method for Conic Quadratic Optimization, Mathematical Programming 95, 2003.

Sturm, Jos F., Using SeDuMi 1.02, a MATLAB Toolbox for Optimization over Symmetric Cones, Optimization Methods and Software 11, 1999.

factored_w_squared(x, s)[source]#

Compute the factored form of W^2 for the second-order cone.

Returns the factors of

W^2 = scale (2 y y^T - J),

where

scale = φ / ψ,   y = (x_hat + J s_hat) / (2 γ)

are the Nesterov-Todd intermediates defined in nt_scaling. Because det(y) = 1 by construction, the quadratic representation simplifies to

Q(y) = 2 y y^T - J,

and the Jordan algebra identity Q(w)^2 = Q(w ∘ w) gives

W^2 = scale Q(w_bar)^2 = scale Q(y).

So W^2, unlike W itself, needs no Jordan square root.

Parameters:

  • x (npt.NDArray[np.float64]) – Primal point, strictly inside the cone.

  • s (npt.NDArray[np.float64]) – Dual slack, strictly inside the cone.

Returns:

  • scale (float) – The ratio φ / ψ of the Lorentz norms of x and s.

  • y (npt.NDArray[np.float64]) – The unit-Lorentz-norm vector of length n defining the rank-one term.

Raises:

  • ValueError – When x or s does not have length n.

  • ConeBoundaryError – When x or s is not strictly inside the cone.

prepare_w_squared(x, s)[source]#

Precompute the second-order cone factors of W^2.

Packages the result of factored_w_squared.

Parameters:

  • x (npt.NDArray[np.float64]) – Primal point, strictly inside the cone.

  • s (npt.NDArray[np.float64]) – Dual slack, strictly inside the cone.

Returns:

factors – The factors (scale, y), with W^2 = scale (2 y y^T - J).

Return type:

SocNtFactors

Raises:

  • ValueError – When x or s does not have length n.

  • ConeBoundaryError – When x or s is not strictly inside the cone.

apply_w_squared(precomputed, v)[source]#

Multiply a vector by W^2 using the factored second-order cone form.

Computes

W^2 v = scale (2 (y^T v) y - J v)

in O(n) time.

Parameters:

  • precomputed (object) – The SocNtFactors returned by prepare_w_squared.

  • v (npt.NDArray[np.float64]) – Vector of length n to multiply.

Returns:

W2v – The product W^2 v, of length n.

Return type:

npt.NDArray[np.float64]

accumulate_normal_equations(precomputed, A_block, c_block)[source]#

Compute the second-order cone block’s contributions to the reduced system.

With W^2 = scale (2 y y^T - J),

A_b W^2 A_b^T = scale (2 (A_b y) (A_b y)^T - a_0 a_0^T + A_r A_r^T),

where a_0 is the first column of A_b and A_r holds the remaining columns; the contractions with c_b follow the same pattern. The cost is O(p^2 n) per block, dominated by A_r A_r^T, and avoids forming W^2 and the O(p n^2) product A_b W^2.

Parameters:

  • precomputed (object) – The SocNtFactors returned by prepare_w_squared.

  • A_block (npt.NDArray[np.float64]) – The block’s columns of the constraint matrix, of shape (p, n).

  • c_block (npt.NDArray[np.float64]) – The block’s entries of the cost vector, of length n.

Returns:

  • M_block (npt.NDArray[np.float64]) – A_b W^2 A_b^T, of shape (p, p).

  • AW2c_block (npt.NDArray[np.float64]) – A_b W^2 c_b, of length p.

  • cW2c_block (float) – c_b^T W^2 c_b.

solve_jordan_arrow(v, rhs)[source]#

Solve the second-order cone arrow system L(v) z = r in O(n) time.

The arrow matrix is

        _               _
L(v) = |  v_0   v_{1:}^T |
       |_ v_{1:}  v_0 I _|

Eliminating the lower-right block v_0 I (a Schur complement) gives a closed form that avoids a dense O(n^3) solve. With the cone determinant

d = v_0^2 - ‖v_{1:}‖^2,

the solution is

z_0    = (v_0 r_0 - v_{1:}^T r_{1:}) / d,
z_{1:} = r_{1:} / v_0 - v_{1:} r_0 / d + v_{1:} (v_{1:}^T r_{1:}) / (v_0 d).
Parameters:

  • v (npt.NDArray[np.float64]) – Vector strictly inside the cone, so that d > 0 and v_0 > 0.

  • rhs (npt.NDArray[np.float64]) – Right-hand side r, of length n.

Returns:

z – The solution, of length n.

Return type:

npt.NDArray[np.float64]

class cvxium.problems.cone.NtScaling[source]#

Bases: ABC

Abstract Nesterov-Todd scaling at one primal-dual pair (x, s).

Declares the W-action interface that every concrete scaling shares. Each subclass chooses its own representation (an explicit matrix or an O(n) factored form) and supplies its own constructor; this base defines none, so the subclass signatures are not incompatible overrides. Every implementation satisfies

W^2 s = x

for the pair (x, s) at which it was constructed.

abstractmethod apply_W(v)[source]#

Multiply a vector by W.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

Wv – The product W v.

Return type:

npt.NDArray[np.float64]

abstractmethod apply_W_inv(v)[source]#

Multiply a vector by W^{-1}.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

Winv_v – The product W^{-1} v.

Return type:

npt.NDArray[np.float64]

abstractmethod apply_W_squared(v)[source]#

Multiply a vector by W^2.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

W2v – The product W^2 v.

Return type:

npt.NDArray[np.float64]

class cvxium.problems.cone.DenseNtScaling(W)[source]#

Bases: NtScaling

Nesterov-Todd scaling stored as an explicit matrix.

Stores W as a dense matrix and precomputes W^{-1} and W^2 at construction, so each apply is a matrix-vector product. Construction costs O(n^3), so the cones in this module return factored scalings instead; this class serves a caller that already holds an explicit W.

Parameters:

W (npt.NDArray[np.float64]) – Symmetric positive definite scaling matrix with W^2 s = x for the intended pair (x, s). The constructor raises ValueError when W is not square, and NumPy’s LinAlgError when it is singular.

apply_W(v)[source]#

Multiply a vector by W.

A matrix-vector product with the matrix stored at construction.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

Wv – The product W v.

Return type:

npt.NDArray[np.float64]

apply_W_inv(v)[source]#

Multiply a vector by W^{-1}.

A matrix-vector product with the matrix stored at construction.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

Winv_v – The product W^{-1} v.

Return type:

npt.NDArray[np.float64]

apply_W_squared(v)[source]#

Multiply a vector by W^2.

A matrix-vector product with the matrix stored at construction.

Parameters:

v (npt.NDArray[np.float64]) – Vector to multiply, of the cone’s dimension.

Returns:

W2v – The product W^2 v.

Return type:

npt.NDArray[np.float64]