fast-minimum-variance: Solving Minimum Variance Portfolios Fast

Overview

fast-minimum-variance solves the long-only minimum variance portfolio without ever forming the sample covariance matrix. The key observation is that the KKT stationarity condition $2\Sigma w = \lambda\mathbf{1}$ immediately gives $w \propto \Sigma^{-1}\mathbf{1}$: the entire problem reduces to one symmetric positive definite linear system $\Sigma v = \mathbf{1}$, solved matrix-free by conjugate gradients. The budget constraint is recovered by a single rescaling $w = v / (\mathbf{1}^\top v)$.

Working directly with the returns matrix $X \in \mathbb{R}^{T \times N}$ — rather than the assembled covariance $X^\top X$ — has two consequences. First, each conjugate gradient iteration costs $O(TN)$ rather than $O(N^2)$, and $X^\top X$ is never stored. Second, Ledoit-Wolf shrinkage enters as a simple row-augmentation of $X$: stacking $[\sqrt{1-\alpha},X;,\sqrt{\gamma},I]$ yields a matrix whose Gram matrix equals $\Sigma_{\text{LW}}$. The same CG code handles both the plain and shrunk problem without modification.

Quick Start

import numpy as np
from fast_minimum_variance import Problem

# 500 daily returns, 20 assets
X = np.random.default_rng(42).standard_normal((500, 20))

w, outer, inner = Problem(X).solve_cg()   # matrix-free CG — recommended

assert abs(w.sum() - 1.0) < 1e-8
assert (w >= 0).all()

Ledoit-Wolf Shrinkage

Ledoit-Wolf shrinkage plays a dual role: statistically it reduces estimation error; numerically it compresses the eigenvalue spectrum and directly cuts CG iteration counts. Use alpha = N / (N + T) as a simple analytical estimate of the optimal shrinkage intensity:

T, N = X.shape
w, outer, inner = Problem(X, alpha=N / (N + T)).solve_cg()

On S&P 500 equity data (495 assets, 1192 days), shrinkage cuts CG iterations from 685 to 205 and makes the matrix-free solver the fastest option by a wide margin.

Solvers

All solvers are methods on Problem and return (w, outer_steps, inner_iters) where $w \in \mathbb{R}^N$, $\sum_i w_i = 1$, $w_i \geq 0$.

Method	Approach	When to use
solve_cg()	Matrix-free conjugate gradients on the SPD reduced system	Default — fastest for large $N$, especially with shrinkage
solve_pcg()	Matrix-free PCG with an RMT (low-rank) preconditioner	Eigenvalue-cleaned shrinkage targets (requires pcg_lr)

solve_cg — matrix-free conjugate gradients

The inner step builds a LinearOperator that applies

$$v ;\mapsto; (1-\alpha),X_a^\top(X_a v) + \gamma v, \qquad \gamma = \frac{\alpha|X|_F^2}{N}$$

to a vector using two matrix-vector products with the active-asset submatrix $X_a$, without ever forming $\Sigma_a = X_a^\top X_a$. Standard CG then solves $\Sigma_a v = \mathbf{1}$. Ledoit-Wolf shrinkage ($\alpha > 0$) compresses the eigenvalue spectrum and reduces iteration counts dramatically — from nearly 2000 iterations at $\alpha \approx 0$ to single digits at $\alpha \approx 1$ in rank-deficient settings.

solve_pcg — preconditioned CG

Runs the same matrix-free active-set iteration as solve_cg, but preconditions the inner CG solve with the RMT low-rank target T0 (applied via the Woodbury identity at $O(n_a k)$ per step). Requires pcg_lr = (bar_lam, U_k, delta_k) from RMT preprocessing and returns the oracle-LW minimum-variance portfolio in far fewer iterations when the shrinkage target has been eigenvalue-cleaned.

Build pcg_lr either from the dense rmt_target_and_alpha (full eigh) or, for large $N$, from rmt_preconditioner_rsvd, which recovers the top eigenpairs via a randomized SVD of $X$ — matrix-free at $O(TNk)$, never forming $X^\top X$. Because a preconditioner only affects the iteration count and not the solution, the randomized factors match the dense ones in CG iterations while cutting setup cost by up to ~10× at $N \sim 2000$:

from fast_minimum_variance.shrinkage.util import rmt_preconditioner_rsvd

pcg_lr = rmt_preconditioner_rsvd(X, n_components=16)
w, outer, inner = Problem(X, alpha=N / (N + T), pcg_lr=pcg_lr).solve_pcg()

The Primal-Dual Active-Set Loop

Long-only weights are enforced by an outer loop that wraps any inner solver:

1. Primal step. Solve the budget-only equality system over the current active asset set. Drop any asset with weight below $-\varepsilon$ (multiple assets at once if violations are large).
2. Dual step. Once all active weights are non-negative, compute the gradient $\nabla_i f(w) = 2[(1-\alpha)(X^\top X w)_i + \gamma w_i] - \rho\mu_i$ for every excluded asset. If any excluded asset has $\nabla_i f(w) < \lambda$ (the budget multiplier), it would decrease variance if added — re-insert the most-violated asset and repeat.
3. Termination. The loop exits when primal and dual feasibility hold simultaneously. Combined with stationarity from the inner solve, this is sufficient for global optimality.

With Ledoit-Wolf shrinkage at the analytically optimal $\alpha$, the loop typically converges in 2–4 outer iterations on real equity data.

Problem Variants

The same solver handles a range of portfolio construction problems by choosing $\alpha$, $\rho$, $\mu$:

Problem	alpha	rho	mu
Minimum variance	$0$	$0$	—
Mean-variance (Markowitz)	any	$> 0$	expected returns
Minimum tracking error to benchmark $b$	any	$2$	X.T @ (X @ b)
LW-regularised minimum variance	$N/(N+T)$	$0$	—

# Mean-variance
mu = np.random.default_rng(0).standard_normal(N)  # expected returns, shape (N,)
w, *_ = Problem(X, rho=1.0, mu=mu).solve_cg()

# Minimum tracking error to benchmark b
b = np.ones(N) / N  # equal-weight benchmark
mu_te = X.T @ (X @ b)
w, *_ = Problem(X, rho=2.0, mu=mu_te).solve_cg()

When rho != 0, two SPD solves are performed per outer step: $\Sigma_a v_1 = \mathbf{1}$ and $\Sigma_a v_2 = \mu_a$. The budget multiplier $\lambda$ is recovered analytically from the budget constraint, avoiding the full saddle-point system.

Balance Systems

To replace the default budget constraint $\mathbf{1}^\top w = 1$ with a general set of linear equality constraints $B w = c$ (e.g. sleeve budgets, factor-exposure targets), pass a balance system (B, c):

B = np.zeros((2, N)); B[0, :N // 2] = 1.0; B[1, N // 2:] = 1.0  # each half holds...
c = np.array([0.5, 0.5])                                        # ...half of the budget
w, *_ = Problem(X, B=B, c=c).solve_cg()

Long-only ($w \ge 0$) is still enforced. B must have full row rank on every active set the shrinking loop visits. Use this path only when you need it — the default path (no B, c) is faster for the standard budget + long-only problem.

Benchmarks

All timings on Apple M4 Pro, Python 3.12, NumPy 2.4, SciPy 1.17.

Universe	$N$	$T$	solve_cg time (s)
Synthetic i.i.d. Gaussian	1000	2000	0.019
S&P 500 (Jul 2021–Apr 2026)	495	1192	0.0091

Both with Ledoit-Wolf shrinkage ($\alpha = 0.333$ synthetic / $0.293$ S&P), 56 and 205 CG iterations respectively.

Installation

pip install fast-minimum-variance

For development:

git clone https://github.com/Jebel-Quant/fast_minimum_variance
cd fast_minimum_variance
make install

Requirements

• Python 3.11+
• numpy
• scipy
• scikit-learn
• cvx-linalg

Citing

If you use this library in academic work or research, please cite:

@software{fast_minimum_variance,
  author  = {Schmelzer, Thomas},
  title   = {fast-minimum-variance: Solving Minimum Variance Portfolios Fast},
  url     = {https://github.com/Jebel-Quant/fast_minimum_variance},
  year    = {2026},
  license = {MIT}
}

License

MIT License — see LICENSE for details.