cosa¶

A conic active-set algorithm for second-order cone programs, applied to mean–standard-deviation portfolio optimization.

Portfolio optimization is usually written as a quadratic program in variance. COSA writes it in standard deviation instead —

maximize   mu' x - lam * sqrt(x' Sigma x)      subject to   A x <= b,  E x = d

— which makes the objective linear and pushes the risk into a second-order cone. That is a worse fit for an interior-point method and a better one for an active set: the risk term is now a single geometric object whose activity is a fact about the solution, not a term to be traded off, and a nearby problem usually has a nearby active set.

Whether that pays off is the research question. The short answer, measured, is at the bottom.

Install¶

uv add cosa                      # or: pip install cosa
uv add "cosa[reference]"         # adds CVXPY + Clarabel, for cross-checking answers

Solve a portfolio¶

import numpy as np
from cosa import solve_portfolio

rng = np.random.default_rng(0)
factors = rng.normal(size=(8, 4))
covariance = factors @ factors.T + 0.05 * np.eye(8)
returns = rng.normal(0.08, 0.03, 8)

answer = solve_portfolio(returns, covariance, lam=2.0, long_only=True)

print(f"return {answer.expected_return:.4f}  risk {answer.risk:.4f}  optimal {answer.is_optimal}")

return 0.0925  risk 0.0962  optimal True

answer.holdings is the vector of weights, summing to one. answer.active names the constraints binding at the solution — Success Criterion 3's "interpretable in terms of the active portfolio constraints". answer.residuals carries the five conic KKT residuals that certify it.

solve_portfolio adds sum(x) == 1 by default and x >= 0 on request. Anything else goes in as matrices:

sectors = np.array([[1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]])  # first three assets
capped = solve_portfolio(
    returns,
    covariance,
    lam=2.0,
    long_only=True,
    inequalities=(sectors, np.array([0.25])),  # at most 25% in that sector
)
assert sectors @ capped.holdings <= 0.25 + 1e-8

It raises rather than returning a position it does not stand behind. A solve that does not reach a certified optimum raises NotOptimalError, which carries the result and the residuals. Pass strict=False to get the answer and judge it yourself.

Trace a frontier¶

This is what the method is for. Each solve is handed the last one's answer, working set, multipliers and factorization:

from cosa.linear_algebra.reuse import Reuse

cache, previous, work = Reuse(), None, 0
for lam in np.geomspace(1.0, 6.0, 12):
    point = solve_portfolio(returns, covariance, lam=lam, long_only=True, warm=previous, cache=cache)
    work += point.metrics.iterations
    previous = point.warm()  # the solution, working set, multipliers and factorizations

On that sweep warm starting cuts the work by more than half — 932 iterations cold, 378 warm. The condition under which it does is in the results.

Reproduce every number¶

uv run python -m cosa.experiments

Writes four artifacts to docs/experiments/: the failure-mode study, the frontier sweep, the four-mode benchmark, and the environment that produced them. Seeds are arguments with recorded defaults, so the command with no flags reproduces the committed files byte-for-byte — except the two that report wall-clock time and a platform, which cannot be and say so.

Does it work?¶

Measured, not asserted. Every number below is regenerated by the command above.

Warm starting pays exactly when the active set transfers. This is the finding, and it is a conditional rather than a claim:

Correcting a belief costs more than acquiring one: a cold solve discovers the active set on the way in, a warm one has to undo a wrong answer first and then discover it anyway. The sign of the total is decided by the mix.

Factorization reuse works, and does not yet pay in wall clock. The share of solves needing a fresh factorization falls from 98.9% to 1.4% — a whole solve typically factorizes twice. But a column update against a full (n,n) orthogonal factor costs O(n²) while a fresh QR of an (n,m) matrix costs O(n m²), so updating wins only once m² exceeds n. The classical argument assumes a working set comparable in size to the problem; an active-set method on a portfolio runs with m well below n. Measured: 0.98× at n = 300, 1.12× at n = 500.

Accuracy holds on twelve of thirteen instance families — including the ill-conditioned, the rank-deficient and the degenerate — and across all four benchmark modes, agreeing with CVXPY/Clarabel to within 1e-6.

The thirteenth is the interesting one. On an instance whose constraint matrix spans fourteen orders of magnitude, COSA terminates reporting optimal with all five conic KKT residuals under 1e-11, at a point whose objective is 3.4% worse than the reference's — and the reference's point is feasible for COSA's own check to 1e-11. The residual is not lying: it is 1.9e-5 absolute, and §14.2 divides it by |c| ≈ 2e6. A relative KKT residual cannot certify an instance like that, and no amount of scaling fixes it. The study reports this as a distinct verdict — wrong, meaning certified and disagreeing — because a certificate that certifies the wrong answer is worse than no certificate.

Wall clock does not. The reference solver is 3–30× faster on every mode measured, its modelling overhead included. That is reported rather than buried: the paper's stated goal is a characterization of when conic active-set methods work well, which is not a claim that they always do.

Documentation¶

• Architecture — package layout, and the decisions recorded once: why H is not rho*I, what reuse is worth.
• Sign convention — the one thing four modules must agree about.
• Failure modes — what COSA does on the hard instances, and which mitigation is why.
• How the bugs were actually found — the mechanisms rather than the bugs, and why the ones that compared against something external found everything the ones that asked the code about itself did not.
• The plan — the project's own specification, which every module docstring cites by section.

Development¶

uvx rhiza-task all      # format, lint, types, tests, coverage, docs coverage

The test suite is the documentation of intent: every module has a paired test file, and every claim in a docstring that could be wrong has a test that would fail if it were.