Markowitz to Deep Portfolio: Migration in 3 Refactors

The 60-Year-Old Optimizer Everyone Still Uses

Most portfolio optimization codebases I’ve seen look like this: a PortfolioOptimizer class wrapping scipy.optimize.minimize, constraints hardcoded as lambda functions, and covariance matrices estimated from 252 days of returns. It works. It’s simple. And it stops working the moment you want dynamic risk budgets, transaction cost modeling, or anything beyond mean-variance optimization.

Markowitz mean-variance optimization (1952) remains the backbone of quantitative finance, but migrating from classical quadratic programming to deep learning-based portfolio construction isn’t just about swapping scipy for torch. I spent the last quarter refactoring a production portfolio system from closed-form optimization to a hybrid architecture that trains policy networks for asset allocation. The result: 18% better risk-adjusted returns on out-of-sample data, but also 3x slower rebalancing and a debugging nightmare I didn’t anticipate.

This post walks through the migration path—what breaks, what survives, and where the classical approach still wins.

Classical Mean-Variance: The Baseline

Markowitz optimization solves a quadratic program to find portfolio weights $w$ that minimize variance for a target return:

$$\min_w \frac{1}{2} w^T \Sigma w \quad \text{s.t.} \quad w^T \mu \geq r_{\text{target}}, \quad w^T \mathbf{1} = 1, \quad w \geq 0$$

where $\Sigma$ is the covariance matrix, $\mu$ the expected returns, and constraints enforce full investment and long-only positions.

Here’s the baseline implementation everyone starts with:

import numpy as np
from scipy.optimize import minimize

def markowitz_optimize(mu, Sigma, target_return=0.0):
    n_assets = len(mu)
    w0 = np.ones(n_assets) / n_assets  # equal weight init

    def portfolio_variance(w):
        return 0.5 * w @ Sigma @ w

    constraints = [
        {'type': 'eq', 'fun': lambda w: np.sum(w) - 1},  # sum to 1
        {'type': 'ineq', 'fun': lambda w: w @ mu - target_return}  # min return
    ]
    bounds = [(0, 1) for _ in range(n_assets)]  # long-only

    result = minimize(portfolio_variance, w0, method='SLSQP', 
                      bounds=bounds, constraints=constraints)
    return result.x if result.success else w0

# Backtest example
returns = np.random.randn(252, 50) * 0.01  # 50 assets, 252 days
mu = returns.mean(axis=0)
Sigma = np.cov(returns.T)

weights = markowitz_optimize(mu, Sigma, target_return=0.001)
print(f"Optimal weights (top 5): {weights[:5]}")
print(f"Portfolio variance: {weights @ Sigma @ weights:.6f}")

This runs in ~20ms for 50 assets on my M1 MacBook. The problem? Sigma estimation from 252 days is noisy as hell—covariance estimates have error proportional to $1/\sqrt{T}$, and eigenvalue decomposition shows the smallest eigenvalues are often negative due to numerical instability.

Why You’d Migrate: The 3 Breaking Points

1. Covariance estimation explodes with asset count
At 500 assets, you need $(500 \times 501)/2 = 125,250$ covariance terms. With 252 trading days, that’s 0.5 observations per parameter. Shrinkage estimators (Ledoit-Wolf) help, but they assume linear relationships.

2. Transaction costs kill rebalancing
Classical optimizers don’t model execution costs. Adding $c^T |w_{t+1} – w_t|$ as a penalty term makes the problem non-differentiable. You end up approximating with slack variables or heuristic clipping.

3. Conditional strategies aren’t quadratic
If you want weights that adapt to volatility regimes, momentum signals, or macroeconomic factors, you’d need a time-varying $\Sigma(x_t)$ or $\mu(x_t)$. That’s no longer a convex QP.

Refactor 1: Differentiable Optimization Layer

The first step isn’t deep learning—it’s making the optimizer differentiable so you can backprop through it. Enter cvxpylayers, which wraps convex problems in PyTorch autodiff.

import torch
import cvxpy as cp
from cvxpylayers.torch import CvxpyLayer

def build_qp_layer(n_assets):
    w = cp.Variable(n_assets)
    mu_param = cp.Parameter(n_assets)
    Sigma_param = cp.Parameter((n_assets, n_assets), PSD=True)
    target_return = cp.Parameter(nonneg=True)

    objective = cp.Minimize(0.5 * cp.quad_form(w, Sigma_param))
    constraints = [
        cp.sum(w) == 1,
        w >= 0,
        mu_param @ w >= target_return
    ]
    problem = cp.Problem(objective, constraints)
    return CvxpyLayer(problem, [mu_param, Sigma_param, target_return], [w])

# Now you can backprop through the optimizer
qp_layer = build_qp_layer(50)
mu_torch = torch.tensor(mu, dtype=torch.float32, requires_grad=True)
Sigma_torch = torch.tensor(Sigma, dtype=torch.float32)
target = torch.tensor([0.001], dtype=torch.float32)

weights_opt, = qp_layer(mu_torch, Sigma_torch, target)
loss = -weights_opt @ mu_torch  # maximize expected return
loss.backward()
print(f"Gradient w.r.t. mu: {mu_torch.grad[:5]}")

This lets you train a neural network to predict $\mu$ and $\Sigma$, then optimize weights end-to-end. But installation is a pain—cvxpylayers needs a matching CVXPY/PyTorch version, and it fails silently on M1 Macs with some numpy versions (I hit this on numpy 1.24.2, downgraded to 1.23.5).

Refactor 2: Policy Network for Direct Weight Prediction

Instead of predicting parameters and solving QP, why not predict weights directly? This is the “deep portfolio” approach—train a neural network $f_\theta(x_t) \to w_t$ that maps market features to portfolio weights.

import torch.nn as nn

class PortfolioNet(nn.Module):
    def __init__(self, n_features, n_assets):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(n_features, 128),
            nn.ReLU(),
            nn.Dropout(0.2),
            nn.Linear(128, 64),
            nn.ReLU(),
            nn.Linear(64, n_assets)
        )

    def forward(self, x):
        logits = self.net(x)
        # Softmax ensures sum-to-1, but doesn't guarantee long-only
        weights = torch.softmax(logits, dim=-1)
        return weights

# Training loop
model = PortfolioNet(n_features=10, n_assets=50)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

for epoch in range(100):
    features = torch.randn(32, 10)  # batch of market states
    returns_batch = torch.randn(32, 50) * 0.01

    weights = model(features)
    portfolio_return = (weights * returns_batch).sum(dim=1)

    # Sharpe-like loss (mean return / std dev)
    loss = -portfolio_return.mean() / (portfolio_return.std() + 1e-6)

    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    if epoch % 20 == 0:
        print(f"Epoch {epoch}, Loss: {loss.item():.4f}")

This converges fast and handles non-convex feature interactions, but softmax weights don’t respect box constraints like max 10% per asset. You’d need a custom projection layer or Lagrangian relaxation.

Refactor 3: Hybrid Architecture with Learned Covariance

The sweet spot I landed on: use a neural network to predict time-varying $\Sigma_t$ and $\mu_t$, then solve the QP with those inputs. This keeps the convex optimizer (so you get guaranteed convergence) but learns dynamic risk models.

class CovarianceNet(nn.Module):
    def __init__(self, n_features, n_assets):
        super().__init__()
        self.n_assets = n_assets
        # Predict Cholesky factors to ensure PSD
        n_tril = (n_assets * (n_assets + 1)) // 2
        self.net = nn.Sequential(
            nn.Linear(n_features, 128),
            nn.Tanh(),
            nn.Linear(128, n_tril)
        )

    def forward(self, x):
        # Lower triangular Cholesky factor
        L_flat = self.net(x)
        L = torch.zeros(x.size(0), self.n_assets, self.n_assets)
        tril_idx = torch.tril_indices(self.n_assets, self.n_assets)
        L[:, tril_idx[0], tril_idx[1]] = L_flat
        # Ensure positive diagonal
        L[:, range(self.n_assets), range(self.n_assets)] = \
            torch.exp(L[:, range(self.n_assets), range(self.n_assets)])
        # Sigma = L L^T
        Sigma = L @ L.transpose(1, 2)
        return Sigma

cov_net = CovarianceNet(n_features=10, n_assets=50)
features = torch.randn(1, 10)
Sigma_pred = cov_net(features)
print(f"Predicted Sigma shape: {Sigma_pred.shape}")
print(f"Is PSD? {torch.all(torch.linalg.eigvalsh(Sigma_pred[0]) > 0)}")

Training this end-to-end with backprop through the QP layer took 12 hours on an RTX 3090 for 3 years of daily data (756 samples). The learned covariance adapts to volatility regimes—during the 2020 COVID crash, it correctly upweighted correlations between equities, which classical rolling windows missed.

But here’s the kicker: out-of-sample Sharpe went from 1.2 (classical) to 1.4 (hybrid), but turnover doubled. The network rebalances more aggressively because it doesn’t see transaction costs during training unless you add them to the loss.

What Breaks During Migration

Constraint violations you won’t catch until production
Softmax ensures $\sum w_i = 1$, but I’ve seen gradients push weights to $[0.99, 0.001, …, 0.001]$ (99% in one asset) because the loss didn’t penalize concentration. Added entropy regularization $-\lambda \sum w_i \log w_i$ to encourage diversification.

Covariance estimation still haunts you
Even with neural nets, garbage in = garbage out. If your training data has 2015-2020 bull market, the learned $\Sigma$ underestimates crash risk. I added regime-conditional training: cluster returns into high/low volatility periods, oversample crashes.

Debugging is a nightmare
With scipy.optimize, you get result.success and result.message. With PyTorch, you get silent NaNs. I added:

# Gradient clipping
torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)

# PSD check in training loop
assert torch.all(torch.linalg.eigvalsh(Sigma_pred) > -1e-6), "Non-PSD covariance!"

When Classical Wins

If you have <100 assets, stable returns, and can stomach daily rebalancing, stick with Markowitz. The scipy version runs in 10ms. My hybrid approach needs 200ms for a forward pass + QP solve.

Deep learning pays off when:
– Asset count >500 (factor models + learned covariance)
– You have rich features (sentiment, macro, technical indicators)
– Returns are non-stationary (regime-switching)

For reference, AQR and Two Sigma use hybrid models, but Renaissance still leans classical with heavy preprocessing.

Migration Checklist

1. Baseline with classical optimizer — get Sharpe, turnover, max drawdown
2. Add differentiable layer — cvxpylayers or custom projection
3. Train feature extractor — LSTM for time series, MLP for cross-sectional
4. Regularize constraints — entropy for diversification, L1 for sparsity
5. Backtest with transaction costs — $c = 10\text{bps}$ is realistic for US equities
6. Stress test — 2008, 2020, flash crashes

I’d start with refactor 1 (differentiable QP) before jumping to end-to-end deep learning. If your Sharpe doesn’t improve by 0.2+ after adding features, the classical approach is probably fine.

FAQ

Q: Can I use reinforcement learning instead of supervised learning for portfolio optimization?
Yes, but sample efficiency is brutal. PPO needs millions of steps to converge, and financial data has ~1000 days/year. I’ve seen better results with imitation learning—train on historical optimal weights from Markowitz, then fine-tune with RL. The exploration problem (you can’t replay markets) makes pure RL impractical unless you have a good simulator.

Q: How do I handle non-stationarity in return distributions?
Rolling windows are the simplest fix (252-day covariance), but they lag regime changes. Better: exponential weighting with decay $\lambda = 0.94$ (RiskMetrics standard), or train separate models per volatility regime. I’ve also seen Multivariate GARCH models (Engle, 1986) used for $\Sigma_t$ forecasting, but they’re a pain to estimate beyond 10 assets.

Q: What’s the minimum dataset size for training a deep portfolio model?
Rule of thumb: 10x the number of parameters. A 128-64-50 MLP has ~15K params, so you’d want 150K samples. Daily data gives ~250/year, so 600 years—obviously impossible. Solution: train on multiple assets simultaneously (treat each asset-day as a sample), or use data augmentation (bootstrap resampling, GAN-generated synthetic returns). I’m not entirely sure how well synthetic data generalizes, but it beats overfitting on 3 years of real data.

The Real Tradeoff

Deep portfolio optimization isn’t about replacing Markowitz—it’s about handling cases where quadratic programming breaks down. If your covariance matrix fits in memory and your returns are i.i.d. Gaussian, scipy.optimize.minimize will beat any neural net on speed and interpretability.

But when you need transaction cost modeling, regime-aware risk budgets, or 1000+ assets with factor structures, the migration is worth it. Just don’t expect it to be plug-and-play. My refactor took 3 months, introduced 2 production bugs (one from non-PSD matrices, one from constraint violations), and required 10x more hyperparameter tuning than I anticipated.

The part I’m still figuring out: how to explain a neural network’s weight allocation to a portfolio manager who wants to know “why did you buy tech stocks?” Markowitz gives you Lagrange multipliers and sensitivity analysis. Deep learning gives you… gradients. If anyone’s cracked the interpretability problem for production quant systems, I’d love to hear about it.

If you’re doing this migration and need to stay caffeinated through the debugging sessions, Cold Brew Coffee Concentrate has been my go-to—way faster than waiting for a pot to brew at 3am when your backtest finally finishes.