SDP_weeds

`SDP.py`: `solve_sdp`

Source: src/SDP.py · User guide: Semidefinite Programs · Theory: Scenario Theory

solve_sdp solves a scenario program whose scenario constraints are linear matrix inequalities (LMIs). The objective can be linear or convex quadratic. As in solve_lp and solve_qp, it then computes the complexity $k$ and a degeneracy flag with get_support. The interface differs from the LP/QP solvers in several ways:

matrices are passed as dicts keyed by variable index;
n is read from the length of c;
x is 1-D, and x_ref is flattened to match;
Q is optional.

Problem solved

$$\begin{aligned} \min_{x \in \mathbb{R}^{n},\ \zeta \in \mathbb{R}^{N}} \quad & \tfrac{1}{2}\, x^\top Q\, x \;+\; c^\top x \;+\; \tau \lVert x - \bar{x} \rVert_p \;+\; \rho \sum_{i=1}^{N} \zeta_i \\\ \text{s.t.} \quad & F_0(\delta_i) + \sum_{j=1}^{n} x_j\, F_j(\delta_i) \;\preceq\; \zeta_i\, I_m, \qquad i = 1, \dots, N, \\\ & E_0 + \sum_{j=1}^{n} x_j\, E_j \;\preceq\; 0, \\\ & \zeta_i \;\ge\; 0, \qquad i = 1, \dots, N, \end{aligned}$$

where $F_j(\delta) \in \mathbb{S}^{m}$ and $E_j \in \mathbb{S}^{m_E}$ are symmetric, and $\preceq$ is the semidefinite order. A single scalar slack $\zeta_i$ shifts every eigenvalue of scenario $i$'s LMI. So at an optimum with $\rho > 0$:

$$\zeta_i^\ast = \max\big(0,\ \lambda_{\max}\big(F(x^\ast, \delta_i)\big)\big),$$

which was confirmed numerically on the example problem below, solved with $\rho = 0.2$. The formulations come from tau, rho and include_slack as in the LP mapping table, and the regularization term is added only when $\tau > 0$.

Signature

from src.SDP import solve_sdp

x_out, zeta_out, cost_out, N, complexity, constraints, degeneracy = solve_sdp(
    deltas, F_d, E, c, Q,
    tau=0.0, x_ref=np.array([0.0]), rho=0.0, norm_type=2,
    solver=None, include_slack=None,
)

Parameters

Name	Type / shape	Default	Meaning
`deltas`	`numpy.ndarray`, `(N, q)` or `(N,)`	required	The scenarios, one per row. Each row goes unchanged to `F_d`. An empty `deltas` raises `ValueError`.
`F_d`	callable, `δ → dict (m, m)`	required	Returns the LMI coefficients for one scenario, as `{'0': F0, '1': F1, ..., 'n': Fn}` (integer keys also work). See Key convention. Every matrix must be exactly symmetric. `m` is read from the first matrix in `F_d(deltas[0])`.
`E`	`dict (m_E, m_E)` or `None`	required	The hard LMI $E_0 + \sum_j x_j E_j \preceq 0$, with the same key convention. `None` or `{}` means no hard LMI. `m_E` can differ from `m`. Every matrix must be exactly symmetric.
`c`	`numpy.ndarray`, `(n,)` or `(n, 1)`	required	Linear cost vector. It sets `n`, the number of entries of `c`. A row vector `(1, n)` fails with `ValueError: Incompatible dimensions`.
`Q`	`numpy.ndarray` `(n, n)`, or `None`	required (may be `None`)	Quadratic cost matrix. `None` or an empty array means no quadratic term (a zero matrix is used). Otherwise it is checked as in `solve_qp`: symmetric and positive semidefinite up to a small tolerance. A non-empty 1-D array raises `LinAlgError`.
`tau`	`float`	`0.0`	Regularization weight $\tau$.
`x_ref`	scalar, or `n` entries (`(n,)` or `(n, 1)`)	`np.array([0.0])`	Regularization centre $\bar{x}$. It is flattened: one entry is used for every coordinate, `n` entries give a length-`n` vector, and any other size raises `ValueError`.
`rho`	`float`	`0.0`	Slack penalty $\rho$.
`norm_type`	any number `p >= 1`, `"inf"` / `np.inf`, `"fro"`	`2`	Order of the vector norm of `x - x_ref`, as for LP/QP; `"fro"` is the 2-norm. Used only when `tau > 0`. (`"nuc"` is not a vector norm and raises `RuntimeError`.)
`solver`	`str` or `None`	`None`	CVXPY solver. It must support PSD cones (see Solver compatibility).
`include_slack`	`bool` or `None`	`None`	`True` forces slack, `False` forbids it, and `None` means slack exactly when `rho != 0`.

Key convention for `F_d` and `E`

The keys are '0', '1', …, 'n', as strings or as integers (0, 1, …, n); each key k is read with int(k). 0 is the constant term. j is the coefficient of $x_j$, which is the code's x[int(j) - 1], so variables are numbered from 1.
Order does not matter, and keys can be left out (a missing key means a zero coefficient). Keys are read from the dict returned for each scenario, so different scenarios may use different key sets.
A key larger than n raises IndexError when x is indexed.
generate_tensor_function and generate_numeric_F both produce dicts with string keys.

Returns

Name	Type / shape	Meaning
`x_out`	`numpy.ndarray`, `(n,)`	Optimal $x$. This is 1-D, unlike the `(n, 1)` columns returned by `solve_lp` and `solve_qp`.
`zeta_out`	`numpy.ndarray`, `(N,)`	Optimal slacks, or `np.zeros(N)` without slack.
`cost_out`	`float`	`prob.value`, the full objective including the $\tau$ and $\rho$ terms. The web interface shows $c^\top x + \tfrac12 x^\top Q x$ instead.
`N`	`int`	`len(deltas)`.
`complexity`	`int`	$k$, the number of scenario LMIs in the support list.
`constraints`	`list` of `N` `cvxpy.constraints.PSD`	One PSD constraint per scenario, in the order of `deltas`. The hard LMI is excluded. Variable values and duals are overwritten by later re-solves, as described in the LP caveat.
`degeneracy`	`bool`	`True` when the lower risk bound is not certified.

Exceptions

Exception	When
`AssertionError`	`Q` is not symmetric, or has a negative eigenvalue, beyond the tolerance of `solve_qp`.
`AssertionError`	A matrix from `F_d(δ)` is not exactly symmetric: `"$F_j(\delta)$ need to be symmetric"`. A matrix in `E` is not symmetric: `"$F_j$ needs to be symmetric"`. The messages contain MathJax delimiters (`$…$`) for display in the web interface.
`numpy.linalg.LinAlgError`	`Q` is a non-empty 1-D array.
`IndexError`	An `F_d`/`E` key is larger than `n`.
`ValueError`	`deltas` is empty, `x_ref` has neither 1 nor `n` entries, `c` has an incompatible shape, or the status is not `optimal` or `optimal_inaccurate` (`"SDP did not solve to optimality. ..."`).
`cvxpy.error.SolverError`	The solver is missing or has no PSD support.
`ValueError`	`get_support` cannot recover a support list.

How it works

Read the dimensions and validate Q. n is the number of entries of c. An empty Q (or None) is replaced by np.zeros((n, n)); otherwise Q = check_psd(Q) runs the symmetry and PSD checks of QP. Then N = check_scenarios(deltas) and m = list(F_d(deltas[0]).values())[0].shape[0].
Variables. x = cp.Variable(n) is 1-D. zeta is cp.Variable(N, nonneg=True) with slack, or np.zeros(N) without.
Scenario LMIs. For each scenario i, the code calls F_d(deltas[i]) once. It asserts that each matrix is symmetric, and adds up the terms:
```
expr = F['0'] + x[0] * F['1'] + ... + x[n-1] * F['n']      # missing keys are skipped
constraints.append(expr << zeta[i] * np.eye(m))
```
In CVXPY, A << B creates a PSD constraint on B - A. So constraints[i].args[0] is the affine matrix expression $\zeta_i I_m - F(x, \delta_i)$ (checked numerically). One constraint object per scenario means $k$ counts scenarios.
Hard LMI. If E is truthy (not None, not {}), the same sum over E gives non_risk_constraints = [expr << 0]. Otherwise the list is empty. The hard LMI is always enforced and is never a support-list candidate.
Objective.
```
quadratic = (1/2) * cp.quad_form(x, cp.psd_wrap(Q)) if np.any(Q) else 0
objective = cp.Minimize(quadratic + c.T @ x
                        + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))
```
Before this, x_ref is flattened to a scalar or a length-n vector. As in the QP, cp.psd_wrap tells CVXPY to trust the check already made on Q, and the norm term is left out when tau == 0. The quadratic term is also left out when $Q = 0$: CVXPY would still add a cone for it, which solvers without second-order cones (e.g. SDPA) cannot handle.
Solve and check the status. prob.solve(solver=solver) runs with default settings. The status must be optimal or optimal_inaccurate.
Support list and degeneracy. get_support(constraints, non_risk_constraints, prob, objective, solver=solver) works exactly as for the LP (see Misc_weeds). Only the way an LMI is screened and tested differs, as described next.

How support and violation of an LMI scenario are measured

Step	What is computed for scenario $i$	Rule
Dual screen	`dual_value` is an `m × m` PSD matrix $Z_i$	candidate if $\max_{jk} \lvert (Z_i)_{jk} \rvert > 10^{-8}$
Violated under relaxation	stationarity in $\zeta_i$ gives $\rho - \operatorname{tr}(Z_i) - \mu_i = 0$, and $\mu_i = 0$ when $\zeta_i > 0$	$\operatorname{tr}(Z_i) = \rho$. For the example problem with `rho=0.2`, every violated scenario had $\operatorname{tr}(Z_i) = 0.2$.
Activity margin (`_is_slack`)	$\text{margin} = \lambda_{\min}\big(\tfrac12 (S + S^\top)\big)$ with $S = \zeta_i I - F(x^\ast, \delta_i)$, and $\text{scale} = \max \lvert S_{jk} \rvert$	strictly slack (a spurious screen) if $\text{margin} > 10^{-6} + 10^{-6}\cdot\text{scale}$
Support membership	re-solve without the scenario	of support if the optimal value changes by more than `rtol=1e-6` (`atol=1e-9`)

The margin is recorded before any re-solve, because re-solves overwrite x.value.

Parallelism

There is none. F_d is evaluated and the LMIs are built in a sequential loop, and all get_support re-solves run in sequence.

Solver compatibility

The problem always contains PSD cones. According to CVXPY 1.7.2's solver metadata, the solvers that accept PSD constraints are CLARABEL, SCS, MOSEK, CVXOPT, SDPA and COPT. get_solvers offers only these for SDPs (COPT through CVXPY's metadata). Choosing any other solver in Python ends in SolverError. SDPA has no second-order cones, so it needs $Q = 0$ and $\tau = 0$ (or $p = 1$ or inf).

Accuracy matters here too: on the example below, SCS returns the right $\lambda$ but $k = 5$ with degeneracy=True, where CLARABEL returns $k = 1$. MOSEK is the recommended solver.

The MOSEK path

In the web interface, uploading a MOSEK licence runs the solve in a separate process, src/mosek_solve.py. That process runs the same solve_form as every other solve, with the solver set to "MOSEK". solve_sdp itself has no MOSEK-specific code. MOSEK is simply passed as solver, which you can also do from Python if MOSEK and a licence are installed.

How the web interface calls it

solve_form in src/solve_form.py calls

solve_sdp(scenarios, F_d, E, c, Q, tau, theta_bar, rho, p, run_solver, include_slack=include_slack)

In symbolic mode, F_d comes from generate_tensor_function. In numeric mode it comes from generate_numeric_F.
E comes from generate_tensor, or is {} when the field is empty.
Q comes from generate_matrix, or is an empty array (no quadratic term) when the field is empty.
theta_bar (that is, $\bar{x}$) is parsed into a 1-D array.

Example

This example computes a scenario upper bound on the largest eigenvalue of a random symmetric matrix. The single decision variable is $\lambda$, and each scenario LMI is $M(\delta_i) - \lambda I \preceq 0$. Run it from the repository root.

import numpy as np
from src.SDP import solve_sdp
from src.Risk import quantify_risk
from src.parsing import generate_tensor_function

# Smallest lambda with M(delta_i) - lambda*I <= 0 (in the PSD order) for every sample,
# i.e. a scenario upper bound on the largest eigenvalue of a random symmetric matrix.
F_d = generate_tensor_function(
    '{"0": [["delta[0]", "delta[1]"], ["delta[1]", "delta[2]"]],'
    ' "1": [["-1", "0"], ["0", "-1"]]}'
)
rng = np.random.default_rng(2)
deltas = rng.standard_normal((100, 3))                # N = 100 scenarios, q = 3

c = np.array([1.0])                                    # (n,), n = 1 (x = lambda)
Q = None                                               # no quadratic term

x, zeta, cost, N, k, cons, degenerate = solve_sdp(deltas, F_d, None, c, Q, solver="CLARABEL")
lam_max = max(np.linalg.eigvalsh(F_d(d)["0"]).max() for d in deltas)
print("x =", x, " shape", x.shape, " k =", k, " degenerate =", degenerate)
print("largest sampled eigenvalue =", lam_max)
print("constraint type:", type(cons[0]).__name__, " dual shape:", cons[0].dual_value.shape)
epsL, epsU = quantify_risk(k, N, 1e-6)
print(f"risk in [{epsL:.4f}, {epsU:.4f}] with confidence 1 - 1e-6")

# Hard LMI E_0 + x_1 E_1 <= 0, i.e. lambda <= 2.5 (binds), plus relaxation so the problem stays feasible
E = {"0": -2.5 * np.eye(2), "1": np.eye(2)}
x, zeta, cost, N, k, cons, degenerate = solve_sdp(deltas, F_d, E, c, Q, rho=0.2, solver="CLARABEL")
print("relaxed + hard LMI: x =", x, " k =", k, " violated =", int((zeta > 1e-7).sum()), " degenerate =", degenerate)

Output (CVXPY 1.7.2, CLARABEL):

x = [4.22844486]  shape (1,)  k = 1  degenerate = False
largest sampled eigenvalue = 4.228444854277633
constraint type: PSD  dual shape: (2, 2)
risk in [0.0000, 0.1858] with confidence 1 - 1e-6
relaxed + hard LMI: x = [2.5]  k = 7  violated = 7  degenerate = False

What the output shows:

Robust run. The single scenario with the largest $\lambda_{\max}$ is the support list, so $k = 1$.
Relaxed run. The hard LMI caps $\lambda$ at 2.5. The 7 scenarios whose $\lambda_{\max}$ exceeds 2.5 are violated, and they make up the support list. The hard LMI is not counted in $k$.

Navigation Menu

Repository navigation

SDP_weeds

`SDP.py`: `solve_sdp`

Problem solved

Signature

Parameters

Key convention for `F_d` and `E`

Returns

Exceptions

How it works

How support and violation of an LMI scenario are measured

Parallelism

Solver compatibility

The MOSEK path

How the web interface calls it

Example

See also

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Using The Tool

Starting Up

Convex Optimization

Graphical Results

Getting into the Weeds

Scenario Theory

Implementation Details

Clone this wiki locally

SDP_weeds

SDP.py: solve_sdp

Problem solved

Signature

Parameters

Key convention for F_d and E

Returns

Exceptions

How it works

How support and violation of an LMI scenario are measured

Parallelism

Solver compatibility

The MOSEK path

How the web interface calls it

Example

See also

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Using The Tool

Starting Up

Convex Optimization

Graphical Results

Getting into the Weeds

Scenario Theory

Implementation Details

Clone this wiki locally

`SDP.py`: `solve_sdp`

Key convention for `F_d` and `E`