"""Schur Complementary Allocation estimator."""
# Copyright (c) 2023-2026
# Author: Hugo Delatte <hugo.delatte@skfoliolabs.com>
# Credits: Peter Cotto, Michal Kaszubski.
# SPDX-License-Identifier: BSD-3-Clause
# Implementation derived from:
# Precise, Copyright (c) 2021, Peter Cotton.
# Algorithm considerations:
# We apply Tail Call Optimization (TCO), replacing recursion with iteration
# and updating the covariance in place, to reduce call-stack depth and memory usage.
# To ensure portfolio variance decreases monotonically with the regularization
# factor gamma, we identify the variance turning point and cap gamma at its maximum
# permissible value. See https://github.com/skfolio/skfolio/discussions/3
from __future__ import annotations
import math
from collections.abc import Callable
from typing import Any
import numpy as np
import skfolio.typing as skt
from skfolio.cluster import HierarchicalClustering
from skfolio.distance import BaseDistance
from skfolio.exceptions import OptimizationError
from skfolio.optimization.hierarchical._seriation._base import (
_BaseSeriatedOptimization,
)
from skfolio.optimization.hierarchical._utils import (
_apply_weight_constraints_to_split_factor,
)
from skfolio.prior import BasePrior, ReturnDistribution
from skfolio.seriation import BaseSeriation
from skfolio.typing import ArrayLike, FloatArray, IntArray
from skfolio.utils.stats import (
cov_nearest,
inverse_multiply,
is_cholesky_dec,
multiply_by_inverse,
symmetric_step_up_matrix,
symmetrize,
)
from skfolio.utils.tools import bisection
[docs]
class SchurComplementary(_BaseSeriatedOptimization):
r"""Schur Complementary Allocation estimator.
Schur Complementary Allocation is a portfolio allocation method developed by Peter
Cotton [1]_.
It uses Schur-complement-inspired augmentation of sub-covariance matrices,
revealing a link between Hierarchical Risk Parity (HRP) and minimum-variance
portfolios (MVP).
By tuning the regularization factor `gamma`, which governs how much off-diagonal
information is incorporated into the augmented covariance blocks, the method
smoothly interpolates from the heuristic divide-and-conquer allocation of HRP
(`gamma = 0`) to the MVP solution (`gamma -> 1`).
The algorithm begins by computing a distance matrix and performing hierarchical
clustering, then applies seriation to reorder assets in the dendrogram so that
adjacent leaves have minimal distance.
Next, it uses recursive bisection: starting with the top-level cluster, each cluster
is split into two sub-clusters in a top-down traversal.
For each sub-cluster, an augmented covariance matrix is built based on the Schur
complement to incorporate off-diagonal block information. From this matrix, the
total cluster variance under an inverse-variance allocation is computed, and a
weighting factor derived from the variances of the two sub-clusters is used to
update their cluster weights.
Notes
-----
A poorly conditioned covariance matrix can prevent convergence to the MVP solution
as gamma approaches one. Setting `keep_monotonic=True` (the default) ensures that
the portfolio variance decreases monotonically with respect to gamma and remains
bounded by the variance of the HRP portfolio (`variance(Schur) <= variance(HRP)`),
even in the presence of ill-conditioned covariance matrices. Additionally, you can
apply shrinkage or other conditioning techniques via the `prior_estimator` parameter
to improve numerical stability and estimation accuracy.
Parameters
----------
gamma : float
Regularization factor in [0, 1].
When gamma is zero, no off-diagonal information is used (equivalent to HRP).
As gamma approaches one, the allocation moves toward the minimum variance
solution. The better the conditioning of the initial covariance matrix, the
closer the allocation will get to the MVP solution when gamma is near one.
keep_monotonic : bool, default=True
If True, ensures that the portfolio variance decreases monotonically with
respect to gamma. This is achieved by capping gamma at its maximum permissible
value (`effective_gamma_`). This constraint guarantees that the solution remains
variance-bounded by the HRP portfolio (`variance(Schur) <= variance(HRP)`),
even in the presence of ill-conditioned covariance matrices.
If False, no monotonicity enforcement or gamma capping is applied.
For more details, see: https://github.com/skfolio/skfolio/discussions/3
prior_estimator : BasePrior, optional
:ref:`Prior estimator <prior>`.
The prior estimator is used to estimate the :class:`~skfolio.prior.ReturnDistribution`
containing estimates of expected asset returns, covariance matrix and
returns. The moments and returns estimations are used for the risk computation
and the returns estimation are used by the distance matrix estimator.
The default (`None`) is to use :class:`~skfolio.prior.EmpiricalPrior`.
distance_estimator : BaseDistance, optional
:ref:`Distance estimator <distance>`.
The distance estimator is used to estimate the codependence and the distance
matrix used by the seriation estimator.
The default (`None`) is to use :class:`~skfolio.distance.PearsonDistance`.
hierarchical_clustering_estimator : HierarchicalClustering, optional
Deprecated. Use `seriation_estimator=HierarchicalSeriation(
hierarchical_clustering_estimator=...)`. Supplying both parameters is an error.
Will be removed in version 2.0.
seriation_estimator : BaseSeriation, optional
Asset ordering estimator. The default is
:class:`~skfolio.seriation.HierarchicalSeriation` with Ward linkage
and optimal leaf ordering. Use :class:`~skfolio.seriation.SpectralSeriation`
for spectral ordering with orientation carried across online updates.
distance_from_prior : bool, default=True
If True, fit the distance estimator on the return scenarios produced by
the prior. With `CovarianceDistance(covariance_estimator="precomputed")`,
use the prior's covariance instead.
If False, use the `X` argument supplied to this optimizer's `fit(X)` or
`partial_fit(X)`, before the prior processes it.
During online learning, True refits the distance estimator on the prior's
current scenarios or covariance at each update. False updates it from
the new observations in `X` and requires support for `partial_fit`.
Portfolio allocation always uses the prior's moments and return scenarios.
See :ref:`asset_seriation`.
min_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
Minimum assets weights (weights lower bounds). The default is 0.0 (no short
selling). Negative weights are not allowed. If a float is provided, it is
applied to each asset. `None` is equivalent to the default `0.0`. If a
dictionary is provided, its (key/value) pair must be the (asset name/asset
minimum weight) and the input `X` of the `fit` methods must be a DataFrame with
the asset names in columns. When using a dictionary, assets values that are not
provided are assigned the default minimum weight of `0.0`.
Example:
* `min_weights = 0.0` --> long only portfolio (default).
* `min_weights = {"SX5E": 0.1, "SPX": 0.2}`
* `min_weights = [0.1, 0.2]`
max_weights : float | dict[str, float] | array-like of shape (n_assets, ), default=1.0
Maximum assets weights (weights upper bounds). The default is 1.0 (each asset
is below 100%). Weights above 1.0 are not allowed. If a float is provided, it is
applied to each asset. `None` is equivalent to the default `1.0`. If a
dictionary is provided, its (key/value) pair must be the (asset name/asset
maximum weight) and the input `X` of the `fit` method must be a DataFrame with
the asset names in columns. When using a dictionary, assets values that are not
provided are assigned the default maximum weight of `1.0`.
Example:
* `max_weights = 1.0` --> each weight must be below 100% (default).
* `max_weights = 0.5` --> each weight must be below 50%.
* `max_weights = {"SX5E": 0.8, "SPX": 0.9}`
* `max_weights = [0.8, 0.9]`
transaction_costs : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
Transaction costs of the assets.
If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset cost) and the input `X` of the `fit` method must be a
DataFrame with the asset names in columns.
The default value is `0.0`.
management_fees : float | dict[str, float] | array-like of shape (n_assets, ), default=0.0
Management fees of the assets.
If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset fee) and the input `X` of the `fit` method must be a
DataFrame with the asset names in columns.
The default value is `0.0`.
previous_weights : float | dict[str, float] | array-like of shape (n_assets, ), optional
Previous weights of the assets. Previous weights are used to compute the
portfolio total cost. If a float is provided, it is applied to each asset.
If a dictionary is provided, its (key/value) pair must be the
(asset name/asset previous weight) and the input `X` of the `fit` method must
be a DataFrame with the asset names in columns.
The default (`None`) means no previous weights.
Additionally, when `fallback="previous_weights"`, failures will fall back to
these weights if provided.
portfolio_params : dict, optional
Portfolio parameters forwarded to the resulting `Portfolio` in `predict`.
If not provided and if available on the estimator, the following
attributes are propagated to the portfolio by default: `name`,
`transaction_costs`, `management_fees`, `previous_weights` and `risk_free_rate`.
fallback : BaseOptimization | "previous_weights" | list[BaseOptimization | "previous_weights"], optional
Fallback estimator or a list of estimators to try, in order, when the primary
optimization raises during `fit`. Alternatively, use `"previous_weights"` (alone
or in a list) to fall back to the estimator's `previous_weights`. When a
fallback succeeds, its fitted `weights_` are copied back to the primary
estimator so that `fit` still returns the original instance. For traceability,
`fallback_` stores the successful estimator (or the string `"previous_weights"`)
and `fallback_chain_` stores each attempt with the associated outcome. With
`partial_fit`, only None or `"previous_weights"` is supported because fallback
estimators have not accumulated the primary model's online history. See
:ref:`optimization_fallbacks`.
raise_on_failure : bool, default=True
Controls error handling when fitting fails and no fallback succeeds. If True,
the estimator raises the final error. If False, the estimator emits a warning
and sets `weights_` to None, so subsequent calls to `predict` return a
:class:`~skfolio.portfolio.FailedPortfolio`. During `fit`, `raise_on_failure`
applies to any fitting error, including errors raised by the prior estimator.
During `partial_fit`, `raise_on_failure` applies only to optimization failures
after learning completes. Input validation failures and errors from the prior or
other learning estimators are always raised. See
:ref:`optimization_failure_handling` for batch recovery and
:ref:`online_failure_handling` for online continuation and restart rules.
Attributes
----------
weights_ : ndarray of shape (n_assets,)
Weights of the assets.
effective_gamma_ : float or None
If `keep_monotonic` is True, the highest permissible `gamma` that preserves
monotonic variance decrease. Otherwise, it equals the input `gamma`.
It is None after an optimization failure or when a fallback supplies the
allocation.
distance_estimator_ : BaseDistance or None
Fitted `distance_estimator`. None when `distance_from_prior=True` and
fewer than two assets are investable.
seriation_estimator_ : BaseSeriation
Fitted ordering estimator. Its `ordering_` contains positions in the
full input schema. Hierarchical linkage diagnostics are available through
`seriation_estimator_.hierarchical_clustering_estimator_`.
hierarchical_clustering_estimator_ : HierarchicalClustering or None
Deprecated alias of `seriation_estimator_.hierarchical_clustering_estimator_`.
Access emits a FutureWarning. Will be removed in version 2.0.
Unavailable with spectral seriation.
investable_mask_ : ndarray of shape (n_assets,) or None
Mask of investable assets from the fitted prior. None when all assets are
investable. May be absent after batch fallback if fitting failed before
the prior's investable mask was determined.
n_features_in_ : int
Number of assets seen during `fit`.
feature_names_in_ : ndarray of shape (`n_features_in_`,)
Names of assets seen during `fit`. Defined only when `X`
has asset names that are all strings.
fallback_ : BaseOptimization | "previous_weights" | None
The fallback estimator instance, or the string `"previous_weights"`, that
produced the final result. `None` if no fallback was used.
fallback_chain_ : list[tuple[str, str]] | None
Sequence describing the optimization fallback attempts. Each element is a
pair `(estimator_repr, outcome)` where `estimator_repr` is the string
representation of the primary estimator or a fallback (e.g. `"EqualWeighted()"`,
`"previous_weights"`), and `outcome` is `"success"` if that step produced
a valid solution, otherwise the stringified error message. For successful
fits without any fallback, this is `None`.
error_ : str | list[str | None] | None
For a single portfolio, this is the recorded error message, or None after a
successful allocation or fallback. For multiple portfolios, it is a list with
one entry per row of `weights_`, containing an error message for each failed
portfolio and None for each successful portfolio.
References
----------
.. [1] "Schur Complementary Allocation: A Unification of Hierarchical Risk Parity
and Minimum Variance Portfolios". Peter Cotton (2024).
.. [2] "Portfolio Optimization. Theory and Application".
Chapter 12.3.4 "From Portfolio Risk Minimization to Hierarchical Portfolios"
Daniel P. Palomar (2025).
.. [3] "Building diversified portfolios that outperform out of sample",
The Journal of Portfolio Management,
Marcos López de Prado (2016).
.. [4] "A robust estimator of the efficient frontier",
SSRN Electronic Journal,
Marcos López de Prado (2019).
.. [5] "Machine Learning for Asset Managers",
Elements in Quantitative Finance. Cambridge University Press,
Marcos López de Prado (2020).
.. [6] "A review of two decades of correlations, hierarchies, networks and
clustering in financial markets",
Gautier Marti, Frank Nielsen, Mikołaj Bińkowski, Philippe Donnat (2020).
Examples
--------
For a full tutorial on Schur Complementary Allocation, see
:ref:`sphx_glr_auto_examples_clustering_plot_6_schur.py`.
For online updates with spectral seriation, late listings, delistings,
holidays and asset warm-up, see
:ref:`sphx_glr_auto_examples_online_learning_plot_online_schur_changing_universe.py`.
>>> from skfolio import RiskMeasure
>>> from skfolio.cluster import HierarchicalClustering, LinkageMethod
>>> from skfolio.datasets import load_sp500_dataset
>>> from skfolio.distance import KendallDistance
>>> from skfolio.moments import LedoitWolf
>>> from skfolio.optimization import SchurComplementary
>>> from skfolio.preprocessing import prices_to_returns
>>> from skfolio.prior import EmpiricalPrior
>>> from skfolio.seriation import HierarchicalSeriation
>>>
>>> prices = load_sp500_dataset()
>>> X = prices_to_returns(prices)
>>>
>>> # Default Schur Complementary allocation
>>> model = SchurComplementary(gamma=0.5)
>>> model.fit(X)
SchurComplementary()
>>> print(model.weights_)
[0.0358 0.0061 0.0262 ... 0.0426 0.1129 0.048 ]
>>>
>>> # Advanced model:
>>> # * Ledoit-Wolf covariance shrinkage
>>> # * Kendall's tau distance (absolute) for asset co-dependence
>>> # * Hierarchical clustering with Ward's linkage
>>> model = SchurComplementary(
... gamma=0.5,
... prior_estimator=EmpiricalPrior(covariance_estimator=LedoitWolf()),
... distance_estimator=KendallDistance(absolute=True),
... seriation_estimator=HierarchicalSeriation(
... hierarchical_clustering_estimator=HierarchicalClustering(
... linkage_method=LinkageMethod.WARD,
... ),
... ),
... )
>>> model.fit(X)
SchurComplementary(...)
>>> print(model.weights_)
[0.0323 0.0095 0.0234 ... 0.0402 0.0515 0.0605]
"""
effective_gamma_: float | None
def __init__(
self,
gamma: float = 0.5,
keep_monotonic: bool = True,
prior_estimator: BasePrior | None = None,
distance_estimator: BaseDistance | None = None,
# TODO remove deprecated hierarchical_clustering_estimator in v2.0
hierarchical_clustering_estimator: HierarchicalClustering | None = None,
min_weights: skt.MultiInput | None = 0.0,
max_weights: skt.MultiInput | None = 1.0,
transaction_costs: skt.MultiInput = 0.0,
management_fees: skt.MultiInput = 0.0,
previous_weights: skt.MultiInput | None = None,
portfolio_params: dict | None = None,
fallback: skt.Fallback = None,
raise_on_failure: bool = True,
*,
seriation_estimator: BaseSeriation | None = None,
distance_from_prior: bool = True,
) -> None:
super().__init__(
prior_estimator=prior_estimator,
distance_estimator=distance_estimator,
# TODO remove deprecated hierarchical_clustering_estimator in v2.0
hierarchical_clustering_estimator=hierarchical_clustering_estimator,
min_weights=min_weights,
max_weights=max_weights,
transaction_costs=transaction_costs,
management_fees=management_fees,
previous_weights=previous_weights,
portfolio_params=portfolio_params,
fallback=fallback,
raise_on_failure=raise_on_failure,
seriation_estimator=seriation_estimator,
distance_from_prior=distance_from_prior,
)
self.gamma = gamma
self.keep_monotonic = keep_monotonic
[docs]
def fit(
self, X: ArrayLike, y: ArrayLike | None = None, **fit_params: Any
) -> SchurComplementary:
"""Fit the Schur Complementary estimator.
Parameters
----------
X : array-like of shape (n_observations, n_assets)
Price returns of the assets.
y : array-like, optional
Targets passed to the prior and, when distance_from_prior=False,
to the distance estimator.
**fit_params : dict
Metadata routed to the underlying estimators. Metadata supplied to
a distance fitted on prior scenarios must align with those scenarios.
Returns
-------
self : SchurComplementary
Fitted estimator.
"""
self._reset()
self.effective_gamma_ = None
self._warn_deprecated_clustering_estimator(stacklevel=4)
return self._fit(X, y, method="fit", **fit_params)
[docs]
def partial_fit(
self, X: ArrayLike, y: ArrayLike | None = None, **fit_params: Any
) -> SchurComplementary:
"""Update the prior and compute a Schur allocation.
The prior must support incremental learning. With
`distance_from_prior=False`, the distance must also support it.
Supply only new observations on each call. See
:ref:`online_failure_handling` for failure handling and restarts.
Parameters
----------
X : array-like of shape (n_observations, n_assets)
New returns with the same full asset schema as previous calls.
y : array-like, optional
Targets passed to the prior and, when distance_from_prior=False,
to the distance estimator.
**fit_params : dict
Metadata routed to the underlying estimators. Requires metadata
routing to be enabled. Observation metadata is not forwarded to
distances fitted on the prior's return scenarios or covariance.
Returns
-------
self : SchurComplementary
Updated estimator.
Raises
------
OptimizationError
If optimization fails with `raise_on_failure=True` and `fallback=None`.
"""
self._warn_deprecated_clustering_estimator(stacklevel=3)
return self._fit(X, y, method="partial_fit", **fit_params)
def _fit(
self, X: ArrayLike, y: ArrayLike | None, method: str, **fit_params: Any
) -> SchurComplementary:
"""Publish gamma only for a successful primary Schur allocation."""
try:
super()._fit(X, y, method, **fit_params)
except OptimizationError:
self.effective_gamma_ = None
raise
if self.weights_ is None or self.fallback_ is not None:
self.effective_gamma_ = None
elif len(self.seriation_estimator_.ordering_) == 1:
self.effective_gamma_ = 0.0
return self
def _validate_params(self) -> None:
"""Validate parameters."""
if not 0.0 <= self.gamma <= 1.0:
raise ValueError(f"gamma must be between 0 and 1. Got {self.gamma}")
def _compute_weights(
self,
return_distribution: ReturnDistribution,
ordering: IntArray,
min_weights: FloatArray,
max_weights: FloatArray,
) -> FloatArray:
"""Allocate on the ordered covariance using the existing Schur recursion."""
covariance = cov_nearest(return_distribution.covariance)
if self.keep_monotonic:
weights, self.effective_gamma_ = _compute_monotonic_weights(
max_gamma=self.gamma,
sorted_assets=ordering,
covariance=covariance,
min_weights=min_weights,
max_weights=max_weights,
)
else:
weights = _compute_weights(
gamma=self.gamma,
sorted_assets=ordering,
covariance=covariance,
min_weights=min_weights,
max_weights=max_weights,
force_spd=True,
)
self.effective_gamma_ = self.gamma
if weights is None:
raise OptimizationError("Schur allocation could not produce valid weights.")
return weights
def _compute_monotonic_weights(
max_gamma: float,
sorted_assets: IntArray,
covariance: FloatArray,
max_weights: FloatArray,
min_weights: FloatArray,
step: float = 0.1,
tol: float = 1e-4,
) -> tuple[FloatArray | None, float]:
"""Find the gamma value corresponding to the turning point where portfolio risk
(variance) stops decreasing monotonically.
This method exploits the assumed smooth (i.e., continuously differentiable)
functional dependence of portfolio variance on the regularization parameter
gamma in Schur-complement-based optimization. It searches for the first gamma
value (up to `max_gamma`) beyond which variance no longer decreases. Within
the interval selected by the initial sweep, the search assumes a decrease
followed by stabilization or an increase.
Parameters
----------
max_gamma : float
Maximum gamma value to sweep up to.
sorted_assets : ndarray of shape (n_assets,)
Array of ordered asset indices.
covariance : FloatArray
Covariance matrix of asset returns.
max_weights : FloatArray
Maximum allowable weights for each asset.
min_weights : FloatArray
Minimum allowable weights for each asset.
step : float, default=0.1
Step size for incrementing gamma during the initial sweep.
tol : float, default=1e-4
Gamma interval tolerance after finding a decreasing point, also used as
the backward-difference step. Refinement continues below this tolerance
when no decreasing point has been found, up to the iteration limit.
Returns
-------
weights : FloatArray or None
Asset weights at the identified turning point, or `None` if no weights can
be computed.
effective_gamma : float
Estimated gamma value corresponding to the first variance turning point,
or `max_gamma` when variance keeps decreasing. The returned weights are
computed at this value.
"""
if max_gamma == 0:
weights = _compute_weights(
gamma=0,
sorted_assets=sorted_assets,
covariance=covariance,
max_weights=max_weights,
min_weights=min_weights,
force_spd=True,
)
return weights, 0.0
def objective(x: float) -> tuple[float, FloatArray | None]:
"""Return the portfolio variance and weights for `gamma=x`.
The variance is `inf` and the weights are `None` if no weights can be computed.
"""
w = _compute_weights(
gamma=x,
sorted_assets=sorted_assets,
covariance=covariance,
max_weights=max_weights,
min_weights=min_weights,
force_spd=False,
)
risk = np.inf if w is None else w @ covariance @ w.T
return risk, w
# Evenly spaced gamma vector in [0, max_gamma]
n = int(np.ceil(max_gamma / step)) + 1
gammas = np.linspace(0, max_gamma, n)
variances = np.full_like(gammas, np.nan)
# Initial sweep of the discrete gamma vector in [0, max_gamma] to find the range
# of the variance turning point if any.
variance, weights = objective(gammas[0])
variances[0] = variance
previous_weights = weights
for i in range(1, n):
# Retain the two preceding sweep allocations for the refinement interval.
lower_weights, previous_weights = previous_weights, weights
variance, weights = objective(gammas[i])
variances[i] = variance
if variance >= variances[i - 1]:
# Refine the turning point in [gammas[i-2], gammas[i]], or
# [0, gammas[1]] on the first step. If variance is non-decreasing
# throughout the interval, retain the lower endpoint's allocation.
low_index = max(0, i - 2)
return _binary_search(
objective,
low_gamma=gammas[low_index],
high_gamma=gammas[i],
low_variance=variances[low_index],
low_weights=lower_weights,
tol=tol,
)
# No turning point found in sweep
# 1) Check local derivative at the terminal gamma
# Keep the backward difference within [0, max_gamma], reusing variance at zero.
variance_h = variances[0] if max_gamma <= tol else objective(max_gamma - tol)[0]
if variance <= variance_h:
# monotonically decreasing up to max_gamma --> we return the terminal gamma
return weights, max_gamma
# 2) Turning point lies between last two gammas, we refine its location by
# binary search
return _binary_search(
objective,
low_gamma=gammas[-2],
high_gamma=max_gamma,
low_variance=variances[-2],
low_weights=previous_weights,
tol=tol,
)
def _binary_search(
objective: Callable,
low_gamma: float,
high_gamma: float,
low_variance: float,
low_weights: FloatArray | None,
tol: float = 1e-4,
) -> tuple[FloatArray | None, float]:
"""Locate the variance turning point by binary search.
Search the interval [low_gamma, high_gamma] for the point where portfolio
variance stops decreasing monotonically. This method assumes that variance
decreases smoothly with gamma up to a point, after which it stabilizes or
increases. It evaluates the objective at midpoints and estimates the local
slope with a backward difference to identify this transition with gamma
precision controlled by `tol`.
After accepting a decreasing point, stop at a feasible midpoint once the
interval width is at most `tol`. Return the last accepted lower endpoint and
its weights. If no decreasing point has been found, continue refining up to
the iteration limit to detect improvements near the initial lower endpoint.
If subsequent midpoints remain infeasible, return the accepted lower endpoint
when the iteration limit is reached. If no midpoint is accepted, return the
initial lower endpoint and its supplied weights.
Parameters
----------
objective : callable
A function that takes a float gamma value and returns a tuple:
(variance: float, weights: FloatArray or None). Infeasible allocations
return infinite variance and None for the weights.
low_gamma : float
Lower bound of the gamma search interval.
high_gamma : float
Upper bound of the gamma search interval.
low_variance : float
The variance value corresponding to `low_gamma`.
low_weights : FloatArray or None
Weights at `low_gamma`, or None if that allocation is infeasible.
tol : float, default=1e-4
Gamma interval tolerance after finding a decreasing point, also used as
the backward-difference step, bounded at zero.
Returns
-------
weights : FloatArray or None
Asset weights corresponding to the returned gamma, or None if neither
the initial lower endpoint nor any accepted midpoint is feasible.
gamma : float
Gamma near the variance turning point or feasibility boundary within
the search interval. The returned weights are computed at this value.
"""
max_iter = max(1, math.ceil(math.log2((high_gamma - low_gamma) / tol) * 2 + 1))
initial_gamma, initial_variance = low_gamma, low_variance
for _ in range(max_iter):
mid_gamma = 0.5 * (low_gamma + high_gamma)
variance, weights = objective(mid_gamma)
# A midpoint with higher variance or an infeasible allocation already
# selects the left subinterval. Estimate the local derivative only when
# the midpoint passes the variance comparison.
is_decreasing = weights is not None and variance <= low_variance
if is_decreasing:
gamma_h = max(0.0, mid_gamma - tol)
variance_h = (
initial_variance if gamma_h == initial_gamma else objective(gamma_h)[0]
)
is_decreasing = variance <= variance_h
if is_decreasing:
low_gamma = mid_gamma
low_variance = variance
low_weights = weights
else:
high_gamma = mid_gamma
# Keep refining near the starting point until a new point is accepted,
# even when the interval is already below tol.
if (
low_gamma > initial_gamma
and weights is not None
and (high_gamma - low_gamma) <= tol
):
return low_weights, low_gamma
# At a feasibility boundary, all subsequent midpoints can be infeasible.
# Retain the last accepted allocation, or the supplied starting allocation,
# after exhausting the refinement budget.
return low_weights, low_gamma
def _compute_weights(
gamma: float,
sorted_assets: IntArray,
covariance: FloatArray,
max_weights: FloatArray,
min_weights: FloatArray,
force_spd: bool = False,
) -> FloatArray | None:
"""
Core Schur-complement allocation recursion.
Parameters
----------
gamma : float
Regularization factor.
sorted_assets : ndarray of shape (n_assets,)
Asset indices in dendrogram order.
covariance : ndarray of shape (n_assets, n_assets)
Asset covariance matrix.
min_weights : ndarray of shape (n_assets,)
Minimum weights array.
max_weights : ndarray of shape (n_assets,)
Maximum weights array.
force_spd : bool
If set to False, return None if any augmented block is not SPD; otherwise
convert non-SPD blocks to their nearest SPD during the bisection steps.
Returns
-------
weights : ndarray | None
Final portfolio weights, or None if SPD check fails.
"""
covariance = covariance.copy()
n_assets = len(covariance)
weights = np.ones(n_assets)
items = [sorted_assets]
while len(items) > 0:
new_items = []
for left_cluster, right_cluster in bisection(items):
new_items += [left_cluster, right_cluster]
a = covariance[np.ix_(left_cluster, left_cluster)]
d = covariance[np.ix_(right_cluster, right_cluster)]
if len(left_cluster) <= 1:
a_aug, d_aug = a, d
else:
b = covariance[np.ix_(left_cluster, right_cluster)]
a_aug = _schur_augmentation(a, b, d, gamma=gamma)
d_aug = _schur_augmentation(d, b.T, a, gamma=gamma)
if not force_spd:
if not is_cholesky_dec(a_aug) or not is_cholesky_dec(d_aug):
return None
else:
try:
# A block with non-positive variances cannot be repaired by
# correlation clipping. Report the failure before NaNs spread.
with np.errstate(invalid="raise", divide="raise", over="raise"):
if not is_cholesky_dec(a_aug):
a_aug = cov_nearest(a_aug)
if not is_cholesky_dec(d_aug):
d_aug = cov_nearest(d_aug)
except (ValueError, FloatingPointError, np.linalg.LinAlgError) as error:
raise OptimizationError(
f"Schur complement failed with gamma={gamma:0.4f}. Choose a "
"smaller gamma or set `keep_monotonic=True`"
) from error
# Subsequent splits must use the repaired blocks too.
covariance[np.ix_(left_cluster, left_cluster)] = a_aug
covariance[np.ix_(right_cluster, right_cluster)] = d_aug
left_variance = _naive_portfolio_variance(a_aug)
right_variance = _naive_portfolio_variance(d_aug)
alpha = 1 - left_variance / (left_variance + right_variance)
# Weights constraints
alpha = _apply_weight_constraints_to_split_factor(
alpha=alpha,
weights=weights,
max_weights=max_weights,
min_weights=min_weights,
left_cluster=left_cluster,
right_cluster=right_cluster,
)
weights[left_cluster] *= alpha
weights[right_cluster] *= 1 - alpha
items = new_items
return weights
def _naive_portfolio_variance(covariance: FloatArray) -> float:
"""Portfolio variance of an inverse variance allocation.
Parameters
----------
covariance : ndarray of shape (n_assets, n_assets)
Covariance matrix.
Returns
-------
variance : float
Portfolio variance of an inverse variance allocation.
"""
weights = 1 / np.diag(covariance)
weights /= weights.sum()
variance = weights @ covariance @ weights.T
return variance
def _schur_augmentation(
a: FloatArray, b: FloatArray, d: FloatArray, gamma: float
) -> FloatArray:
"""Compute an augmented covariance matrix `A` inspired by the
Schur complement [1]_.
Parameters
----------
a : ndarray of shape (n1, n1)
Upper left block matrix `A`
b : ndarray of shape (n1, n2)
Upper right block matrix `B`
d : ndarray of shape (n2, n2)
Lower right block matrix `D`
gamma : float
Regularization factor between 0 and 1.
A value of 0 means that no additional information is used from off-diagonal
matrix blocks and is equivalent to a Hierarchical Risk Parity.
As the value increases to 1, the allocation tends to the Minimum Variance
Optimization allocation.
Returns
-------
a_aug : ndarray of shape (n1, n1)
Augmented covariance matrix `A`.
References
----------
.. [1] "Schur Complementary Portfolios - A Unification of Machine Learning and
Optimization-Based Allocation".
Peter Cotton (2022).
"""
n_a = a.shape[0]
n_d = d.shape[0]
if gamma == 0 or n_a == 1 or n_d == 1:
return a
a_aug = a - gamma * b @ inverse_multiply(d, b.T)
m = symmetric_step_up_matrix(n1=n_a, n2=n_d)
r = np.eye(n_a) - gamma * multiply_by_inverse(b, d) @ m.T
a_aug = inverse_multiply(r, a_aug)
# make it symmetric
symmetrize(a_aug)
return a_aug
