Structured Linear Algebra#

Matrix Systems#

Composable structured linear systems for Cvxium’s Newton steps.

A MatrixSystem represents a structured symmetric positive-definite matrix H that can solve H x = b and multiply H y efficiently by exploiting its structure: zero, diagonal, banded, dense, low-rank updated, block diagonal, permuted, or arrow (a block bordered by a few dense rows and columns). MatrixSystems compose: a low-rank update can sit on top of any base, the block-diagonal, permuted and arrow systems take other MatrixSystems as parts, and plus sums two systems while keeping their structure, which is how Problem.build_hessian_system assembles a Hessian from its contributions. The arithmetic lives in the kernels of numerical_helpers; this module arranges them into a type-safe tree.

The indefinite KKT matrix

 _        _
|  H  A^T  |
|_ A   0  _|

of an equality-constrained Newton step is not a MatrixSystem (it is not positive definite, and has no multiply or low_rank_update). SaddleSystem and its subclasses solve it, holding a MatrixSystem for the H block; a SaddleFactory chooses the subclass for an equality constraint. CoupledSaddleSystem and EqualityCoupledSaddleSystem solve the two-block saddle-point systems of a stationarity solver.

class cvxium.linalg.systems.Definiteness(*values)[source]#

Bases: Enum

Declared definiteness of a MatrixSystem.

A MatrixSystem is positive definite by contract, the shape the barrier and primal-dual solvers factor. Maximizing a concave objective by minimizing its negation produces a negative-definite Hessian contribution, which MatrixSystem.negated declares by flipping this value.

This is a declaration (the contract a system advertises), distinct from the numerical check MatrixSystem.is_positive_definite.

property negated#

Return the definiteness of -H, given that of H.

Returns:: definiteness – Positive and negative definite swap; indefinite stays indefinite.
Return type:: Definiteness

class cvxium.linalg.systems.MatrixSystem[source]#

Bases: ABC

Structured symmetric positive-definite matrix that can solve and multiply itself.

Each subclass stores its matrix H in a structured form and implements solve, multiply and scaled in that form. The base supplies plus (a structured sum, through each shape’s _absorb), low_rank_update, negated, and a dense is_positive_definite fallback.

abstract property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

abstractmethod solve(b)[source]#

Solve H x = b.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side, a vector of length n or a matrix with n rows; a matrix is solved column by column.
Returns:: x – The solution, of the same shape as b.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – Implementations raise when H is not numerically positive definite.

abstractmethod multiply(y)[source]#

Compute the product H y.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows; a matrix is multiplied column by column.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

abstractmethod scaled(t)[source]#

Return t H, the matrix scaled by a nonnegative factor.

Scaling preserves structure. It weights a Hessian contribution, for example the objective’s curvature by the barrier parameter, when the Newton system is assembled.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – A system of the same structure representing t H.
Return type:: MatrixSystem

plus(other)[source]#

Return H + other, the structured sum of two systems.

Assembling a Hessian from its contributions is repeated plus. A ZeroSystem is the additive identity; a low-rank-updated operand re-nests its update on top of the summed bases; two pure bases combine through _absorb, tried in both orders so the sum does not depend on operand order. Two bases of the same shape compose into that shape (diagonal plus diagonal is diagonal). This is closure within a shape, not a general cross-shape combine.

Parameters:: other (MatrixSystem) – The system to add, of the same order.
Returns:: system – The sum, in the shared structure.
Return type:: MatrixSystem
Raises:: IncompatibleStructureError – When the orders differ, or the two bases cannot be combined without densifying.

low_rank_update(kappa, d=None, coalesce=True)[source]#

Return the system H + kappa diag(d) kappa^T.

Defined once here and inherited by every MatrixSystem, so a low-rank update can sit on top of any structure. The result solves through the Woodbury identity.

Parameters:

kappa (npt.NDArray[np.float64]) – A vector of length n (rank one) or a matrix of shape (n, p) (rank p).
d (npt.NDArray[np.float64] | float | None) – A scalar or a vector of length p of weights: positive entries are updates, negative entries downdates; None means unit weights. Zero-weight columns are dropped.
coalesce (bool) – Ignored here, and accepted for a signature uniform with LowRankUpdatedSystem.low_rank_update: a first update always nests a fresh layer, since there is nothing to merge into.

Returns:

system – A LowRankUpdatedSystem wrapping this system, or this system itself when every weight is zero.

Return type:

MatrixSystem

Raises:

ValueError – Propagated from _normalize_low_rank when kappa or d is malformed.

is_positive_definite()[source]#

Check whether the matrix is numerically positive definite.

A MatrixSystem is positive definite by contract, but curvature assembled from caller-supplied signs can violate that. This check lets a consumer, such as the stationarity solvers’ check on each diagonal block, fail cleanly on wrong-sign input without densifying the block by hand.

The default forms the dense matrix by multiplying the identity (n products) and attempts a Cholesky factorization, O(n^3). Structured subclasses override it with a cheaper test that reuses the factorization they cache.

Returns:: is_pd – True when the Cholesky factorization succeeds.
Return type:: bool

property definiteness#

Return the declared definiteness of the matrix, positive definite by default.

This is the contract the system advertises, distinct from the numerical check is_positive_definite. Every concrete base in this module is positive definite; a _NegatedSystem flips the declaration, so a concave objective’s Hessian can be carried as negative definite until NegatedObjective negates it back.

Returns:: definiteness – Definiteness.POSITIVE_DEFINITE in the default.
Return type:: Definiteness

negated()[source]#

Return -H as a sign-carrying view, flipping the declared definiteness.

scaled takes only a nonnegative factor, so negation is the one sign change it does not cover. It exists to maximize a concave objective g by minimizing -g: the objective advertises its negative-definite Hessian as H.negated(), and NegatedObjective negates it once more, recovering the positive-definite H the solver factors.

The default wraps this system in a lazy _NegatedSystem, whose own negated collapses the double negation back to this system, so a maximize round trip never asks a factorization to handle a negative-definite matrix.

Returns:: system – A system representing -H.
Return type:: MatrixSystem

class cvxium.linalg.systems.ZeroSystem(dimension)[source]#

Bases: MatrixSystem

Zero matrix of a given order, the additive identity for plus.

A component with no base curvature contributes a ZeroSystem: a linear objective, or an inequality layer whose barrier curvature is purely a low-rank update (the update sits on the zero base through low_rank_update). It multiplies to zero and drops out of plus. The zero matrix is singular, so solve raises; a ZeroSystem reaches a solve only after being summed with a positive-definite term.

Parameters:: dimension (int) – Order of the matrix.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Refuse to solve, since the zero matrix is singular.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side; unused.
Returns:: x – Never returned.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – Always.

multiply(y)[source]#

Compute 0 y, the zero array shaped like y.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows.
Returns:: zeros – Zeros of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t 0 = 0, this same system.

Parameters:: t (float) – Scale factor; unused.
Returns:: system – This system.
Return type:: MatrixSystem

negated()[source]#

Return -0 = 0, this same system.

Overrides the lazy default so that the negated additive identity stays the additive identity and keeps dropping out of plus (a _NegatedSystem wrapping zero would not). A negated linear objective’s zero curvature relies on this during assembly.

Returns:: system – This system.
Return type:: MatrixSystem

is_positive_definite()[source]#

Check positive definiteness, true only for the empty (order 0) matrix.

Returns:: is_pd – True when the order is 0.
Return type:: bool

class cvxium.linalg.systems.DiagonalSystem(eta, eta_inverse=None)[source]#

Bases: MatrixSystem

Diagonal matrix H = diag(eta).

Positive definiteness is not checked at construction: a diagonal with zero or negative entries may be built and composed through plus, and only solve enforces it. That lets a component declare a zero-curvature coordinate, such as a slack variable with no objective curvature of its own, that a later barrier term makes positive in the assembled sum.

Parameters:

eta (npt.NDArray[np.float64]) – The diagonal entries, a 1-D array of length n, or the constructor raises ValueError.
eta_inverse (npt.NDArray[np.float64] | None) – The reciprocals 1 / eta, for a caller that can compute them more stably than dividing. For example, a barrier term with curvature 1 / (x - lb)^2 already holds the slack s = x - lb and can pass s^2, avoiding two round trips through tiny quantities. It is trusted, not checked against 1 / eta, which would defeat its purpose; it must have the shape of eta, or the constructor raises ValueError. None computes 1 / eta (inf for a zero entry).

property eta#

Return the diagonal entries of the matrix.

Returns:: eta – The diagonal, of length n.
Return type:: npt.NDArray[np.float64]

property eta_inverse#

Return the reciprocal diagonal, carrying any stable value supplied at build.

Exposed so a caller that reorders a diagonal (for example, under a coordinate permutation) can carry the stable reciprocal forward rather than recompute it from eta and lose precision near a barrier boundary.

Returns:: eta_inverse – The reciprocals of the diagonal entries, of length n.
Return type:: npt.NDArray[np.float64]

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Solve diag(eta) x = b.

Positive definiteness is enforced here, not at construction. The check is on eta rather than eta_inverse, because a zero entry’s reciprocal is inf, which the kernel’s own check (eta_inverse > 0) would let through.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector of length n or a matrix with n rows.

Returns:

x – The solution eta_inverse ∘ b, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

NewtonStepError – When some entry of eta is not positive, or, propagated from solve_diagonal_eta_inverse, some caller-supplied reciprocal is not positive.
ValueError – Propagated from solve_diagonal_eta_inverse when b does not have n rows.

multiply(y)[source]#

Compute diag(eta) y.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H, scaling the diagonal and carrying the stable reciprocal.

Parameters:: t (float) – Positive scale factor; t = 0 makes the reciprocals infinite.
Returns:: system – A DiagonalSystem with diagonal t eta, its reciprocals divided by t.
Return type:: MatrixSystem

is_positive_definite()[source]#

Check positive definiteness: every diagonal entry strictly positive.

Returns:: is_pd – True when every entry of eta is positive.
Return type:: bool

class cvxium.linalg.systems.BandedSystem(ab, lower=False)[source]#

Bases: MatrixSystem

Symmetric banded matrix in SciPy’s compact banded storage.

The Cholesky factor is computed on the first solve (or is_positive_definite) and cached; the matrix is immutable, since every transform returns a new instance, so the factor cannot go stale and the several solves of one Newton step share it.

Parameters:

ab (npt.NDArray[np.float64]) – The band storage, of shape (p + 1, n) for bandwidth p, in the convention of scipy.linalg.cholesky_banded (see numerical_helpers.solve_banded). Must be 2-D, or the constructor raises ValueError.
lower (bool) – True when ab stores the lower triangle, False for the upper.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Solve H x = b with the cached banded Cholesky factor.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side, a vector of length n or a matrix with n rows.
Returns:: x – The solution, of the same shape as b.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – Propagated from factor_banded when H is not positive definite or has non-finite entries.

multiply(y)[source]#

Compute H y, exploiting the band structure.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H, preserving the band structure.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – A BandedSystem storing t ab.
Return type:: MatrixSystem

is_positive_definite()[source]#

Check positive definiteness by attempting the banded Cholesky factorization.

Reuses and populates the cached factor, so a later solve pays nothing.

Returns:: is_pd – True when the factorization succeeds.
Return type:: bool

class cvxium.linalg.systems.DenseSystem(matrix)[source]#

Bases: MatrixSystem

Dense symmetric positive-definite matrix.

A fallback for blocks with no exploitable structure. The Cholesky factor is computed on the first solve (or is_positive_definite) and cached for the remaining solves of the same instance, at O(n^3) for the factorization; prefer a structured MatrixSystem where one applies.

Parameters:: matrix (npt.NDArray[np.float64]) – The square symmetric matrix, of shape (n, n). The constructor raises ValueError when it is not square, and NewtonStepError when it has non-finite entries.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Solve H x = b with the cached dense Cholesky factor.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side, a vector of length n or a matrix with n rows.
Returns:: x – The solution, of the same shape as b.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – When the Cholesky factorization fails, so H is not positive definite.

multiply(y)[source]#

Compute H y without factoring.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – A DenseSystem storing t H.
Return type:: MatrixSystem

is_positive_definite()[source]#

Check positive definiteness by attempting the dense Cholesky factorization.

Reuses and populates the cached factor, so a later solve pays nothing.

Returns:: is_pd – True when the factorization succeeds (or a factor is already cached).
Return type:: bool

class cvxium.linalg.systems.BlockDiagonalSystem(blocks)[source]#

Bases: MatrixSystem

Block-diagonal matrix whose diagonal blocks are themselves MatrixSystems.

Solving and multiplying dispatch block by block, so the blocks may exploit entirely different structures.

Parameters:: blocks (Sequence[MatrixSystem]) – The diagonal blocks, in order. Must be non-empty, or the constructor raises ValueError.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Solve H x = b one diagonal block at a time.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side, a vector of length n or a matrix with n rows.
Returns:: x – The solution, of the same shape as b.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – Propagated from a block’s solve.

multiply(y)[source]#

Compute H y one diagonal block at a time.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H, scaling each diagonal block.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – A BlockDiagonalSystem of the scaled blocks.
Return type:: MatrixSystem

is_positive_definite()[source]#

Check positive definiteness: every diagonal block positive definite.

Returns:: is_pd – True when every block is positive definite.
Return type:: bool

class cvxium.linalg.systems.PermutedSystem(inner, to_grouped, to_interleaved=None)[source]#

Bases: MatrixSystem

Matrix whose exploitable structure appears only in a reordered basis.

Structure is a property of a matrix together with a variable ordering. A Hessian can be block diagonal in one ordering and merely banded, or dense-looking, in another, and the decision vector’s layout is fixed by the problem rather than chosen for the linear algebra. A typical case is a decision matrix whose row-major flattening interleaves its columns: a column-separable Hessian is block diagonal only once each column’s strided indices are gathered (grouped).

inner is the system in the grouped basis, and the wrapper presents it in the caller’s interleaved basis, representing

H = P^T H_grouped P,

where P gathers coordinates (P v = v[to_grouped]), so that inner equals H[np.ix_(to_grouped, to_grouped)]. Because it presents itself in the decision-variable basis, it is an ordinary MatrixSystem: a Problem sums it with the constraint layers’ contributions correctly, and no permutation-aware solver is needed. solve and multiply accept a matrix of right-hand sides, which the Woodbury path needs, so a low-rank update composes on top of a permuted base.

Parameters:

inner (MatrixSystem) – The system in the grouped (structured) basis.
to_grouped (npt.NDArray[np.intp]) – A permutation of range(n), n the order of inner, gathering interleaved coordinates into the grouped basis. The constructor raises ValueError when it is not such a permutation.
to_interleaved (npt.NDArray[np.intp] | None) – The inverse permutation, or None to derive it. The constructor raises ValueError when a supplied one is not the inverse of to_grouped.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

solve(b)[source]#

Solve H x = b: gather, solve in the grouped basis, scatter back.

Parameters:: b (npt.NDArray[np.float64]) – Right-hand side in the interleaved basis, a vector of length n or a matrix with n rows.
Returns:: x – The solution in the interleaved basis, of the same shape as b.
Return type:: npt.NDArray[np.float64]
Raises:: NewtonStepError – Propagated from the inner system’s solve.

multiply(y)[source]#

Compute H y: gather, multiply in the grouped basis, scatter back.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length n or a matrix with n rows, in the interleaved basis.
Returns:: Hy – The product in the interleaved basis, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H; scaling commutes with the permutation.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – A PermutedSystem over the scaled inner system.
Return type:: MatrixSystem

property definiteness#

Return the inner system’s declared definiteness; permuting preserves it.

Returns:: definiteness – The inner system’s definiteness.
Return type:: Definiteness

is_positive_definite()[source]#

Check positive definiteness of the inner system, which reordering preserves.

Returns:: is_pd – The inner system’s answer.
Return type:: bool

class cvxium.linalg.systems.ArrowSystem(upper_left, border, corner)[source]#

Bases: MatrixSystem

Bordered (arrow) matrix: a structured block with a small dense fringe.

The matrix is

      _            _
H =  |  A11    A12  |
     |_ A12^T  A22 _|

where the upper-left block A11 is any MatrixSystem of order M, and the border A12 (M × p) and corner A22 (p × p) are a small dense fringe of p rows and columns. Solving eliminates A11 through its Schur complement. When A11 is a DiagonalSystem and p = 1, the matrix has a literal arrow sparsity pattern, and a dedicated, numerically stable kernel is used.

The whole matrix is assumed positive definite; the indefinite saddle-point case (a zero corner) is a SaddleSystem.

Parameters:

upper_left (MatrixSystem) – The block A11, of order M.
border (npt.NDArray[np.float64]) – The block A12: a vector of length M (for p = 1) or a matrix of shape (M, p). The constructor raises ValueError when it does not have M rows.
corner (npt.NDArray[np.float64] | float) – The block A22: a scalar (for p = 1) or a matrix of shape (p, p). The constructor raises ValueError when its shape does not match the border.

property dimension#

Return the order of the matrix.

Returns:: n – The number of rows (and columns) of the square matrix.
Return type:: int

property upper_left#

Return the upper-left block A11.

Returns:: upper_left – The block A11, of order M.
Return type:: MatrixSystem

property border#

Return the border block A12.

Returns:: border – The block A12, of shape (M, p).
Return type:: npt.NDArray[np.float64]

property corner#

Return the corner block A22.

Returns:: corner – The block A22, of shape (p, p).
Return type:: npt.NDArray[np.float64]

diagonal_plus_low_rank()[source]#

Express the arrow as a diagonal plus a rank-two term, when possible.

When A11 is a DiagonalSystem and the corner is a single positive finite scalar c, the arrow is

H = D_full + K diag(1, -1) K^T,

where D_full is the diagonal of A11 with c appended, and the rank-two term encodes the border u. With u_ext = [u; 0] of length M + 1 and e the last unit vector,

 _        _
|  0    u  |  =  u_ext e^T + e u_ext^T  =  K diag(1, -1) K^T,
|_ u^T  0 _|
K = [u_ext + e, u_ext - e] / sqrt(2).

Returns None instead when A11 is not diagonal, p > 1, or the corner is not positive and finite.

Returns:

D_full (DiagonalSystem) – The diagonal of A11 with c appended, of order M + 1, carrying the stable reciprocals.
K (npt.NDArray[np.float64]) – The rank-two factor, of shape (M + 1, 2).
d_lr (npt.NDArray[np.float64]) – The weights [1, -1].

solve(b)[source]#

Solve H x = b through the Schur complement on the upper-left block.

Uses the arrow-sparsity kernel for a diagonal A11 with p = 1, the block-plus-one kernel for another A11 with p = 1, and the general Schur-complement kernel otherwise.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector of length M + p or a matrix with M + p rows.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

NewtonStepError – Propagated from the kernels or the upper-left solve when H is not positive definite.
ValueError – Propagated from the kernels when b has the wrong number of rows.

multiply(y)[source]#

Compute H y block by block.

Parameters:: y (npt.NDArray[np.float64]) – A vector of length M + p or a matrix with M + p rows.
Returns:: Hy – The product, of the same shape as y.
Return type:: npt.NDArray[np.float64]

scaled(t)[source]#

Return t H, scaling the upper-left block, border and corner.

Parameters:: t (float) – Nonnegative scale factor.
Returns:: system – An ArrowSystem of the scaled blocks.
Return type:: MatrixSystem

Saddle Systems#

class cvxium.linalg.systems.SaddleSystem[source]#

Bases: ABC

Indefinite saddle-point system of an equality-constrained Newton step.

The matrix is

 _        _
|  H  A^T  |
|_ A   0  _|

with the Hessian H a MatrixSystem of order M and A the p × M equality-constraint matrix. A Newton step from an equality-feasible iterate has a zero lower right-hand side; an infeasible-start (primal-dual) step carries the equality residual there.

Not a MatrixSystem: the matrix is indefinite, so the positive-definite contract (multiply, low_rank_update, plus) does not apply. It holds a MatrixSystem for the H block, so that block keeps all its structure.

This is the abstract base. SchurComplementSaddleSystem eliminates H through dense Schur-complement algebra; SparseSaddleSystem and its subclasses handle a sparse A. The Newton, barrier and primal-dual solvers all solve their KKT systems through this interface.

abstract property hessian#

Return the Hessian block H.

Returns:: hessian – The H block, of order M.
Return type:: MatrixSystem

abstract property num_equality_constraints#

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

Returns:: p – The number of equality constraints.
Return type:: int

abstractmethod solve(b_top, b_bottom=None)[source]#

Solve the saddle-point system.

Solves

 _        _     _         _     _          _
|  H  A^T  |   |  delta_x  |   |  b_top     |
|_ A   0  _|   |_   nu    _| = |_ b_bottom _|

for the primal step and the equality multipliers.

Parameters:

b_top (npt.NDArray[np.float64]) – Top block of the right-hand side (the negative gradient or dual residual): a vector of length M, or a matrix of shape (M, k) solved column by column.
b_bottom (npt.NDArray[np.float64] | None) – Bottom block (the equality residual), of length p. None means zero, the equality-feasible Newton step; a nonzero value arises in primal-dual methods, whose iterates need not satisfy A x = b mid-iteration.

Returns:

delta_x (npt.NDArray[np.float64]) – The primal step, the top block of the solution.
nu (npt.NDArray[np.float64]) – The equality multipliers, the bottom block of the solution.

Raises:

NewtonStepError – Implementations raise when H, or the reduced system, cannot be solved.

class cvxium.linalg.systems.SchurComplementSaddleSystem(hessian, A)[source]#

Bases: SaddleSystem

Saddle-point system solved through the Schur complement on H.

solve eliminates H with solve_kkt_system, touching the Hessian only through H.solve: p + 1 solves with H, O(p^2 M) to form the p × p Schur complement, and O(p^3) to factor it.

Parameters:

hessian (MatrixSystem) – The Hessian block H, of order M.
A (npt.NDArray[np.float64]) – The equality-constraint matrix, of shape (p, M). The constructor raises ValueError when it is not 2-D or its column count is not M.

property hessian#

Return the Hessian block H.

Returns:: hessian – The H block, of order M.
Return type:: MatrixSystem

property num_equality_constraints#

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

Returns:: p – The number of equality constraints.
Return type:: int

property A#

Return the equality-constraint matrix A.

Returns:: A – The matrix A, of shape (p, M).
Return type:: npt.NDArray[np.float64]

solve(b_top, b_bottom=None)[source]#

Solve the saddle-point system through the Schur complement on H.

Solves

 _        _     _         _     _          _
|  H  A^T  |   |  delta_x  |   |  b_top     |
|_ A   0  _|   |_   nu    _| = |_ b_bottom _|

for the primal step and the equality multipliers.

Parameters:

b_top (npt.NDArray[np.float64]) – Top block of the right-hand side (the negative gradient or dual residual): a vector of length M, or a matrix of shape (M, k) solved column by column.
b_bottom (npt.NDArray[np.float64] | None) – Bottom block (the equality residual), of length p. None means zero, the equality-feasible Newton step; a nonzero value arises in primal-dual methods, whose iterates need not satisfy A x = b mid-iteration.

Returns:

delta_x (npt.NDArray[np.float64]) – The primal step, the top block of the solution.
nu (npt.NDArray[np.float64]) – The equality multipliers, the bottom block of the solution.

Raises:

NewtonStepError – Propagated from solve_kkt_system: H is not positive definite or is singular, or the reduced system has no solution.
ValueError – Propagated from solve_kkt_system when b_top or b_bottom has the wrong length.

class cvxium.linalg.systems.SparseSaddleSystem(hessian, A)[source]#

Bases: SaddleSystem

Saddle-point system for a sparse equality-constraint matrix A.

Like SchurComplementSaddleSystem, this eliminates H through the Schur complement

S = A H^{-1} A^T,

but it keeps A in CSR form and never forms the dense p × M constraint matrix. When H is a DiagonalSystem (the Hessian of a barrier-regularized linear program), H^{-1} A^T is A^T with its rows rescaled, so S is assembled entirely with sparse products and nothing of size p M is allocated. That makes the system usable when A has millions of rows and columns but only a handful of nonzeros per column, as in a large bipartite network-flow linear program. When H is an ArrowSystem with a diagonal core and a single border column, a Woodbury path keeps the reduced solve on the diagonal-Hessian route.

The reduced p × p system is solved densely, by Cholesky with an SVD fallback for a rank-deficient A, so this class scales in M (the number of variables) but not in p (the number of constraints). The subclasses BandedSchurSaddleSystem and ArrowSaddleSystem solve a banded or arrow-shaped S faster.

Parameters:

hessian (MatrixSystem) – The Hessian block H, of order M.
A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – The equality-constraint matrix, of shape (p, M), stored as CSR. A dense array is accepted and converted, but the point of this class is a sparse A. The constructor raises ValueError when its column count is not M.

property hessian#

Return the Hessian block H.

Returns:: hessian – The H block, of order M.
Return type:: MatrixSystem

property num_equality_constraints#

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

Returns:: p – The number of equality constraints.
Return type:: int

property A#

Return the sparse equality-constraint matrix A.

Returns:: A – The matrix A in scipy.sparse CSR form, of shape (p, M).
Return type:: Any

solve(b_top, b_bottom=None)[source]#

Solve the saddle-point system through the Schur complement on H.

Takes the same block-elimination steps as solve_kkt_system (Boyd and Vandenberghe, 2004, appendix C), but forms the reduced system with sparse products, so the dense p × M matrix is never built. The Hessian is touched only through its own solve.

Solves

 _        _     _         _     _          _
|  H  A^T  |   |  delta_x  |   |  b_top     |
|_ A   0  _|   |_   nu    _| = |_ b_bottom _|

for the primal step and the equality multipliers.

Parameters:

b_top (npt.NDArray[np.float64]) – Top block of the right-hand side (the negative gradient or dual residual): a vector of length M, or a matrix of shape (M, k) solved column by column.
b_bottom (npt.NDArray[np.float64] | None) – Bottom block (the equality residual), of length p. None means zero, the equality-feasible Newton step; a nonzero value arises in primal-dual methods, whose iterates need not satisfy A x = b mid-iteration.

Returns:

delta_x (npt.NDArray[np.float64]) – The primal step, the top block of the solution.
nu (npt.NDArray[np.float64]) – The equality multipliers, the bottom block of the solution.

Raises:

ValueError – When b_bottom does not have length p.
NewtonStepError – Propagated from the Hessian’s solve or from _reduced_solve.

References

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

class cvxium.linalg.systems.BandedSchurSaddleSystem(hessian, A)[source]#

Bases: SparseSaddleSystem

Sparse saddle-point system whose Schur complement is banded.

Forms the same complement S = A H^{-1} A^T as SparseSaddleSystem, but solves the reduced p × p system by banded Cholesky, O(p b^2) for bandwidth b, rather than dense O(p^3) Cholesky. This is the fast path when the equality rows are ordered so that S is banded: for example, when each row couples only to rows near it in a suitable sort order, the bandwidth stays small and the reduced solve is effectively linear in the number of constraints.

The caller asserts the structure by choosing this system (through BandedSaddleFactory). The bandwidth is read off the assembled sparse S, and the solve falls back to the dense path with SVD when S is too wide for banded storage to pay (2 (b + 1) > p), when H is not diagonal, or when S is not positive definite (a rank-deficient A), so it is always correct, if not always fast.

Parameters:

hessian (MatrixSystem) – The Hessian block H, of order M.
A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – The equality-constraint matrix, of shape (p, M). The constructor raises ValueError when its column count is not M.

class cvxium.linalg.systems.ArrowSaddleSystem(hessian, A, n_border)[source]#

Bases: BandedSchurSaddleSystem

Sparse saddle-point system whose reduced Schur complement is an arrow matrix.

When the equality rows split into a core, ordered so that the core complement A_core H^{-1} A_core^T is banded, plus n_border dense border rows (aggregate side constraints that touch many variables), the reduced complement S = A H^{-1} A^T is an arrow matrix: a banded core bordered by n_border dense rows and columns.

The reduced system is solved by wrapping S in an ArrowSystem with a BandedSystem core, so the core is factored once and the border is folded in through an n_border × n_border Schur complement, instead of a dense O(p^3) factorization of all of S. The reduced solve is thus exactly an ArrowSystem solve.

The border rows are the last n_border rows of A. With n_border = 0 this is the BandedSchurSaddleSystem solve. It falls back to the dense solve with SVD when the structure is absent (a non-diagonal Hessian, a core that is not narrowly banded, or a reduced system that is not positive definite), so the solution is the same either way.

Parameters:

hessian (MatrixSystem) – The Hessian block H, of order M.
A (npt.NDArray[np.float64] | scipy.sparse.spmatrix | scipy.sparse.sparray) – The equality-constraint matrix, of shape (p, M), whose last n_border rows form the dense border. The constructor raises ValueError when its column count is not M.
n_border (int) – Number of trailing rows of A that form the border. Must satisfy 0 ≤ n_border < p, or the constructor raises ValueError.

property num_border#

Return the number of trailing rows of A treated as the dense border.

Returns:: n_border – The border size given at construction.
Return type:: int

class cvxium.linalg.systems.SaddleFactory(*args, **kwargs)[source]#

Bases: Protocol

Protocol for building a SaddleSystem from a Hessian block and a matrix A.

A factory is how a caller chooses which structure of A to exploit (a dense Schur complement, a sparse one, a banded one, an arrow) without the equality layer branching on it: an EqualityConstraintLayer holds a factory and delegates to it, rather than being subclassed. The sparse saddle classes already have this call signature, so a factory can be one of the classes below or a configured object that carries extra structural metadata, as ArrowSaddleFactory carries its border size.

class cvxium.linalg.systems.SchurComplementSaddleFactory[source]#

Bases: object

Default saddle factory: a dense Schur-complement saddle-point system.

class cvxium.linalg.systems.SparseSaddleFactory[source]#

Bases: object

Saddle factory that exploits a sparse constraint matrix A.

class cvxium.linalg.systems.BandedSaddleFactory[source]#

Bases: object

Saddle factory that exploits a banded Schur complement.

class cvxium.linalg.systems.ArrowSaddleFactory(n_border)[source]#

Bases: object

Saddle factory for a reduced Schur complement shaped as a banded-core arrow.

Carries the border size, which the plain SaddleFactory call signature cannot convey. See ArrowSaddleSystem.

Parameters:: n_border (int) – Number of trailing rows of the constraint matrix that form the dense border.

Numerical Kernels#

Structured linear-system solvers and matrix-vector products for Cvxium.

Provides solve_* and multiply_* helpers that exploit the special structure of the Hessians arising in Cvxium’s interior-point Newton steps, so that those steps run in near-linear rather than cubic time. The supported structures are diagonal, rank-one and rank-p updates and downdates (through the Sherman-Morrison and Woodbury identities), arrow sparsity, block and block-plus-one Schur complements, symmetric banded, block diagonal, and the full KKT system (solve_kkt_system).

The helpers compose: a rank-p update can wrap an arrow solve, a block-diagonal solve can dispatch to a different structure per block, and so on. A structured solve takes the solve of its base matrix as a callable, so any helper here can serve as another’s base. Several solves also have a prepare/apply split (_prepare_rank_one_update, _prepare_rank_p_update, factor_banded) so that systems can factor once and reuse the factorization across the right-hand sides of one Newton step.

Solves raise NewtonStepError when the purported Hessian is not positive definite, is singular, or is non-finite, signaling callers to fall back rather than trust a corrupt Newton step. aspect_ratio_svd and tikhonov_filter support pseudo-inverses of rank-deficient constraint matrices.

cvxium.linalg.numerical_helpers.multiply_diagonal(y, eta)[source]#

Multiply a vector or matrix by a diagonal matrix.

Computes H y with H = diag(η), in linear time.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column.
eta (npt.NDArray[np.float64]) – Diagonal η of H, of length M.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When y is not 1D or 2D, or its leading dimension does not match eta.

cvxium.linalg.numerical_helpers.solve_diagonal(b, eta)[source]#

Solve a diagonal linear system.

Solves H x = b with H = diag(η), in linear time.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column.
eta (npt.NDArray[np.float64]) – Diagonal η of H, of length M.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

NewtonStepError – When some η_i ≤ 0, so H is not positive definite. Checked before the shapes.
ValueError – When b is not 1D or 2D, or its leading dimension does not match eta.

cvxium.linalg.numerical_helpers.solve_diagonal_eta_inverse(b, eta_inverse)[source]#

Solve a diagonal linear system given the reciprocal of its diagonal.

Solves H x = b with H = diag(η), given 1 / η, so the solve is a multiplication rather than a division. Linear time.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column.
eta_inverse (npt.NDArray[np.float64]) – The reciprocals 1 / η_i of the diagonal of H, of length M.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

NewtonStepError – When some 1 / η_i ≤ 0, so H is not positive definite. Checked before the shapes.
ValueError – When b is not 1D or 2D, or its leading dimension does not match eta_inverse.

cvxium.linalg.numerical_helpers.multiply_rank_one_update(y, kappa, A_multiply, d=None, **kwargs)[source]#

Multiply by a rank-one update of a structured matrix.

Computes H y with

H = A + d κ κ^T,

where A has structure that makes A z cheap to compute and κ is a vector. H is a rank-one update of A when d > 0 and a downdate when d < 0; when d is None, the weight is 1.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column.
kappa (npt.NDArray[np.float64]) – The vector κ, of length M.
A_multiply (Callable[..., npt.NDArray[np.float64]]) – Computes A z, called as A_multiply(z, **kwargs). Must accept a vector or a matrix and preserve its shape.
d (float | None) – Weight on the rank-one term; None means 1.
**kwargs (Any) – Extra arguments passed to A_multiply.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When y is not 1D or 2D, its leading dimension does not match kappa, or A_multiply does not preserve the shape of its input.

Notes

For a vector y, computing

H y = A y + κ (d κ^T y)

takes t + 4M flops, where t is the cost of A y:

Compute A y (t flops).
Compute the scalar d κ^T y (2M flops).
Add κ times that scalar to A y (2M flops).

cvxium.linalg.numerical_helpers.solve_rank_one_update(b, kappa, A_solve, d=None, **kwargs)[source]#

Solve a rank-one update of a structured system by Sherman-Morrison.

Solves H x = b with

H = A + d κ κ^T,

where A has structure that makes A y = c cheap to solve and κ is a vector. H is a rank-one update of A when d > 0 and a downdate when d < 0; when d is None, the weight is 1. Composes _prepare_rank_one_update and _apply_rank_one_update.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column.
kappa (npt.NDArray[np.float64]) – The vector κ, of length M.
A_solve (Callable[..., npt.NDArray[np.float64]]) – Solves A y = c, called as A_solve(c, **kwargs). Must accept a vector or a matrix and preserve its shape.
d (float | None) – Weight on the rank-one term; must be nonzero. None means 1.
**kwargs (Any) – Extra arguments passed to A_solve.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When d is zero, b is not 1D or 2D, its leading dimension does not match kappa, or A_solve does not preserve the shape of its input.
NewtonStepError – When the Sherman-Morrison scalar G is zero to within a relative tolerance (H singular), or has the wrong sign for a downdate (H indefinite). Propagated from _prepare_rank_one_update.

Notes

By the Sherman-Morrison formula, with the scalar

G = 1 / d + κ^T A^{-1} κ,

the solution is

x = A^{-1} b - ((κ^T A^{-1} b) / G) A^{-1} κ,

valid whenever G ≠ 0. When A is positive definite, the sign of G reveals the definiteness of H:

d > 0 (update): G > 0 always, and H is positive definite.
d < 0 (downdate): G < 0 means H is positive definite, G = 0 means H is singular, and G > 0 means H is indefinite (the downdate is too large).

The solve raises NewtonStepError unless H is positive definite.

For a vector b, the cost is 2t + 6M flops, where t is the cost of one solve with A:

Solve A ξ = κ (t flops).
Form G = 1 / d + κ^T ξ and its reciprocal (2M flops).
Solve A x’ = b (t flops).
Compute x = x’ - ((κ^T x’) / G) ξ (4M flops).

cvxium.linalg.numerical_helpers.multiply_rank_p_update(y, kappa, A_multiply, d=None, **kwargs)[source]#

Multiply by a rank-p update of a structured matrix.

Computes H y with

H = A + κ D κ^T,   D = diag(d),

where A has structure that makes A z cheap to compute and κ is M × p. Positive entries of d are updates and negative entries downdates; mixed signs are allowed. When d is None, D = I.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column.
kappa (npt.NDArray[np.float64]) – The matrix κ, of shape (M, p).
A_multiply (Callable[..., npt.NDArray[np.float64]]) – Computes A z, called as A_multiply(z, **kwargs). Must accept a vector or a matrix and preserve its shape.
d (npt.NDArray[np.float64] | None) – Diagonal of D, of length p; None means all ones.
**kwargs (Any) – Extra arguments passed to A_multiply.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When y is not 1D or 2D, its leading dimension does not match kappa, d does not have length p, or A_multiply does not preserve the shape of its input.

Notes

For a vector y, computing

H y = A y + κ (D (κ^T y))

takes t + 4Mp + p flops, where t is the cost of A y:

Compute A y (t flops).
Compute κ^T y (2Mp flops).
Scale by D (p flops).
Multiply by κ and add to A y (2Mp flops).

cvxium.linalg.numerical_helpers.solve_rank_p_update(b, kappa, A_solve, d=None, check_inertia=True, **kwargs)[source]#

Solve a rank-p update of a structured system by the Woodbury identity.

Solves H x = b with

H = A + κ D κ^T,   D = diag(d),

where A has structure that makes A y = c cheap to solve and κ is M × p with p ≪ M. Positive entries of d are updates and negative entries downdates; mixed signs are allowed. When d is None, D = I. Composes _prepare_rank_p_update and _apply_rank_p_update.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column.
kappa (npt.NDArray[np.float64]) – The matrix κ, of shape (M, p).
A_solve (Callable[..., npt.NDArray[np.float64]]) – Solves A y = c, called as A_solve(c, **kwargs). Must accept a vector or a matrix and preserve its shape.
d (npt.NDArray[np.float64] | None) – Diagonal of D, of length p, with no zero entries; None means all ones.
check_inertia (bool) – When True and d has mixed signs, verify that the inertia of the Woodbury matrix G matches the signs of d, which ensures H is positive definite. Set it to False when the caller detects a non-positive-definite Hessian by other means.
**kwargs (Any) – Extra arguments passed to A_solve.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When b is not 1D or 2D, kappa is not 2D, their leading dimensions disagree, d does not have length p or has a zero entry, or A_solve does not preserve the shape of its input.
NewtonStepError – When the Woodbury matrix G is singular, non-finite, or has inertia that does not match the signs of d, so H is not positive definite. Propagated from _prepare_rank_p_update and _apply_rank_p_update.

Notes

By the Woodbury identity, with the p × p Woodbury matrix

G = D^{-1} + κ^T A^{-1} κ,

the solution is

x = A^{-1} b - A^{-1} κ G^{-1} κ^T A^{-1} b,

valid whenever G is invertible. The code forms the equilibrated matrix

|D|^{1/2} G |D|^{1/2}

instead (see _prepare_rank_p_update), which has the same inertia. When A is positive definite, the signs of d determine how the positive definiteness of H is checked through G:

All d > 0 (pure update): G and H are always positive definite; G is factored by Cholesky.
All d < 0 (pure downdate): H is positive definite exactly when G is negative definite; -G is factored by Cholesky, and a failure means H is indefinite.
Mixed signs: H is positive definite exactly when G has as many positive eigenvalues as d has positive entries (see _verify_woodbury_inertia); G is factored by LU.

With q right-hand sides and t the cost of one solve with A, the cost is

(p + q) t + 2Mp (p + 2q)

plus O(p^3 + p^2 q) for G, which is negligible when p, q ≪ M:

Solve A ξ = κ, with ξ of shape (M, p) (p t).
Form G = D^{-1} + κ^T ξ (2Mp^2).
Factor G (O(p^3)).
Solve A x’ = b, with x’ of shape (M, q) (q t).
Compute z = κ^T x’, of shape (p, q) (2Mpq).
Solve G y = z (O(p^2 q)).
Compute x = x’ - ξ y (2Mpq).

cvxium.linalg.numerical_helpers.multiply_block_plus_one(y, A12, A22, A11_multiply, **kwargs)[source]#

Multiply by a matrix with a block-plus-one structure.

Computes H y for

      _              _
     |   A11    A12   |
H =  |                |,
     |_  A12^T  A22  _|

where A11 is M × M with structure that makes A11 z cheap to compute, A12 is a vector of length M, and A22 is a scalar.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column. Must have M + 1 rows.
A12 (npt.NDArray[np.float64]) – Last column of H above the diagonal, of length M.
A22 (float) – The bottom-right element of H.
A11_multiply (Callable[..., npt.NDArray[np.float64]]) – Computes A11 z, called as A11_multiply(z, **kwargs); must accept a vector or a matrix.
**kwargs (Any) – Extra arguments passed to A11_multiply.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When A12 is not 1D, y does not have M + 1 rows, or y is not 1D or 2D.

cvxium.linalg.numerical_helpers.solve_block_plus_one(b, A12, A22, A11_solve, **kwargs)[source]#

Solve a block-plus-one system by its scalar Schur complement.

Solves H x = b for

      _              _
     |   A11    A12   |
H =  |                |,
     |_  A12^T  A22  _|

where A11 is M × M with structure that makes A11 y = c cheap to solve, A12 is a vector of length M, and A22 is a scalar.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column. Must have M + 1 rows.
A12 (npt.NDArray[np.float64]) – Last column of H above the diagonal, of length M.
A22 (float) – The bottom-right element of H.
A11_solve (Callable[..., npt.NDArray[np.float64]]) – Solves A11 y = c, called as A11_solve(c, **kwargs); must accept a vector or a matrix.
**kwargs (Any) – Extra arguments passed to A11_solve.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When A12 is not 1D, b does not have M + 1 rows, or b is not 1D or 2D.
NewtonStepError – When the Schur complement s is not positive, so H is not positive definite. Errors from A11_solve propagate unchanged.

Notes

Let t be the cost of one solve with A11, and q the number of right-hand sides. Split b into its first M rows b1 and last row b2, and x likewise into x1 and x2.

Form A12’ = A11^{-1} A12 (t flops) and b1’ = A11^{-1} b1 (at most q t flops; a batched solve may be cheaper).
Form the scalar Schur complement
```
s = A22 - A12^T A12'
```
(2M flops). When H and A11 are positive definite, s > 0.
Compute x2 = (b2 - A12^T b1’) / s (2Mq flops and q divisions); x2 is a scalar or a vector of length q.
Compute x1 = b1’ - A12’ x2 (2Mq flops).

In total, that is (q + 1) t + 2M (2q + 1) + q flops.

cvxium.linalg.numerical_helpers.multiply_block(y, A12, A22, A11_multiply, **kwargs)[source]#

Multiply by a symmetric 2 × 2 block matrix with a structured leading block.

Computes H y for

      _              _
     |   A11    A12   |
H =  |                |,
     |_  A12^T  A22  _|

where A11 is M × M with structure that makes A11 z cheap to compute, A12 is M × p, and A22 is p × p, with p ≥ 2 (use multiply_block_plus_one when p = 1).

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column. Must have M + p rows.
A12 (npt.NDArray[np.float64]) – Off-diagonal block of H, of shape (M, p).
A22 (npt.NDArray[np.float64]) – Lower-right block of H, of shape (p, p).
A11_multiply (Callable[..., npt.NDArray[np.float64]]) – Computes A11 z, called as A11_multiply(z, **kwargs); must accept a vector or a matrix.
**kwargs (Any) – Extra arguments passed to A11_multiply.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When A12 is not a matrix with at least two columns, A22 is not p × p, y does not have M + p rows, or y is not 1D or 2D.

cvxium.linalg.numerical_helpers.multiply_arrow_sparsity_pattern(y, eta, zeta, theta)[source]#

Multiply by a matrix with an arrow sparsity pattern.

Computes H y for

      _               _
     |  diag(η)    ζ   |
H =  |                 |,
     |_   ζ^T      θ  _|

where η and ζ are vectors of length M and θ is a scalar, in linear time.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column. Must have M + 1 rows.
eta (npt.NDArray[np.float64]) – Diagonal η of the upper-left block of H, of length M.
zeta (npt.NDArray[np.float64]) – Last column ζ of H above the diagonal, of length M.
theta (float) – The bottom-right element θ of H.

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When eta and zeta differ in shape, y does not have M + 1 rows, or y is not 1D or 2D.

Notes

For a vector y = (y1, y2), with y1 of length M and y2 a scalar,

z[:M] = η ∘ y1 + ζ y2,
z[M]  = ζ^T y1 + θ y2,

which takes 3M + 1 multiplications and 2M additions, or 5M + 1 flops.

cvxium.linalg.numerical_helpers.solve_with_schur(b, A12, A22, A11_solve, **kwargs)[source]#

Solve a symmetric 2 × 2 block system by its Schur complement.

Solves H x = b for

      _              _
     |   A11    A12   |
H =  |                |,
     |_  A12^T  A22  _|

where A11 is M × M with structure that makes A11 y = c cheap to solve, A12 is M × p, and A22 is p × p, with p ≥ 2 (use solve_block_plus_one when p = 1).

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column. Must have M + p rows.
A12 (npt.NDArray[np.float64]) – Off-diagonal block of H, of shape (M, p).
A22 (npt.NDArray[np.float64]) – Lower-right block of H, of shape (p, p).
A11_solve (Callable[..., npt.NDArray[np.float64]]) – Solves A11 y = c, called as A11_solve(c, **kwargs); must accept a vector or a matrix.
**kwargs (Any) – Extra arguments passed to A11_solve.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When A12 is not a matrix with at least two columns, A22 is not p × p, b does not have M + p rows, or b is not 1D or 2D. Also propagated from scipy.linalg when the Schur complement is finite but the right-hand side is not.
NewtonStepError – When the Schur complement is not positive definite or is non-finite, so H is not positive definite. Errors from A11_solve propagate unchanged.

Notes

Let t be the cost of one solve with A11, and q the number of right-hand sides. Split b into its first M rows b1 and last p rows b2, and x likewise into x1 and x2.

Form A12’ = A11^{-1} A12 (at most p t flops; a batched solve may be cheaper) and b1’ = A11^{-1} b1 (at most q t flops).
Form the p × p Schur complement
```
S = A22 - A12^T A12'
```
(2Mp^2 flops), symmetrized to remove rounding skew. When H and A11 are positive definite, so is S.
Solve S x2 = b2 - A12^T b1’: forming the right-hand sides takes 2Mpq flops, the Cholesky factorization of S p^3 / 3, and the triangular solves 2p^2 q. x2 is p × q.
Compute x1 = b1’ - A12’ x2 (2Mpq flops).

In total, that is

(p + q) t + 2Mp (p + 2q) + p^3 / 3 + 2p^2 q

flops.

cvxium.linalg.numerical_helpers.solve_arrow_sparsity_pattern(b, eta, zeta, theta)[source]#

Solve a system with an arrow sparsity pattern in linear time.

Solves H x = b for

      _               _
     |  diag(η)    ζ   |
H =  |                 |,
     |_   ζ^T      θ  _|

where η and ζ are vectors of length M and θ is a scalar.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column. Must have M + 1 rows.
eta (npt.NDArray[np.float64]) – Diagonal η of the upper-left block of H, of length M.
zeta (npt.NDArray[np.float64]) – Last column ζ of H above the diagonal, of length M.
theta (float) – The bottom-right element θ of H.

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

NewtonStepError – When some η_i ≤ 0 (checked first, before the shapes), or when the Schur residual ψ^2 (see Notes) is not positive to within a relative tolerance of 10^-12, so H is not positive definite.
ValueError – When eta and zeta differ in shape, b does not have M + 1 rows, or b is not 1D or 2D.

Notes

The Cholesky factorization H = L L^T is available in closed form:

      _                               _
     |   diag(η)^{1/2}            0    |
L =  |                                 |
     |_  ζ^T diag(η)^{-1/2}       ψ   _|

with

ψ^2 = θ - ζ^T diag(η)^{-1} ζ.

Write b = (b1, b2) and x = (x1, x2), with b1 and x1 of length M. Forward substitution in L y = b gives

y1 = diag(η)^{-1/2} b1,   y2 = (b2 - ζ^T diag(η)^{-1} b1) / ψ,

and back substitution in L^T x = y gives

x2 = y2 / ψ = (b2 - (diag(η)^{-1} ζ)^T b1) / ψ^2,
x1 = diag(η)^{-1} b1 - x2 diag(η)^{-1} ζ.

So no square root is needed. Forming diag(η)^{-1} ζ takes M divisions and ψ^2 another M multiplications and M additions (summed with math.fsum, since θ and the inner product nearly cancel as the barrier term shrinks); x2 takes M multiplications, M additions, and one division; and x1 takes M divisions, M multiplications, and M additions. In total, for a vector b, that is 2M + 1 divisions, 3M multiplications, and 3M additions, or 8M + 1 flops.

cvxium.linalg.numerical_helpers.multiply_banded(y, ab, lower=False)[source]#

Multiply by a symmetric banded matrix in compact banded storage.

Computes H y for a symmetric banded H stored as for scipy.linalg.cholesky_banded.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column. Must have n rows.
ab (npt.NDArray[np.float64]) – Banded storage of the symmetric matrix H, of shape (p + 1, n), where n is the order of H and p its bandwidth (the number of superdiagonals, or of subdiagonals when lower). With lower False, ab[i, j] = H[j - p + i, j]; with lower True, ab[i, j] = H[j + i, j]. This is the layout of scipy.linalg.cholesky_banded.
lower (bool) – Whether ab stores the lower band (True) or the upper band (False).

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When y is not 1D or 2D, or its leading dimension does not match n.

Notes

Iterates over the p + 1 stored diagonals, using symmetry for the other triangle, so the product costs O(n p) per column rather than O(n^2) for a dense matrix.

cvxium.linalg.numerical_helpers.solve_banded(b, ab, lower=False)[source]#

Solve a symmetric positive definite banded system by banded Cholesky.

Solves H x = b for a symmetric positive definite banded H stored as for scipy.linalg.cholesky_banded. Composes factor_banded and solve_banded_factored.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column. Must have n rows.
ab (npt.NDArray[np.float64]) –
Banded storage of the symmetric matrix H, of shape (p + 1, n), where n is the order of H and p its bandwidth (the number of superdiagonals, or of subdiagonals when lower). With lower False, ab[i, j] = H[j - p + i, j]; with lower True, ab[i, j] = H[j + i, j]. This is the layout of scipy.linalg.cholesky_banded. To build the upper form from a dense H,
```
ab = np.zeros((p + 1, n))
for d in range(p + 1):
    ab[p - d, d:] = np.diag(H, d)
```
and the lower form,
```
ab = np.zeros((p + 1, n))
for d in range(p + 1):
    ab[d, : n - d] = np.diag(H, -d)
```
lower (bool) – Whether ab stores the lower band (True) or the upper band (False).

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When b is not 1D or 2D, or its leading dimension does not match n; also propagated from solve_banded_factored for a non-finite b.
NewtonStepError – Propagated from factor_banded when H is not positive definite or is non-finite.

Notes

For a matrix of order n and bandwidth p, the factorization costs O(n p^2) and each right-hand side O(n p), against O(n^3) and O(n^2) for a dense matrix.

cvxium.linalg.numerical_helpers.multiply_block_diagonal(y, block_sizes, block_multipliers)[source]#

Multiply by a block-diagonal matrix whose blocks have their own structure.

Computes H y for

      _                        _
     |  H_0                     |
     |      H_1                 |
H =  |            ...           |,
     |_                H_{k-1} _|

where each block H_i may have its own structure and block_multipliers[i] computes H_i z.

Parameters:

y (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the product is taken column by column.
block_sizes (Sequence[int]) – Order of each diagonal block; must sum to the number of rows of y.
block_multipliers (Sequence[Callable[..., npt.NDArray[np.float64]]]) – One callable per block; block_multipliers[i](z) computes H_i z. Each must have its block’s parameters already bound (for example, with functools.partial).

Returns:

z – The product H y, of the same shape as y.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When y is not 1D or 2D, block_sizes and block_multipliers differ in length, or block_sizes does not sum to the number of rows of y. Errors from the block callables propagate unchanged.

Notes

Because H is block diagonal,

H y = (H_0 y_0, H_1 y_1, ..., H_{k-1} y_{k-1}),

where y_i is the segment of y for block i. The blocks are multiplied independently, so the cost is the sum of the per-block costs, and each block may use a different structure (diagonal, arrow, rank-p update, and so on). For example, with block 0 diagonal and block 1 an arrow,

mult0 = functools.partial(multiply_diagonal, eta=eta0)
mult1 = functools.partial(
    multiply_arrow_sparsity_pattern, eta=eta1, zeta=zeta1, theta=theta1
)
z = multiply_block_diagonal(y, [M0, M1 + 1], [mult0, mult1])

cvxium.linalg.numerical_helpers.solve_block_diagonal(b, block_sizes, block_solvers)[source]#

Solve a block-diagonal system whose blocks have their own structure.

Solves H x = b for

      _                        _
     |  H_0                     |
     |      H_1                 |
H =  |            ...           |,
     |_                H_{k-1} _|

where each block H_i may have its own structure and block_solvers[i] solves H_i x_i = b_i.

Parameters:

b (npt.NDArray[np.float64]) – Right-hand side, a vector or a matrix; for a matrix, the system is solved for each column.
block_sizes (Sequence[int]) – Order of each diagonal block; must sum to the number of rows of b.
block_solvers (Sequence[Callable[..., npt.NDArray[np.float64]]]) – One callable per block; block_solvers[i](c) solves H_i x_i = c. Each must have its block’s parameters already bound (for example, with functools.partial).

Returns:

x – The solution, of the same shape as b.

Return type:

npt.NDArray[np.float64]

Raises:

ValueError – When b is not 1D or 2D, block_sizes and block_solvers differ in length, or block_sizes does not sum to the number of rows of b.
NewtonStepError – Propagated from a block solver when that block is not positive definite, is singular, or is non-finite.

Notes

Because H is block diagonal, H x = b splits into the k independent systems

H_i x_i = b_i,   i = 0, ..., k - 1,

where b_i is the segment of b for block i. The cost is the sum of the per-block costs, and each block may use a different structure (diagonal, arrow, rank-p update, and so on). For example, with block 0 diagonal and block 1 an arrow,

solver0 = functools.partial(solve_diagonal, eta=eta0)
solver1 = functools.partial(
    solve_arrow_sparsity_pattern, eta=eta1, zeta=zeta1, theta=theta1
)
x = solve_block_diagonal(b, [M0, M1 + 1], [solver0, solver1])

cvxium.linalg.numerical_helpers.solve_kkt_system(A, g, hessian_solve, g_bot=None, **kwargs)[source]#

Solve an equality-constrained KKT system by block elimination.

Solves

 _        _    _   _     _       _
|  H   A^T |  |  Δx |   |    g    |
|          |  |     | = |         |
|_ A    0 _|  |_ ν _|   |_ g_bot _|

where H is the Hessian, A the p × M equality-constraint matrix, and g_bot defaults to 0. When systems in H can be solved in O(M) time, the Schur complement gives Δx in O(p^3 + p^2 M) time.

Parameters:

A (npt.NDArray[np.float64]) – The equality-constraint matrix, of shape (p, M).
g (npt.NDArray[np.float64]) – Upper (gradient) block of the right-hand side, of length M.
hessian_solve (Callable[..., npt.NDArray[np.float64]]) – Solves H x = y, called as hessian_solve(y, **kwargs); must accept a vector or a matrix.
g_bot (npt.NDArray[np.float64] | None) – Lower (equality-residual) block of the right-hand side, of length p. None means zero, the equality-feasible Newton step. A nonzero g_bot arises in primal-dual methods, where the iterate need not satisfy A x = b mid-iteration.
**kwargs (Any) – Extra arguments passed to hessian_solve.

Returns:

delta_x (npt.NDArray[np.float64]) – The primal step Δx, of length M.
nu (npt.NDArray[np.float64]) – The dual variable ν, of length p.

Raises:

ValueError – When g does not have one entry per column of A, or g_bot is given without one entry per row of A.
NewtonStepError – Propagated from hessian_solve (H not positive definite, singular, or non-finite) or from _solve_schur_complement (Schur complement non-finite, or the KKT system has no solution).

Notes

Follows the block elimination of (Boyd and Vandenberghe, 2004, Algorithm C.4):

Form B = H^{-1} A^T and b = H^{-1} g: p + 1 solves with H, O((p + 1) M) time when each takes O(M).
Form -S = A B, which takes p^2 dot products of length M, O(p^2 M) time, and -c = A b - g_bot, O(p M) time.
Solve (-S) ν = -c by Cholesky in O(p^3) time (-S is positive definite when A has full row rank). When A is rank deficient, -S is too; the system can usually still be solved by SVD, which is slower, so Cholesky is tried first and SVD is the fallback (see _solve_schur_complement).
Solve H Δx = g - A^T ν: O(p M) time to form the right-hand side and O(M) for the solve.

In total that is O(p^2 M + p^3), dominated by forming S when p ≤ M.

References

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