Hierarchical Risk Parity (HRP)

Hierarchical Risk Parity (HRP), introduced by Marcos López de Prado, builds diversified allocations without ever inverting the covariance matrix. That single design choice is the point: the instability of mean-variance optimization comes almost entirely from Σ⁻¹, which amplifies the least reliable parts of an estimated covariance into enormous, concentrated bets. HRP replaces the inversion with clustering and recursive risk allocation, so it captures correlation structure while staying numerically robust. This article develops the conditioning argument quantitatively, walks the three steps with a full implementation, and is honest about where HRP helps and where it does not.

Why covariance inversion is the problem

Markowitz weights are w ∝ Σ⁻¹μ. The trouble is Σ⁻¹. When assets are correlated the sample covariance is near-singular: its smallest eigenvalues are tiny and biased downward, and inverting divides by those tiny numbers, magnifying estimation noise. The relevant diagnostic is the condition number — the ratio of largest to smallest eigenvalue — which grows with the number of correlated assets. Marchenko-Pastur theory makes this precise: for N assets and T observations with ratio q = N/T, the sample eigenvalue spectrum spreads over a width that grows with q, so the smallest eigenvalues are pushed toward zero exactly as you add the correlated assets you wanted for diversification.

import numpy as np

def condition_number(cov):
    eig = np.linalg.eigvalsh(cov)
    return eig.max() / eig.min()

The symptoms are familiar: extreme long-short weights, allocations that flip when one correlation is nudged, and out-of-sample performance that tracks the noise in the sample rather than the future. Shrinkage and constraints help, but HRP takes a more structural view: mean-variance treats all assets as flat, fully substitutable competitors, ignoring the hierarchy that real markets obviously have (stocks cluster by sector, bonds by duration). HRP exploits that hierarchy instead of fighting it.

Step 1: tree clustering

Convert the correlation matrix to a distance metric and run hierarchical agglomerative clustering. The standard correlation distance is:

d_ij = sqrt( 0.5 * (1 - rho_ij) )

so perfectly correlated assets are distance 0 and perfectly anti-correlated assets are distance 1. Clustering on this distance builds a dendrogram — a tree where similar assets merge first. This recovers the discrete grouping structure that PCA expresses as continuous factors, but in an interpretable hierarchy.

Step 2: quasi-diagonalization

Reorder rows and columns of the covariance matrix according to the dendrogram leaf order (seriation), so that similar assets sit adjacent. After this, large covariances cluster around the diagonal and the matrix becomes "quasi-diagonal." This reordering is what lets the next step allocate between genuinely related groups rather than between arbitrary, index-ordered assets.

Step 3: recursive bisection

Walk down the tree, splitting each cluster into two and allocating capital between the halves inversely to their risk. Each half's risk is computed with an inverse-variance allocation within the cluster, and the parent weight is split so the riskier branch gets less. Because allocation happens between clusters at every level, correlated assets share a risk budget instead of competing for it.

Full implementation

import numpy as np
import pandas as pd
from scipy.cluster.hierarchy import linkage
from scipy.spatial.distance import squareform

def correl_distance(corr):
    return np.sqrt(0.5 * (1 - corr))

def quasi_diag(link):
    link = link.astype(int)
    sort_ix = pd.Series([link[-1, 0], link[-1, 1]])
    num_items = link[-1, 3]
    while sort_ix.max() >= num_items:
        sort_ix.index = range(0, sort_ix.shape[0] * 2, 2)
        df0 = sort_ix[sort_ix >= num_items]
        i, j = df0.index, df0.values - num_items
        sort_ix[i] = link[j, 0]
        df1 = pd.Series(link[j, 1], index=i + 1)
        sort_ix = pd.concat([sort_ix, df1]).sort_index()
        sort_ix.index = range(sort_ix.shape[0])
    return sort_ix.tolist()

def cluster_var(cov, items):
    sub = cov.loc[items, items]
    ivp = 1.0 / np.diag(sub)             # inverse-variance weights
    ivp /= ivp.sum()
    return float(ivp @ sub.values @ ivp)

def recursive_bisection(cov, sort_ix):
    w = pd.Series(1.0, index=sort_ix)
    clusters = [sort_ix]
    while clusters:
        clusters = [c[j:k] for c in clusters
                    for j, k in ((0, len(c) // 2), (len(c) // 2, len(c)))
                    if len(c) > 1]
        for i in range(0, len(clusters), 2):
            c0, c1 = clusters[i], clusters[i + 1]
            v0, v1 = cluster_var(cov, c0), cluster_var(cov, c1)
            alpha = 1 - v0 / (v0 + v1)    # less weight to the riskier branch
            w[c0] *= alpha
            w[c1] *= 1 - alpha
    return w

def hrp(returns: pd.DataFrame) -> pd.Series:
    cov, corr = returns.cov(), returns.corr()
    dist = squareform(correl_distance(corr).values, checks=False)
    link = linkage(dist, method="single")
    order = [corr.index[i] for i in quasi_diag(link)]
    return recursive_bisection(cov, order)

The output is a fully invested, long-only portfolio. The heaviest linear-algebra operation anywhere is inverting a diagonal of variances — numerically trivial and stable. There is no Σ⁻¹, so the conditioning problem that destabilizes Markowitz simply does not arise.

Robustifying the inputs

HRP is more robust than mean-variance but still consumes an estimated covariance, so the same estimation hygiene applies. With N assets you want well over N observations; otherwise denoise the matrix first. Two improvements are worth building in:

from sklearn.covariance import LedoitWolf

def shrunk_cov_corr(returns):
    cov = LedoitWolf().fit(returns.values).covariance_
    cov = pd.DataFrame(cov, index=returns.columns, columns=returns.columns)
    d = np.sqrt(np.diag(cov))
    corr = cov / np.outer(d, d)
    return cov, corr

Ledoit-Wolf shrinkage stabilizes the within-cluster variance estimates, and an eigenvalue-clipping ("constant residual") denoising of the correlation matrix — also from López de Prado — further removes the noisy small-eigenvalue directions before clustering. The linkage method matters too: single linkage (the default above) can "chain" elongated clusters, so test ward and average on your own universe rather than assuming the default is best.

Denoising the correlation matrix first

The single most impactful refinement to HRP is to denoise the correlation matrix before clustering. The empirical eigenvalue spectrum of a sample correlation matrix contains a bulk of noise eigenvalues that, under the Marchenko-Pastur null, carry no information. The constant-residual-eigenvalue method fits the MP distribution to locate the noise band, then replaces all eigenvalues inside it with their average while keeping the signal eigenvalues intact. Clustering on the denoised matrix produces far more stable, meaningful trees.

import numpy as np

def denoise_corr(corr, q, bandwidth=0.25):
    """q = T / N. Replace noise eigenvalues with their average (constant residual)."""
    evals, evecs = np.linalg.eigh(corr)
    lam_plus = (1 + np.sqrt(1 / q)) ** 2          # MP upper edge (var=1)
    noise = evals < lam_plus
    if noise.any():
        evals = evals.copy()
        evals[noise] = evals[noise].mean()        # flatten the noise band
    out = evecs @ np.diag(evals) @ evecs.T
    d = np.sqrt(np.diag(out))
    return out / np.outer(d, d)                   # renormalize to a correlation matrix

The intuition mirrors why HRP works at all: the small, noisy eigenvalues are precisely what mean-variance inversion amplifies into nonsense, and flattening them removes the noise from the structure HRP then clusters on. In practice, denoise-then-HRP is more robust than either alone, and it is a cleaner way to handle the T not-much-larger-than N regime than hoping single-linkage clustering ignores the noise.

HRP vs Markowitz vs naive risk parity

Method	Uses Σ⁻¹?	Uses correlations?	Stability	Concentration risk
Mean-variance (Markowitz)	Yes	Yes	Poor	High
Naive risk parity (inverse-vol)	No	No	Good	Moderate
Hierarchical risk parity	No	Yes (via tree)	Good	Low

Naive risk parity is stable but discards correlation information, so it over-allocates to clusters of redundant assets. Markowitz uses correlations but pays with instability. HRP is the middle path: it uses the correlation structure through clustering without the dangerous inversion. In López de Prado's Monte Carlo studies, HRP delivered lower out-of-sample variance than both alternatives despite never being the in-sample variance minimizer — the signature of a method that generalizes rather than overfits.

Variants worth knowing

The base algorithm has several refinements that address its known weak points. Hierarchical equal risk contribution (HERC) replaces the somewhat arbitrary inverse-variance split at each node with a true equal-risk-contribution allocation within and between clusters, and it uses the dendrogram's natural cut (via a gap statistic) to decide how many clusters to form rather than bisecting blindly. Nested clustered optimization (NCO), also from López de Prado, goes further: it clusters assets, runs a small, well-conditioned mean-variance optimization within each cluster and across clusters, and stitches the results together — combining HRP's robustness with the ability to express return views inside each cluster.

The unifying idea across all these variants is to confine any matrix inversion to small, well-conditioned sub-problems where it is safe, and to use the hierarchy to keep the between-cluster allocation stable. That is the general lesson HRP teaches even beyond the specific algorithm: the danger is not optimization per se, but optimizing a large, ill-conditioned covariance all at once. Break the problem along its natural structure and much of the instability disappears.

Honest limits

HRP is a robustness tool, not a free lunch, and a few caveats deserve emphasis:

• It allocates risk, not return. HRP ignores expected returns entirely. If you have

genuine alpha forecasts, HRP alone leaves them on the table — blend it with a return view (for example via Black-Litterman).

• It is not variance-optimal in sample, and an attacker constructing the worst case

can beat it; its advantage is expected out-of-sample robustness, not a guarantee.

• Clusters are non-stationary. Correlation regimes shift, so re-estimate (and

re-cluster) periodically rather than freezing the tree, and expect the diversification to degrade when correlations spike in crises.

• Turnover and costs are real. HRP weights drift and re-cluster over time; rebalance

with a turnover band and explicit transaction costs rather than chasing every small change in the tree.

• Linkage and distance choices change the answer, so they should be validated like

any other hyperparameter, not assumed.

The core message is narrow and durable: HRP removes the single most damaging step in classical optimization — covariance inversion — and replaces it with a stable, structured risk allocation that needs only scipy clustering and a few dozen lines. It is one of the most practical answers to mean-variance fragility, best deployed with a shrunk, denoised covariance, a sensible rebalancing band, and clear awareness that it manages risk structure, not return.