Long-only efficient frontier benchmark challenge — Shrinkage as Preconditioning paper

minvar_frontier_challenge.py

"""
Long-Only Efficient Frontier Benchmark Challenge
=================================================
"Shrinkage as Preconditioning: Matrix-Free Methods for
 Long-Only Portfolio Optimization"

Generate the exact input data from the paper's efficient frontier experiment
(Section 7, Figure 3) and benchmark your solver against the reference timings.

Challenge: compute all 50 frontier portfolios faster than warm-started CG.

Problem
-------
For each rho in np.linspace(0, 2, 50), solve:

    min  w' Sigma_lw w  -  rho * mu' w
    s.t. sum(w) = 1,  w >= 0

where
    Sigma_lw = (1 - alpha) * X'X/T  +  alpha * bar_lam * I
    bar_lam  = ||X||_F^2 / (n * T)
    alpha    = 0.5

Reference hardware: Apple M4 Pro, 14-core CPU, 48 GB RAM
Reference software: Python 3.12, NumPy 2.4, SciPy 1.17

Reference timings (50-point sweep, best of 3 runs):
    Warm-start CG  (matrix-free, active-set):  0.055 s  (1.1 ms/point)
    Cold-start CG  (matrix-free, active-set):  0.234 s  (4.7 ms/point)
    CVXPY/Clarabel (naive loop):              ~135   s  (~2700 ms/point)

Verification checkpoints:
    rho=0.00  ->  vol=4.85% ann,  ret=5.98% ann,  130 active assets
    rho=1.00  ->  vol=6.40% ann,  ret=7.85% ann,   39 active assets
    rho=2.00  ->  vol=6.89% ann,  ret=7.89% ann,   30 active assets
"""

import numpy as np

# ---------------------------------------------------------------------------
# Reproduce simulate_equity_returns(n=500, T=1250, rng=42) exactly
# ---------------------------------------------------------------------------
N, T, K = 500, 1250, 50        # assets, time steps, factors (max(3, N//10))

rng = np.random.default_rng(42)

factor_vols = np.concatenate([[0.01], np.full(K - 1, 0.005)])
F = rng.standard_normal((T, K)) * factor_vols

B = np.zeros((N, K))
B[:, 0] = rng.uniform(0.4, 0.8, size=N)
for j in range(1, K):
    mask = rng.random(N) < 0.5
    B[mask, j] = rng.standard_normal(int(mask.sum())) * 0.2

idio_vols = rng.uniform(0.005, 0.015, size=N)
E = rng.standard_normal((T, N)) * idio_vols

X = F @ B.T + E
X -= X.mean(axis=0)          # shape (T, N) = (1250, 500), demeaned

# ---------------------------------------------------------------------------
# Expected returns  (same rng seed as paper's rng_ef = np.random.default_rng(42))
# ---------------------------------------------------------------------------
rng_mu = np.random.default_rng(42)
betas  = rng_mu.uniform(0.4, 0.8, N)
mu     = betas * (0.10 / 250)   # 10% annual market premium -> daily units

# ---------------------------------------------------------------------------
# LW shrinkage  (alpha = 0.5, scaled-identity target)
# ---------------------------------------------------------------------------
ALPHA   = 0.5
bar_lam = float(np.linalg.norm(X, "fro") ** 2) / (N * T)
target  = bar_lam * np.eye(N)
Sigma_lw = (1 - ALPHA) * (X.T @ X) / T + ALPHA * target

# ---------------------------------------------------------------------------
# Frontier parameter grid
# ---------------------------------------------------------------------------
rhos = np.linspace(0, 2, 50)   # 50 points, rho in [0, 2]

# ---------------------------------------------------------------------------
# Inputs summary
# ---------------------------------------------------------------------------
if __name__ == "__main__":
    kappa = np.linalg.cond(Sigma_lw)
    print(f"n={N}, T={T}, alpha={ALPHA}, kappa(Sigma_lw)={kappa:.1f}")
    print(f"bar_lam={bar_lam:.6f},  ||X||_F^2/(nT) check: {float(np.linalg.norm(X,'fro')**2/(N*T)):.6f}")
    print()
    print("Inputs:")
    print(f"  X          {X.shape}   demeaned return matrix")
    print(f"  mu         {mu.shape}          expected daily returns")
    print(f"  Sigma_lw   {Sigma_lw.shape}   shrunk covariance (N x N)")
    print(f"  rhos       {rhos.shape}         frontier grid")
    print(f"  ALPHA      {ALPHA}             shrinkage intensity")
    print()
    print("Reference timings (Apple M4 Pro, best of 3):")
    print("  Warm-start CG:  0.055 s for all 50 points")
    print("  Cold-start CG:  0.234 s for all 50 points")
    print("  CVXPY/Clarabel:  ~135 s for all 50 points")