Constrained Risk Budgeting Portfolios: Theory, Algorithms, Applications & Puzzles

TL;DR

This paper advances the theory and algorithms for risk budgeting portfolios with constraints, offering practical solutions for investment management and addressing complex mathematical and computational challenges.

Contribution

It introduces a new optimization framework for constrained risk budgeting portfolios and develops scalable algorithms inspired by machine learning techniques.

Findings

Algorithms effectively solve constrained risk budgeting problems

Application to dynamic turnover control in risk parity portfolios

Improved liquidity and reduced small cap bias in smart beta portfolios

Abstract

This article develops the theory of risk budgeting portfolios, when we would like to impose weight constraints. It appears that the mathematical problem is more complex than the traditional risk budgeting problem. The formulation of the optimization program is particularly critical in order to determine the right risk budgeting portfolio. We also show that numerical solutions can be found using methods that are used in large-scale machine learning problems. Indeed, we develop an algorithm that mixes the method of cyclical coordinate descent (CCD), alternating direction method of multipliers (ADMM), proximal operators and Dykstra's algorithm. This theoretical body is then applied to some investment problems. In particular, we show how to dynamically control the turnover of a risk parity portfolio and how to build smart beta portfolios based on the ERC approach by improving the liquidity of the portfolio or reducing the small cap bias. Finally, we highlight the importance of the homogeneity property of risk measures and discuss the related scaling puzzle.

Figures

Tables

Table 1. Table 1: Computation of ERC and RB portfolios

Asset	ERC portfolio				RB portfolio
	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$
$1$	$41.01$	$7.79$	$3.19$	$25.00$	$45.05$	$8.06$	$3.63$	$30.00$
$2$	$27.34$	$11.68$	$3.19$	$25.00$	$30.04$	$12.09$	$3.63$	$30.00$
$3$	$18.99$	$16.82$	$3.19$	$25.00$	$14.67$	$16.10$	$2.36$	$19.50$
$4$	$12.66$	$25.23$	$3.19$	$25.00$	$10.24$	$24.23$	$2.48$	$20.50$
$\sigma \left(x\right)$			$12.78$				$12.11$

Table 2. Table 2: Computation of ERC and RB portfolios when x i ≤ 30 % subscript 𝑥 𝑖 percent 30 x_{i}\leq 30\%

Asset	ERC portfolio				RB portfolio
	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$
$1$	$30.00$	$7.19$	$2.16$	$15.50$	$30.00$	$7.19$	$2.16$	$15.48$
$2$	$30.00$	$11.60$	$3.48$	$24.98$	$30.00$	$11.60$	$3.48$	$24.96$
$3$	$24.57$	$17.43$	$4.28$	$30.74$	$24.43$	$17.42$	$4.26$	$30.52$
$4$	$15.43$	$25.98$	$4.01$	$28.78$	$15.57$	$26.01$	$4.05$	$29.04$
$\sigma \left(x\right)$			$13.93$				$13.94$

Table 3. Table 3: Computational time using our Matlab implementation (relative value)

Algorithm	$x$-update	(1)	(2)	(3)
ADMM	Newton	$2$	$1$	$1$
ADMM	BFGS	$380$	$280$	$25$
ADMM	QP	$220$	$120$	$110$
ADMM	CCD	$10$	$9$	$8$
CCD		$1$	$1$

Table 4. Table 4: Volatility breakdown (in %) of the current portfolio

Asset	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$
1	$25.00$	$10.00$	$2.50$	$20.21$
2	$25.00$	$15.40$	$3.85$	$31.10$
3	$10.00$	$20.30$	$2.03$	$16.41$
4	$10.00$	$22.24$	$2.22$	$17.98$
5	$30.00$	$5.90$	$1.77$	$14.30$
$\sigma \left(x\right)$			$12.37$

Table 5. Table 5: Volatility breakdown (in %) of the ERC portfolio

Asset	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$
1	$22.40$	$10.61$	$2.38$	$20.00$
2	$16.51$	$14.39$	$2.38$	$20.00$
3	$12.03$	$19.74$	$2.38$	$20.00$
4	$10.51$	$22.60$	$2.38$	$20.00$
5	$38.54$	$6.16$	$2.38$	$20.00$
$\sigma \left(x\right)$			$11.88$

Table 6. Table 6: Volatility breakdown (in %) of the constrained RB portfolio

Asset	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$	${\lambda }_i^{-}$	${\lambda }_i^{+}$
1	$22.89$	$10.28$	$2.35$	$19.39$	$0.00$	$0.00$
2	$20.00$	$14.90$	$2.98$	$24.55$	$3.13$	$0.00$
3	$11.69$	$20.13$	$2.35$	$19.39$	$0.00$	$0.00$
4	$10.42$	$22.57$	$2.35$	$19.39$	$0.00$	$0.00$
5	$35.00$	$6.00$	$2.10$	$17.29$	$0.00$	$0.73$
$\sigma \left(x\right)$			$12.14$		$\lambda =11.76$

Table 7. Table 7: Volatility breakdown (in %) of naive and least squares solutions

	Naive solution				Least squares solution
Asset	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℳ{ℛ}_i$	$ℛ{𝒞}_i$	$ℛ{𝒞}_i^{\star }$
1	$22.84$	$10.25$	$2.34$	$19.30$	$23.13$	$10.32$	$2.39$	$19.70$
2	$20.00$	$14.98$	$3.00$	$24.70$	$20.00$	$14.86$	$2.97$	$24.53$
3	$12.34$	$20.18$	$2.49$	$20.53$	$11.39$	$20.07$	$2.29$	$18.87$
4	$9.83$	$22.46$	$2.21$	$18.20$	$10.48$	$22.55$	$2.36$	$19.51$
5	$35.00$	$5.99$	$2.10$	$17.28$	$35.00$	$6.02$	$2.11$	$17.39$
$\sigma \left(x\right)$			$12.13$				$12.11$

Table 8. Table 8: Volatility and correlation matrix of asset returns (in %)

${\sigma }_i$		1	2	3	4	5	6	7	8
		$5.0$	$5.0$	$7.0$	$10.0$	$15.0$	$15.0$	$15.0$	$18.0$
${\rho }_{i,j}$	1	$100$
	2	$80$	$100$
	3	$60$	$40$	$100$
	4	$-20$	$-20$	$50$	$100$
\cdashline2-10	5	$-10$	$-20$	$30$	$60$	$100$
	6	$-20$	$-10$	$20$	$60$	$90$	$100$
	7	$-20$	$-20$	$20$	$50$	$70$	$60$	$100$
	8	$-20$	$-20$	$30$	$60$	$70$	$70$	$70$	$100$

Table 9. Table 9: The case of inequality constraints

$x_5+x_6+x_7+x_8\ge 30\%$			$✓$		$✓$
$x_2+x_6\ge x_1+x_5+5\%$					$✓$
Asset	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$
1	$26.83$	$12.50$	$25.78$	$8.64$	$24.52$	$8.16$
2	$28.68$	$12.50$	$27.41$	$8.64$	$28.69$	$9.13$
3	$11.41$	$12.50$	$9.51$	$8.64$	$9.52$	$8.61$
4	$9.80$	$12.50$	$7.29$	$8.64$	$7.27$	$8.61$
5	$5.61$	$12.50$	$7.06$	$15.91$	$6.97$	$15.69$
6	$5.90$	$12.50$	$7.71$	$16.58$	$7.80$	$16.82$
7	$6.66$	$12.50$	$9.23$	$18.14$	$9.23$	$18.16$
8	$5.11$	$12.50$	$6.00$	$14.82$	$6.00$	$14.81$
$\sigma \left(x\right)$ (in %)	$4.78$		$5.20$		$5.19$

Table 10. Table 10: Volatility breakdown (in %) of constrained ERC portfolios

Asset	CW		ERC		LC-ERC		LS-ERC		C-ERC
	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$
1	$34.00$	$32.08$	$17.22$	$14.29$	$25.81$	$23.39$	$26.62$	$24.23$	$25.87$	$23.46$
2	$25.00$	$24.82$	$15.90$	$14.29$	$24.06$	$23.44$	$24.20$	$23.63$	$24.07$	$23.46$
3	$20.00$	$20.92$	$14.78$	$14.29$	$22.44$	$23.44$	$22.09$	$23.08$	$22.46$	$23.46$
4	$15.00$	$16.01$	$13.83$	$14.29$	$21.69$	$23.57$	$21.09$	$22.89$	$21.59$	$23.46$
5	$3.00$	$3.10$	$13.17$	$14.29$	$3.00$	$3.10$	$3.00$	$3.10$	$3.00$	$3.10$
6	$2.00$	$2.03$	$12.86$	$14.29$	$2.00$	$2.02$	$2.00$	$2.02$	$2.00$	$2.02$
7	$1.00$	$1.05$	$12.23$	$14.29$	$1.00$	$1.05$	$1.00$	$1.05$	$1.00$	$1.05$
$\sigma \left(x\right)$	$14.50$		$15.23$		$14.68$		$14.66$		$14.68$

Table 11. Table 11: Constrained RB portfolios (in %) with turnover control

Asset	${\tau }^{\star }$
	$0.00$	$10.00$	$20.00$	$30.00$	$40.00$	$50.00$	$60.00$	$70.00$
1	$12.50$	$14.86$	$17.28$	$19.68$	$22.01$	$24.28$	$26.58$	$26.83$
2	$12.50$	$15.14$	$17.72$	$20.32$	$22.99$	$25.72$	$28.42$	$28.68$
3	$12.50$	$12.50$	$12.50$	$12.50$	$12.50$	$12.50$	$11.65$	$11.41$
4	$12.50$	$12.50$	$12.50$	$12.50$	$12.50$	$11.50$	$9.90$	$9.80$
5	$12.50$	$11.20$	$9.70$	$8.49$	$7.27$	$6.28$	$5.66$	$5.61$
6	$12.50$	$12.02$	$10.36$	$9.02$	$7.69$	$6.63$	$5.95$	$5.90$
7	$12.50$	$12.50$	$11.72$	$10.16$	$8.66$	$7.47$	$6.71$	$6.66$
8	$12.50$	$9.28$	$8.22$	$7.33$	$6.39$	$5.62$	$5.14$	$5.11$
$\tau \left(x^{\star }\mid x_0\right)$	$0.00$	$10.00$	$20.00$	$30.00$	$40.00$	$50.00$	$60.00$	$61.02$

Table 12. Table 12: Illustration of the scaling puzzle

$\sum _{i=5}^8{x}_i\ge 30\%$	$✓$
$\sum _{i=1}^4{x}_i\le 70\%$			$✓$
\hdashline$\sum _{i=5}^8{x}_i\ge 40\%$					$✓$
$\sum _{i=1}^4{x}_i\le 60\%$							$✓$
Asset	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$	$x_i$	$ℛ{𝒞}_i^{\star }$
1	$25.78$	$8.64$	$23.39$	$6.50$	$24.09$	$4.35$	$18.73$	$2.01$
2	$27.41$	$8.64$	$24.34$	$6.11$	$25.09$	$4.35$	$19.08$	$1.68$
3	$9.51$	$8.64$	$12.46$	$10.98$	$6.57$	$4.35$	$12.48$	$7.84$
4	$7.29$	$8.64$	$9.81$	$12.07$	$4.24$	$4.35$	$9.71$	$10.43$
\hdashline5	$7.06$	$15.91$	$7.30$	$16.09$	$8.74$	$18.59$	$9.82$	$19.51$
6	$7.71$	$16.58$	$7.66$	$16.09$	$10.75$	$21.87$	$10.27$	$19.51$
7	$9.23$	$18.14$	$8.46$	$16.09$	$14.09$	$27.32$	$11.15$	$19.51$
8	$6.00$	$14.82$	$6.57$	$16.09$	$6.42$	$14.82$	$8.76$	$19.51$
$\sigma \left(x\right)$ (in %)	$5.20$		$5.43$		$5.98$		$6.56$
$ℒ(x^{\star };{\lambda }^{\star })$ (in %)	$13.29$		$20.86$		$10.68$		$28.27$
Global minimum	$✓$				$✓$

Equations

$$\mathcal{R}\left(x\right)=-x^{\top }\left(\mu -r\right)+c\cdot \sqrt{x^{\top }\Sigma x}$$

$${\mathcal{R}\mathcal{C}}_i\left(x\right)=x_i\cdot \left(-\left({\mu }_i-r\right)+c\frac{{\left(\Sigma x\right)}_i}{\sqrt{x^{\top }\Sigma x}}\right)$$

$$\left\{\begin{array}[]{l}\mathcal{RC}_{i}\left(x\right)=b_{i}\mathcal{R}\left(x\right)\\ b_{i}>0\\ x_{i}\geq 0\\ \sum_{i=1}^{n}b_{i}=1\\ \sum_{i=1}^{n}x_{i}=1\end{array}\right.$$

$${\mathrm{SR}}^{+}=\max \left(\underset{x\in {\left[0,1\right]}^n}{\sup }\mathrm{SR}\left(x\mid r\right),0\right)$$

$${\mathcal{R}\mathcal{C}}_i\left(x\right)=\frac{x_i\cdot {\left(\Sigma x\right)}_i}{\sqrt{x^{\top }\Sigma x}}=\frac{1}{n}\sqrt{x^{\top }\Sigma x}$$

$$\frac{1}{b_i}{\mathcal{R}\mathcal{C}}_i\left(x\right)=\frac{1}{b_j}{\mathcal{R}\mathcal{C}}_j\left(x\right)\text{for all }i,j$$

$$x_{\mathrm{RB}}$$

$$\displaystyle\left\{\begin{array}[]{l}\mathbf{1}^{\top}x=1\\ x\geq\mathbf{0}\end{array}\right.$$

$$x_{\mathrm{RB}}$$

$$\displaystyle\left\{\begin{array}[]{l}\sum_{i=1}^{n}b_{i}\ln x_{i}\geq\kappa^{\star}\\ \mathbf{1}^{\top}x=1\\ x\geq\mathbf{0}\end{array}\right.$$

$$x^{\star }\left(\kappa \right)$$

$$\displaystyle\left\{\begin{array}[]{l}\sum_{i=1}^{n}b_{i}\ln x_{i}\geq\kappa\\ x\geq\mathbf{0}\end{array}\right.$$

$$x_{\mathrm{RB}}=\frac{x^{\star }\left(\kappa \right)}{1^{\top }{x}^{\star }\left(\kappa \right)}$$

$$x^{\star }\left(\lambda \right)$$

$$x_{\mathrm{RB}}=\frac{x^{\star }\left(\lambda \right)}{1^{\top }{x}^{\star }\left(\lambda \right)}$$

$$\frac{\partial f\left(x\right)}{\partial x_i}$$

$$\frac{{\partial }^{2}f\left(x\right)}{\partial x_i\partial x_j}$$

$$\frac{{\partial }^{2}f\left(x\right)}{\partial x_i^2}$$

$$x^{\left(k+1\right)}=x^{\left(k\right)}-{\eta }^{\left(k\right)}{\left(\frac{{\partial }^{2}f\left(x^{\left(k\right)}\right)}{\partial x\partial x^{\top }}\right)}^{-1}\frac{\partial f\left(x^{\left(k\right)}\right)}{\partial x}$$

$$x^{\left(k+1\right)}$$

$$x_i^{\left(k+1\right)}$$

$$x_i^{\left(k+1\right)}=\underset{x}{\mathrm{argmin}}f\left(x_1^{\left(k+1\right)},\dots ,x_{i-1}^{\left(k+1\right)},x,x_{i+1}^{\left(k\right)},\dots ,x_n^{\left(k\right)}\right)$$

$$f\left(x\right)=-x^{\top }\pi +c\sqrt{x^{\top }\Sigma x}-\lambda \sum _{i=1}^n{b}_i\ln x_i$$

$$\frac{\partial \mathcal{L}\left(x;\lambda \right)}{\partial x_i}=-{\pi }_i+c\frac{{\left(\Sigma x\right)}_i}{\sigma \left(x\right)}-\lambda \frac{b_i}{x_i}$$

$$c{\sigma }_i^2{x}_i^2+\left(c{\sigma }_i\sum _{j\ne i}{x}_j{\rho }_{i,j}{\sigma }_j-{\pi }_i\sigma \left(x\right)\right)x_i-\lambda b_i\sigma \left(x\right)=0$$

$$x_i=\frac{-c\left({\sigma }_i{\sum }_{j\ne i}{x}_j{\rho }_{i,j}{\sigma }_j\right)+{\pi }_i\sigma \left(x\right)+\sqrt{{\left(c\left({\sigma }_i{\sum }_{j\ne i}{x}_j{\rho }_{i,j}{\sigma }_j\right)-{\pi }_i\sigma \left(x\right)\right)}^2+4\lambda c{b}_i{\sigma }_i^2\sigma \left(x\right)}}{2c{\sigma }_i^2}$$

$$\left\{\begin{array}[]{l}\mathcal{RC}_{i}\left(x\right)=b_{i}\mathcal{R}\left(x\right)\\ x\in\mathcal{S}\\ x\in\Omega\end{array}\right.$$

$$\mathcal{S}=\left\{x_i\ge 0:\sum _{i=1}^n{x}_i=1\right\}$$

$$\left\{\begin{array}[]{l}\mathcal{RC}_{i}\left(x\right)\approx b_{i}\mathcal{R}\left(x\right)\\ x\in\mathcal{S}\\ x\in\Omega\end{array}\right.$$

$$x^{\star }\left(\mathcal{S},\Omega \right)$$

$$\left\{x^{\star }\left(\mathcal{S},\Omega \right),{\theta }^{\star }\right\}$$

$$\rho=\left(\begin{array}[]{cccc}1.00&&&\\ 0.50&1.00&&\\ 0.50&0.50&1.00&\\ 0.50&0.50&0.75&1.00\end{array}\right)$$

$$x^{\star }\left(\mathcal{S},\Omega \right)$$

Full text

Keywords: Risk budgeting, large-scale optimization, Lagrange function, cyclical coordinate descent (CCD), alternating direction method of multipliers (ADMM), proximal operator, Dykstra’s algorithm, turnover, liquidity, risk parity, smart beta portfolio.

JEL classification: C61, G11.

1 Introduction

Since the 2008 Global Financial Crisis, the development of risk budgeting (RB) techniques has marked an important milestone in portfolio management by putting diversification at the center of portfolio construction (Qian, 2005; Maillard et al., 2010). In particular, the equal risk contribution (ERC) portfolio has been very popular and has significantly impacted the asset management industry. For instance, this allocation approach is extensively implemented in risk parity funds, factor investing and alternative risk premia (Roncalli, 2017).

The main advantages of RB portfolios are the stability of the allocation, and the diversification management principle, which appear more robust than the diversification mechanism of mean-variance optimized portfolios (Bourgeron et al., 2018). This is why we don’t need to add constraints in order to regularize the solution. This advantage is also its drawback. Indeed, there are some situations where portfolio managers have to impose constraints. For example, they may want to impose a minimum investment weight, a sector-neutrality or some liquidity thresholds. The goal of this paper is then to define what does a constrained risk budgeting portfolio mean, since adding constraints will change the risk budgets that are targeted, meaning that ex-post risk contributions are not equal to ex-ante risk budgets.

This paper is organized as follows. Section Two illustrates the bridge between risk budgeting and portfolio optimization. In Section Three, we present the right mathematical formulation of constrained risk budgeting portfolios, and develop the corresponding numerical algorithms. In Section Four, we consider some applications in asset allocation, in particular the management of turnover or the consideration of liquidity. Finally, we discuss the compatibility puzzle of the homogeneity property of coherent risk measures.

2 The original risk budgeting portfolio

2.1 Definition of the risk budgeting portfolio

We consider a universe of $n$ risky assets. Let $\mu$ and $\Sigma$ be the vector of expected returns and the covariance matrix of asset returns. We have $\Sigma_{i,j}=\rho_{i,j}\sigma_{i}\sigma_{j}$ where $\sigma_{i}$ is the volatility of asset $i$ and $\rho_{i,j}$ is the correlation between asset $i$ and asset $j$. Following Roncalli (2015), we consider the standard deviation-based risk measure defined as follows:

[TABLE]

where $c$ is a scalar that measures the trade-off between the expected return of the portfolio and its volatility. We deduce that the risk contribution of Asset $i$ is given by:

[TABLE]

Following Maillard et al. (2010), Roncalli (2013) defines the risk budgeting (RB) portfolio using the following non-linear system:

[TABLE]

where $b_{i}$ is the risk budget of Asset $i$ expressed in relative terms. The constraint $b_{i}>0$ implies that we cannot set some risk budgets to zero. This restriction is necessary in order to ensure that the RB portfolio is unique.

Remark 1

Roncalli (2015) shows that the existence of the RB portfolio depends on the value taken by $c$. In particular, the RB portfolio exists and is unique if $c>\mathop{\rm SR}^{+}$ where $\mathop{\rm SR}^{+}$ is the maximum Sharpe ratio of the asset universe:

[TABLE]

Remark 2

The original ERC portfolio is obtained by considering the volatility risk measure and the same risk budgets. It is equivalent to seting $\mu_{i}=r$, $c=1$ and $b_{i}=1/n$. In this case, we have:

[TABLE]

2.2 The associated optimization problem

2.2.1 The wrong formulation

System (2) is equivalent to solving $n$ non-linear equations with $n$ unknown variables. Therefore, we can use the Newton or Broyden methods to find the numerical solution. We also deduce that:

[TABLE]

In order to find the solution, an alternative approach is to solve the optimization problem:

[TABLE]

At the optimum $x_{\mathrm{RB}}$, the objective function $f\left(x\right)$ must be equal to zero. This approach was originally proposed by Maillard et al. (2010) in the case $b_{i}=b_{j}$.

At first sight, Problem (3) seems to be easy to solve because it resembles how a quadratic functions in $n$ variables is defined and we can analytically compute the gradient vector and the hessian matrix of the objective function. In fact, it is not a convex problem (Feng and Palomar, 2015), and the optimization is tricky when the number of assets is large. From a theoretical point of view, the objective function is not well defined because the solution is only valid if the zero can be reached: $f\left(x_{\mathrm{RB}}\right)=0$. But the most important issue is that the optimization problem is driven by the equality constraint: $\mathbf{1}^{\top}x=1$. Indeed, if we remove it, the solution is $x_{\mathrm{RB}}=\mathbf{0}$.

2.2.2 The right formulation

Roncalli (2013) shows that the RB portfolio is the solution of the following optimization problem:

[TABLE]

where $\kappa^{\star}$ is a constant to be determined. This optimization program is equivalent to finding the optimal solution $x^{\star}\left(\kappa\right)$:

[TABLE]

where $\kappa$ is an arbitrary constant and to scale the solution:

[TABLE]

Using the Lagrange formulation, we obtain an equivalent solution:

[TABLE]

where $\lambda$ is an arbitrary positive scalar and:

[TABLE]

$x^{\star}\left(\lambda\right)$ is the solution of a standard logarithmic barrier problem, which has very appealing characteristics. First, it defines a unique solution. Second, the constraint $\mathbf{1}^{\top}x=1$ is removed, meaning that the optimization exploits the scaling property. Finally, the constraint $x\geq\mathbf{0}$ is redundant since the logarithm is defined for strictly positive numbers.

We claim that Problem (2.2.2) is the right risk budgeting problem. For instance, Maillard et al. (2010) used this formulation to show that the ERC portfolio exists and is unique. Roncalli (2013) also noticed that there is a discontinuity when one or more risk budgets $b_{i}$ are equal to zero. In this case, we can find several solutions that satisfy $\mathcal{RC}_{i}\left(x\right)=b_{i}\mathcal{R}\left(x\right)$ or Problem (3), but only one solution if we consider the logarithmic barrier program.

2.3 Numerical solution

2.3.1 The Newton algorithm

Spinu (2013) proposes solving Problem (2.2.2) by using the Newton algorithm222The first and second derivatives are computed using the following analytical expressions:

$\displaystyle\frac{\partial\,f\left(x\right)}{\partial\,x_{i}}$ $\displaystyle=$ $\displaystyle-\left(\mu_{i}-r\right)+c\frac{\left(\Sigma x\right)_{i}}{\sqrt{x^{\top}\Sigma x}}-\lambda\frac{b_{i}}{x_{i}}$

$\displaystyle\frac{\partial^{2}\,f\left(x\right)}{\partial\,x_{i}\,\partial\,x_{j}}$ $\displaystyle=$ $\displaystyle c\frac{\rho_{i,j}\sigma_{i}\sigma_{j}\sqrt{x^{\top}\Sigma x}-\left(\Sigma x\right)_{i}\left(\Sigma x\right)_{j}}{x^{\top}\Sigma x}$

$\displaystyle\frac{\partial^{2}\,f\left(x\right)}{\partial\,x_{i}^{2}}$ $\displaystyle=$ $\displaystyle c\frac{\sigma_{i}^{2}\sqrt{x^{\top}\Sigma x}-\left(\Sigma x\right)_{i}^{2}}{x^{\top}\Sigma x}+\lambda\frac{b_{i}}{x_{i}^{2}}$

:

[TABLE]

where $\eta^{\left(k\right)}\in\left[0,1\right]$ is the step size and $k$ is the iteration index. Generally, we set $\eta^{\left(k\right)}=1$. Spinu (2013) noticed that the Newton algorithm may be improved because the risk measure is self-concordant. In this case, we can use the results of Nesterov (2004) to determine the optimal size $\eta^{\left(k\right)}$ at each iteration.

2.3.2 The CCD algorithm

The descent algorithm is defined by the following rule:

[TABLE]

At the $k^{\mathrm{th}}$ Iteration, the current solution $x^{\left(k\right)}$ is updated by going in the opposite direction to $D^{\left(k\right)}$. For instance, $D^{\left(k\right)}$ is equal respectively to $\partial_{x}\,f\left(x^{\left(k\right)}\right)$ in the gradient algorithm and $\left(\partial_{x,x}^{2}\,f\left(x^{\left(k\right)}\right)\right)^{-1}\partial_{x}\,f\left(x^{\left(k\right)}\right)$ in the Newton algorithm. Coordinate descent is a modification of the descent algorithm by minimizing the function along one coordinate at each step:

[TABLE]

The coordinate descent algorithm becomes a scalar problem, and we know that minimizing a function with respect to one variable is easier than with $n$ variables. Concerning the choice of the variable $i$, there are two approaches: random coordinate descent or RCD (Nesterov, 2012) and cyclical coordinate descent or CCD (Tseng, 2001). In the first case, we assign a random number between $1$ and $n$ to the index $i$. In the second case, we cyclically iterate through the coordinates:

[TABLE]

This ensures that all the indices are selected during one cycle. In the CCD algorithm, $k$ is the cycle index while $i$ is the iteration index within a cycle.

Griveau-Billion et al. (2013) propose applying the CCD algorithm to find the solution of the objective function:

[TABLE]

where $\pi=\mu-r$. The first-order condition is:

[TABLE]

At the optimum, we have $\partial_{x_{i}}\,\mathcal{L}\left(x;\lambda\right)=0$ or:

[TABLE]

By definition of the RB portfolio we have $x_{i}>0$. We notice that the polynomial function is convex because we have $\sigma_{i}^{2}>0$. Since the product of the roots is negative, we always have two solutions with opposite signs. It can be deduced that the solution is the positive root of the second-degree equation. For the cycle $k+1$ and the $i^{\mathrm{th}}$ coordinate, we have:

[TABLE]

In this equation, we have the following correspondence: $x_{i}\rightarrow x_{i}^{\left(k+1\right)}$, $x_{j}\rightarrow x_{j}^{\left(k+1\right)}$ if $j<i$, $x_{j}\rightarrow x_{j}^{\left(k\right)}$ if $j>i$, and $x\rightarrow\left(x_{1}^{\left(k+1\right)},\ldots,x_{i-1}^{\left(k+1\right)},x_{i}^{\left(k\right)},x_{i+1}^{\left(k\right)},\ldots,x_{n}^{\left(k\right)}\right)$. If the values of $\left(x_{1},\ldots,x_{n}\right)$ are strictly positive and if $c>\mathop{\rm SR}^{+}$, $x_{i}^{\left(k+1\right)}$ should be strictly positive. The positivity of the solution is then achieved after each iteration and each cycle if the starting values are positive. Therefore, the coordinate-wise descent algorithm consists in iterating the previous equation, and we can show that it always converges (Roncalli, 2015).

Remark 3

As noted by Griveau-Billion et al. (2013), the previous algorithm can be simplified by setting $\lambda$ equal to $1$ and by rescaling the solution once the convergence is obtained. Our experience shows that it is better to rescale the solution once the CCD algorithm has converged rather than after each cycle. Moreover, Griveau-Billion et al. (2013) derive analytical formulas in order to update $\sigma\left(x\right)$ and $\sum_{j\neq i}x_{j}\rho_{i,j}\sigma_{j}$ at each iteration.

3 Theory of constrained risk budgeting portfolio

3.1 Mathematical issues

If we consider the definition of Roncalli (2013), introducing constraints leads to the following formulation of the constrained risk budgeting portfolio333We assume that $b_{i}>0$ and $\sum_{i=1}^{n}b_{i}=1$.:

[TABLE]

where $\mathcal{S}$ is the standard simplex:

[TABLE]

and $x\in\Omega$ is the set of additional constraints. Let $x^{\star}\left(\mathcal{S}\right)$ be the risk budgeting portfolio, i.e. the solution such that $x\in\mathcal{S}$, and $x^{\star}\left(\mathcal{S},\Omega\right)$ be the constrained risk budgeting portfolio, i.e. the solution such that $x\in\mathcal{S}$ and $x\in\Omega$. Since $x^{\star}\left(\mathcal{S}\right)$ is unique, we deduce that the solution $x^{\star}\left(\mathcal{S},\Omega\right)$ exists only if $x^{\star}\left(\mathcal{S}\right)\in\Omega$, and we have $x^{\star}\left(\mathcal{S},\Omega\right)=x^{\star}\left(\mathcal{S}\right)$. Since we generally have $x^{\star}\left(\mathcal{S}\right)\notin\Omega$, we deduce that there is almost certainly no solution.

This is why professionals generally replace the equality constraint by an approximate equality:

[TABLE]

Therefore, the optimization problem becomes:

[TABLE]

Bai et al. (2016) show that this optimization problem can be simplified as follows:

[TABLE]

Example 1

We consider a universe of four assets. Their volatilities are equal to $10\%$, $15\%$, $20\%$ and $30\%$. The correlation matrix of asset returns is given by the following matrix:

[TABLE]

By using the volatility risk measure, we compute the ERC portfolio and the RB portfolio corresponding to the risk budgets $\left(30\%,30\%,19.5\%,20.5\%\right)$. In Table 1, we report the solution of ERC and RB portfolios. We also indicate the marginal risk $\mathcal{MR}_{i}$, the absolute risk contribution $\mathcal{RC}_{i}$ and the relative risk contribution $\mathcal{RC}_{i}^{\star}$. Let us now introduce the constraint $x_{i}\leq 30\%$ and solve the optimization problem (3.1). Since this constraint is not satisfied by the previous unconstrained portfolio, it has an impact as shown in Table 2. As expected, the relative risk contributions (or ex-post risk budgets) are completely different from the ex-ante risk budgets. The concept of “equal risk contribution” does not make any sense. Moreover, we observe that the ordering relationship between risk budgets are not preserved when we introduce constraints. For example, in the case of the RB portfolio, we have $b_{3}<b_{4}$ ($19.50\%$ versus $20.50\%$) but $\mathcal{RC}_{3}^{\star}>\mathcal{RC}_{4}^{\star}$ ($30.52\%$ versus $29.04\%$) although the allocation in Assets 3 and 4 does not reach the upper bound constraint. We also notice that the constrained ERC portfolio is very close to the constrained RB portfolio. This gives us the feeling that the choice of risk budgets has little impact, and the solution is mainly driven by the constraints.

3.2 Formulation of the optimization problem

Like for the unconstrained risk budgeting portfolio, we argue that the right optimization problem is:

[TABLE]

where $\kappa^{\star}$ is a constant to be determined. We notice that the previous problem can be simplified because:

the logarithmic barrier constraint imposes that $x_{i}\geq 0$; 2. 2.

Roncalli (2013) shows that there is only one constant $\kappa^{\star}$ such that the constraint $\mathbf{1}^{\top}x=1$ is satisfied.

We can then consider the following optimization problem:

[TABLE]

and we have:

[TABLE]

This new formulation is appealing since the constraint $x\in\mathcal{S}$ is not explicit, but it is implicitly embedded in the optimization problem. From a computational point of view, this reduces the complexity of the numerical algorithm. When the constraint $x\in\Omega$ vanishes, we retrieve the previous scaling rule. Otherwise, we consider the Lagrange formulation:

[TABLE]

Again, we have:

[TABLE]

Formulations (19), (23) and (3.2) have the advantage of revealing the true nature of risk budgeting. The objective is to minimize the risk measure subject to a penalization (Richard and Roncalli, 2015). A risk budgeting portfolio is then a minimum risk portfolio subject to hard risk budgeting, constraints, whereas a constrained risk budgeting portfolio is a minimum risk portfolio subject to soft risk budgeting constraints. Because of the convexity of the optimization problem (23), it follows that:

[TABLE]

for a given set of constraints $\Omega$. This property is fundamental since it is the essence of risk budgeting, and it explains the relationships between long-only minimum variance, risk budgeting and equally-weighted portfolios obtained by Maillard et al. (2010) and the relationships between long-only minimum risk, risk budgeting and weight-budgeting portfolios obtained by Roncalli (2013). This property is also necessary to impose the continuity of risk budgeting portfolios in particular when some risk budgets tend to zero. Without this property, it is impossible to show that the RB portfolio is unique, and to determine the true solution in the case where there are several solutions to Problem (2) when $b_{i}=0$.

Remark 4

*Imposing tighter constraints does not necessarily increase the risk measure444Indeed, we have:

$\Omega_{2}\subset\Omega_{1}\Rightarrow\mathcal{R}\left(x^{\star}\left(\Omega_{2},\kappa\right)\right)\geq\mathcal{R}\left(x^{\star}\left(\Omega_{1},\kappa\right)\right)$

for a given value of $\kappa$. However, the value $\kappa^{\star}$ to obtain the optimal portfolio $x^{\star}\left(\mathcal{S},\Omega\right)$ depends on $\Omega$.:*

[TABLE]

In order to illustrate the risk measure issue, we consider Example 1 with the risk budgets $b=\left(10\%,20\%,30\%,40\%\right)$. Moreover, we impose that the weight $x_{1}$ of Asset 1 is greater than a given lower bound $x_{1}^{-}$. In Figure 1, we have reported the portfolio volatility $\sigma\left(x^{\star}\left(\mathcal{S},\Omega\right)\right)$ of the constrained risk budgeting with respect to $x_{1}^{-}$. If we consider the least squares problem (3.1), the volatility is non-monotonous. This is not the case if we consider the logarithmic barrier problem (3.2). Indeed, the least squares formulation only considers the dimension of risk contribution matching, but not the dimension of risk measure minimization.

3.3 Numerical algorithms

The optimization function becomes:

[TABLE]

where $\mathds{1}_{\Omega}\left(x\right)$ is the convex indicator function of $\Omega$, meaning that $\mathds{1}_{\Omega}\left(x\right)=0$ for $x\in\Omega$ and $\mathbf{1}_{\Omega}\left(x\right)=+\infty$ for $x\notin\Omega$. The choice of $\lambda$ is very important, since the constrained RB portfolio is obtained for the optimal value $\lambda^{\star}$ such that the sum of weights is equal to one. This can be done using the Newton-Raphson or the bisection algorithm. If we note $x^{\star}\left(\lambda\right)$ the solution of the minimization problem (29), we obtain Algorithm 1 in the case of the bisection method.

3.3.1 ADMM algorithm

In order to solve Problem (29), we exploit the separability of $\mathcal{L}\left(x;\lambda\right)$. For example, we can write:

[TABLE]

or:

[TABLE]

We notice that we have:

[TABLE]

It follows that we can use the alternative direction method of multipliers (ADMM) to solve this optimization problem. Algorithm 2 describes the different steps.

In the case of the Lagrange function (30), the $x$-update is equivalent to solving a penalized risk budgeting problem whereas the $z$-update corresponds to a proximal operator. Therefore, we can use the Newton algorithm555The first and second derivatives of $f^{\left(k\right)}\left(x\right)=f\left(x\right)+\dfrac{\varphi}{2}\left\|x-z^{\left(k-1\right)}+u^{\left(k-1\right)}\right\|_{2}^{2}$ are equal to:

$\displaystyle\frac{\partial\,f^{\left(k\right)}\left(x\right)}{\partial\,x_{i}}$ $\displaystyle=$ $\displaystyle\frac{\partial\,f\left(x\right)}{\partial\,x_{i}}+\varphi\left(x_{i}-z^{\left(k-1\right)}+u^{\left(k-1\right)}\right)$

$\displaystyle\frac{\partial^{2}\,f^{\left(k\right)}\left(x\right)}{\partial\,x_{i}\,\partial\,x_{j}}$ $\displaystyle=$ $\displaystyle\frac{\partial^{2}\,f\left(x\right)}{\partial\,x_{i}\,\partial\,x_{j}}$

$\displaystyle\frac{\partial^{2}\,f^{\left(k\right)}\left(x\right)}{\partial\,x_{i}^{2}}$ $\displaystyle=$ $\displaystyle\frac{\partial^{2}\,f\left(x\right)}{\partial\,x_{i}^{2}}+\varphi$

where the derivatives of $f\left(x\right)$ are given in Footnote 2 on page 2. to find $x^{\left(k\right)}$. For the $z$-update, we have:

[TABLE]

where $v_{z}^{\left(k\right)}=x^{\left(k\right)}+u^{\left(k-1\right)}$. If we assume that $g\left(z\right)=\mathds{1}_{\Omega}\left(z\right)$ where $\Omega$ is a convex set, we obtain:

[TABLE]

where $\mathcal{P}_{\Omega}\left(v\right)$ is the standard projection. In Appendix A.2 on page A.2, we give the results for the generic constraints that we encounter in portfolio optimization. We develop some special cases in the next section.

For the Lagrange function (31), the $x$-update corresponds to a constrained risk minimization problem whereas the $z$-update is equivalent to solving a penalized logarithmic barrier problem. The $x$-update can be done using constrained non-linear optimization methods666In order to accelerate the convergence, we can implement analytical derivatives, which are the same than those given in Footnote 5 on page 5 by setting $\lambda=0$ in the derivatives of the function $f\left(x\right)$., whereas the $z$-step corresponds to the proximal operator of the logarithmic barrier function777See Appendix A.4 on page A.4.:

[TABLE]

If we consider the volatility risk measure $\mathcal{R}\left(x\right)=\sqrt{x^{\top}\Sigma x}$ instead of the standard deviation-based risk measure given by Equation (1), the ADMM algorithm is simplified as follows888We have:

$\displaystyle f^{\left(k\right)}\left(x\right)$ $\displaystyle=$ $\displaystyle\frac{1}{2}x^{\top}\Sigma x+\frac{\varphi}{2}\left(x-z^{\left(k-1\right)}+u^{\left(k-1\right)}\right)^{\top}\left(x-z^{\left(k-1\right)}+u^{\left(k-1\right)}\right)$

$\displaystyle=$ $\displaystyle\frac{1}{2}x^{\top}\Sigma x+\frac{\varphi}{2}\left(x^{\top}x-2x^{\top}v_{x}^{\left(k\right)}+\left(v_{x}^{\left(k\right)}\right)^{\top}v_{x}^{\left(k\right)}\right)$

where $v_{x}^{\left(k\right)}=z^{\left(k-1\right)}-u^{\left(k-1\right)}$.:

[TABLE]

This algorithm exploits the property that minimizing the portfolio volatility is equivalent to minimizing the portfolio variance, even if this last risk measure does not satisfy the Euler decomposition. If $\Omega$ is a set of linear (equality and inequality) constraints, the $x$-update reduces to a standard QP problem.

3.3.2 CCD algorithm

Another route for solving Problem (30) is to consider the CCD algorithm. Convergence of coordinate descent methods requires that the function is strictly convex and differentiable. However, Tseng (2001) has extended the convergence properties to a non-differentiable class of functions:

[TABLE]

where $f_{0}$ is strictly convex and differentiable and the functions $f_{i}$ are non-differentiable. Dealing with convex constraints is equivalent to writing the constraints in the sum term of $f\left(x\right)$. If we consider the formulation (30) and the standard deviation-based risk measure, we have:

[TABLE]

where $\mathcal{L}_{0}\left(x;\lambda\right)$ is defined as follows:

[TABLE]

The case of separable constraints

If we assume that the set of constraints is separable with respect to all the variables:

[TABLE]

where $\Omega_{i}$ is the constraint on $x_{i}$, we have $f_{i}\left(x_{i}\right)=\mathds{1}_{\Omega_{i}}\left(x_{i}\right)$. We deduce that the CCD algorithm consists in two steps. We first solve the minimization problem of $\mathcal{L}_{0}\left(x;\lambda\right)$ for one coordinate, and then we compute the projection onto $\Omega_{i}$. For the first step, the first-order condition is:

[TABLE]

It follows that $cx_{i}\left(\Sigma x\right)_{i}-\pi_{i}x_{i}\sigma\left(x\right)-\lambda b_{i}\sigma\left(x\right)=0$ or equivalently:

[TABLE]

where:

[TABLE]

We deduce that the coordinate solution is the positive root of the second-degree equation:

[TABLE]

The second step is the projection into the set $\Omega_{i}$:

[TABLE]

Finally, we obtain Algorithm 3.

Remark 5

Our algorithm differs from the one given by Nesterov (2012) and Wright (2015). Let $\eta>0$ be the stepsize of the gradient descent. The coordinate update is:

[TABLE]

where $\xi$ is a positive scalar and:

[TABLE]

In our case, this algorithm is very simple999See Appendix A.5 on page A.5., because it reduces to calculate the proximal of $x_{i}-\eta g_{i}$ associated with the function $\mathds{1}_{\Omega_{i}}\left(x\right)$. However, we prefer to use the previous algorithm in order to exploit the analyticity of the solution.

The case of non-separable constraints

A first idea is to replace the projection step $x_{i}\leftarrow\mathcal{P}_{\Omega_{i}}\left(x_{i}\right)$ by something equivalent that ensures that the constraints are verified. The natural approach is to apply the proximal operator or equivalently the projection: $x\leftarrow\mathcal{P}_{\Omega}\left(x\right)$. In practice, we observe that the CCD solution does not always converge to the true solution. It will depend on how the constraints and the variables are ordered. In fact, the true approach is to use a block-coordinate algorithm if constraints are separable by blocks. This is not always the case. This is why we prefer to use the previous ADMM algorithms or the ADMM-CCD algorithm, which is described in the next paragraph.

3.3.3 Mixed ADMM-CCD algorithm

If we consider the formulation (30), we have already shown that the $x$-update of the ADMM algorithm corresponds to a regularized risk budgeting problem. Therefore, we can use the CCD algorithm to find the solution $x^{\left(k\right)}$. We remind that:

[TABLE]

where $v_{x}^{\left(k\right)}=z^{\left(k-1\right)}-u^{\left(k-1\right)}$. The first-order condition is:

[TABLE]

It follows that:

[TABLE]

or:

[TABLE]

where:

[TABLE]

We deduce that the coordinate solution is the positive root of the previous second-degree equation:

[TABLE]

Finally, we obtain the ADMM-CCD algorithm, which is described on the previous page.

3.3.4 Efficiency of the algorithms

We may investigate the efficiency of the previous algorithms. Firstly, our experience shows that traditional constrained optimization algorithms (SQP, trust region, constrained interior point, etc.) fail to find the solution. It is somewhat surprising because we have the feeling that the optimization problem of constrained risk budgeting portfolios seems to be standard. Unfortunately, this is not the case, because the mixing of the logarithmic barrier and constraints is not usual. This is why it is important to implement the previous algorithms. Secondly, all the algorithms are not equal and the implementation is key in particular when we consider large problems with more than one hundred assets. In Table 3, we have reported the computational time we have obtained for solving the example described on page 3.4.1. For that, we consider five different methods: ADMM-Newton, ADMM-BFGS, ADMM-QP, ADMM-CCD and CCD. For each method, we consider three implementations:

The first one considers that the primal variable $\varphi$ is constant ($\varphi=1$) and we use the classical bisection method described in Algorithm 1. 2. 2.

The second one considers that the penalization variable $\varphi$ is constant ($\varphi=1$), and we use an accelerated bisection method. The underlying idea is to choose a starting value $x^{\left(0\right)}$ for the ADMM/CCD algorithm, which is not constant and depends on the Lagrange coefficient $\lambda$. The corresponding method is described in Appendix A.6 on page A.6. 3. 3.

The third uses the accelerated bisection algorithm, and considers the adaptive method for the variable $\varphi^{\left(k\right)}$, which is given in Appendix A.7 on page A.7. The underlying idea is to accelerate the convergence of the ADMM algorithm by using the right scale of the primal residual variable $u^{\left(k\right)}$.

Results using our Matlab implementation are reported in Table 3. Absolute figures are not interesting, because they depend on the processing power of the computer. Relative figures show that computational times can differ dramatically from one algorithm to another, from one implementation to another. The best algorithms are ADMM-Newton and CCD, followed by ADMM-CCD. Curiously, the ADMM-QP method is less efficient101010This is due to the quadprog procedure of Matlab. Indeed, we do not observe same results in Python.. Finally, the worst algorithm is the ADMM-BFGS algorithm. However, we notice a large improvement in this algorithm if we implement the accelerated bisection and the adaptive method for scaling the regularization parameter. Indeed, the computational time is divided by a factor of 15!

3.4 Special cases

3.4.1 Box constrained optimization

In portfolio optimization, imposing lower and upper bounds is frequent:

[TABLE]

For example, box constraints are used when limiting single exposures because of regulatory constraints or controlling the turnover of the portfolio.

The optimization framework

In the box constrained case, the Lagrange function becomes:

[TABLE]

The first-order condition is:

[TABLE]

We deduce that:

[TABLE]

Since the Kuhn-Tucker conditions are:

[TABLE]

we obtain three cases:

If no bound is reached, we have $\lambda_{i}^{-}=0$ and $\lambda_{i}^{+}=0$, and we retrieve the RB portfolio $x_{\mathrm{RB}}$; 2. 2.

If the lower bound is reached, $\lambda_{i}^{-}>0$ and the risk contribution of Asset $i$ is higher than $b_{i}$; 3. 3.

If the upper bound is reached, $\lambda_{i}^{+}>0$ and the risk contribution of Asset $i$ is lower than $b_{i}$.

Remark 6

When we compare the logarithmic barrier solution with the least squares solution of the ERC portfolio, we observe that the first one preserves the “equal risk contribution” property for all the assets that do not reach lower or upper bounds. This is not the case with the least squares solution, since the risk contribution is different for all the assets.

The previous algorithms require the proximal operator to be computed:

[TABLE]

where $\tilde{x}=x^{\left(k\right)}+u^{\left(k-1\right)}$ is the value of $v_{z}^{k}$ in the ADMM procedure. In Appendix A.2 on page A.2, we show that the proximal operator corresponds to the truncation operator:

[TABLE]

where:

[TABLE]

For the CCD algorithm, the projection $x_{i}\leftarrow\mathcal{P}_{\Omega_{i}}\left(x_{i}\right)$ reduces to apply the truncation operator to the single coordinate $x_{i}$.

At the optimum, we deduce that:

[TABLE]

and:

[TABLE]

where $\lambda^{\star}$ is the solution of the bisection algorithm.

Remark 7

*In order to find the optimal value $\lambda^{\star}$, we need an initial guess $\lambda_{0}$ for the Newton-Raphson algorithm or an interval $\left[a_{\lambda},b_{\lambda}\right]$ for the bisection method. Let $x_{\mathrm{RB}}$ be the RB portfolio without constraints. The optimal Lagrange coefficient associated with $x_{\mathrm{RB}}$ is equal to111111Indeed, the first-order condition is:

$\displaystyle\frac{\partial\,\mathcal{R}\left(x\right)}{\partial\,x_{i}}-\lambda\frac{b_{i}}{x_{i}}=0$ $\displaystyle\Leftrightarrow$ $\displaystyle x_{i}\frac{\partial\,\mathcal{R}\left(x\right)}{\partial\,x_{i}}-\lambda b_{i}=0$

$\displaystyle\Leftrightarrow$ $\displaystyle\sum_{i=1}^{n}x_{i}\frac{\partial\,\mathcal{R}\left(x\right)}{\partial\,x_{i}}-\lambda\sum_{i=1}^{n}b_{i}=0$

$\displaystyle\Leftrightarrow$ $\displaystyle\lambda=\sum_{i=1}^{n}x_{i}\frac{\partial\,\mathcal{R}\left(x\right)}{\partial\,x_{i}}=\mathcal{R}\left(x\right)$

because $\sum_{i=1}^{n}b_{i}=1$.:*

[TABLE]

It follows that a good initial guess is $\lambda_{0}=\mathcal{R}\left(x_{\mathrm{RB}}\right)$. We also notice that:

[TABLE]

We deduce that $a_{\lambda}=m_{a}\cdot\mathcal{R}\left(x_{\mathrm{RB}}\right)$ and $b_{\lambda}=m_{b}\cdot\mathcal{R}\left(x_{\mathrm{RB}}\right)$ where the parameters $m_{a}$ and $m_{b}$ depends on the tightness of constraints121212We have $m_{a}\leq 1$ and $m_{b}\geq 1$.. Generally, we have $\mathcal{R}\left(x^{\star}\left(S,\Omega\right)\right)\approx\mathcal{R}\left(x_{\mathrm{RB}}\right)$, implying that the values $m_{a}=0.5$ and $m_{b}=2.0$ are sufficient.

An example of dynamic allocation

We consider a universe of five assets. Their volatilities are equal to $15\%$, $20\%$, $25\%$, $30\%$ and $10\%$. The correlation matrix of asset returns is given by the following matrix:

[TABLE]

Let us assume that the current portfolio is $x_{0}=\left(25\%,25\%,10\%,15\%,30\%\right)$. In Table 4, we report the volatility breakdown of this portfolio. We notice that there are some large differences in terms of risk contributions. In particular, the second asset has a volatility contribution of $31.1\%$. We would like to obtain a more balanced portfolio. Table 5 shows the results of the ERC portfolio.

We notice that the ERC portfolio is relatively far from the current portfolio. In particular, the turnover is equal to $22.18\%$. In order to obtain a solution closer to the current allocation, we impose that the weights cannot deviate from the current ones by $5\%$:

[TABLE]

The underlying idea is to move from the initial portfolio to a risk budgeting portfolio, which presents the risk parity property as much as possible. In this case, we obtain the results in Table 6. We notice that three assets (#1, #3 and #4) present the same contributions ($2.35\%$) because they do not reach the bounds. On the contrary, the second and fifth assets have respectively a higher and lower risk contribution ($2.98\%$ and $2.10\%$) because the lower and upper bounds are reached. This solution helps to reduce the turnover, since it is equal to $14.22\%$.

A naive solution to obtain a risk parity portfolio that matches the constraints would be to identify the assets that reach the lower and upper bounds and to allocate the remaining weight between the other assets by imposing the same risk contribution. In this example, the second and fifth assets do not satisfy the constraints. Therefore, we have to allocate $45\%$ of the allocation between the first, third and fourth assets. This naive solution is given in Table 7. The ERC portfolio between the three unconstrained assets is equal to $\left(22.84\%,12.34\%,9.83\%\right)$. In this case, the risk contributions are the same and are equal to $2.65\%$. However, the equal risk contribution property does not hold if we consider the full portfolio, when we take into account the constrained assets. Indeed, we obtain $\mathcal{RC}_{1}=2.34\%$, $\mathcal{RC}_{3}=2.49\%$ and $\mathcal{RC}_{4}=2.21\%$. The reason is that risk budgeting portfolios are sensitive to the asset universe definition (Roncalli and Weisang, 2016). This is why this two-step naive procedure does not give the right answer. In Table 7, we also report the least squares solution corresponding to the optimization problem (3.1). Again, no assets verify the equal risk contribution property.

Remark 8

The previous example shows how to take into account a current allocation when building a risk budgeting portfolio. This approach is particularly relevant when considering dynamic rebalancing and tactical asset allocation. One of the main advantages of the RB approach is that it produces a stable allocation. However, it does not enable us to consider investment constraints such as the current allocation. Introducing weight constraints allows the fund manager to better control the portfolio construction.

3.4.2 Risk budgeting with linear constraints

We now consider general linear constraints:

[TABLE]

These constraints generalize the case of lower and upper bounds. For example, it is common to add some (lower and upper) limits in terms of exposures by asset classes, sectors, ratings, etc. These types of limits are implemented thanks to inequality constraints $Cx\leq D$. Equality constraints $Ax=B$ are less common in portfolio optimization. In Appendix A.2 on page A.2, we provide some closed-form formulas to calculate $\mathcal{P}_{\Omega}\left(x\right)$, when $\Omega$ corresponds to $Ax=B$ or $c^{\top}x\leq d$ or $x^{-}\leq x\leq x^{+}$. First, we notice that the analytical formula only exists for the half-space constraint $c^{\top}x\leq d$, but not for multiple inequality constraints $\Omega=\left\{x\in\mathbb{R}^{n}:Cx\leq D\right\}$ when $\mathop{\rm card}\Omega=m>1$. The underlying idea is then to break down $\Omega$ as the intersection of $m$ half-space sets:

[TABLE]

where $\Omega_{j}=\left\{x\in\mathbb{R}^{n}:c_{\left(j\right)}^{\top}x\leq d_{\left(j\right)}\right\}$, $c_{\left(j\right)}^{\top}$ corresponds to the $j^{\mathrm{th}}$ row of $C$ and $d_{\left(j\right)}$ is the $j^{\mathrm{th}}$ element of $D$. In this case, we can apply the Dykstra’s algorithm for computing the proximal operator of $\mathds{1}_{\Omega}\left(x\right)$. This algorithm is given on page 5. Second, mixing the constraints is not straightforward. Again, we can break down $\Omega$ as the intersection of three basic convex sets:

[TABLE]

where $\Omega_{1}=\left\{\in\mathbb{R}^{n}:Ax=B\right\}$, $\Omega_{2}=\left\{x\in\mathbb{R}^{n}:Cx\leq D\right\}$ and $\Omega_{3}=\left\{x\in\mathbb{R}^{n}:x^{-}\leq x\leq x^{+}\right\}$. Since we know how to project each basic set, we calculate the projection $\mathcal{P}_{\Omega}\left(x\right)$ with the Dykstra’s algorithm, which is described on page 6.

Let us consider an example of multi-asset allocation131313This example is taken from Roncalli (2013) on page 287.. We consider a universe of eight asset classes: (1) US 10Y Bonds, (2) Euro 10Y Bonds, (3) Investment Grade Bonds, (4) High Yield Bonds, (5) US Equities, (6) Euro Equities, (7) Japan Equities and (8) EM Equities. In Table 8, we indicate the statistics used to compute the optimal allocation.

Using these figures, we calculate the risk parity portfolio, which corresponds to the second column in Table 9. We notice that the bond allocation is equal to $76.72\%$ whereas the equity allocation is equal to $23.28\%$. In order to increase the equity allocation, we impose that the weight of the last four assets is greater than $30\%$. The solution is given in the fourth column. Finally, we overweight the allocation in European assets with respect to American assets by $5\%$ (sixth column).

Remark 9

Algorithm 3 is no longer valid when the coordinates are coupled via constraints. This is generally the case when we impose $Ax=B$ and $Cx\leq D$. This is why we use the ADMM-Newton or ADMM-CCD algorithms for solving this type of constrained risk budgeting problem.

4 Applications

We consider two applications that are based on our professional experience. The first application is the design of risk-based equity indices. Building an ERC portfolio on the Eurostoxx 50 universe is straightforward. This is not the case if we consider the universe of the Eurostoxx index. Indeed, this universe contains many small cap stocks, and having the same risk contribution for small cap and large cap stocks may induce some liquidity issues. In particular, this type of problem occurs when considering a large universe of non-homogenous stocks. The second application is the control of rebalancing effects. This issue happens when we consider short-term covariance matrices. In this case, the allocation can be very reactive, implying large turnovers. For example, this type of situation is observed when we implement multi-asset risk parity strategies with daily or weekly rebalancing and empirical covariance matrices that are estimated with less than one year of historical data. These two applications are illustrated below.

4.1 Risk-based indexation and smart beta portfolios

We now consider a capitalization-weighted index composed of seven stocks. The weights are equal to $34\%$, $25\%$, $20\%$, $15\%$, $3\%$, $2\%$ and $1\%$. We assume that the volatilities of these stocks are equal to $15\%$, $16\%$, $17\%$, $18\%$, $19\%$, $20\%$ and $21\%$, whereas the correlation matrix of stock returns is given by:

[TABLE]

As shown by Demey et al. (2010), the ERC portfolio defined by Maillard et al (2010) is a good candidate for building a risk-based equity index. However, an ERC index does not take into account liquidity constraints. For instance, we notice that the ERC allocation does not take into account the size of stocks in Table 10. We may assume that the three last assets are small cap stocks. In this case, it can be more realistic to distinguish small cap and large cap stocks. A first idea is to keep the CW weights on the small cap universe and to apply the ERC on the large cap universe. This solution called LC-ERC (for Large Cap ERC) is presented in Table 10. We face an issue here, because the ERC portfolio on large cap stocks does not depend on the full correlation matrix. Therefore, this solution assumes that the two universes of stocks are not correlated. We have the same issue if we consider the least squares solution (LS-ERC). A better approach is to find the ERC portfolio by imposing that the weights of small cap stocks are exactly equal to the corresponding CW weights. Let $\Omega_{\mathcal{SC}}$ be the universe of small cap stocks. We have $x_{i}=x_{\mathrm{cw},i}$ if $i\in\Omega_{\mathcal{SC}}$. This is equivalent to imposing the following weight constraints:

[TABLE]

The result corresponds to the C-ERC portfolio. Again, we notice that the property of equal risk contribution is satisfied at the global level for large cap stocks, contrary to LC-ERC and LS-ERC portfolios.

4.2 Managing the portfolio turnover

On page 3.4.1, we have considered a dynamic allocation and imposed some constraints in order to control the turnover of the portfolio. For that, we have used lower and upper bounds in order to impose a maximum deviation between the current allocation and the new allocation. This is a way of controlling the turnover. However, we can directly impose a turnover control. Let $x_{0}$ be the current allocation. The two-way turnover of Portfolio $x$ with respect to Portfolio $x_{0}$ is defined by:

[TABLE]

It corresponds to the $\boldsymbol{\ell}_{1}$-norm of $x$ with respect to the centroid vector $x_{0}$. Therefore, the corresponding Lagrange function is:

[TABLE]

where $\Omega=\left\{x\in R:\tau\left(x\mid x_{0}\right)\leq\tau^{\star}\right\}$ and $\tau^{\star}$ is the turnover limit. If we use the previous algorithms, the only difficulty is calculating the proximal operator141414We use the properties of proximal operators described in Appendix A.2 on page A.2. of $g\left(x\right)=\mathds{1}_{\Omega}\left(x\right)$:

[TABLE]

where $f\left(x\right)=\mathds{1}_{\Omega^{\prime}}\left(x\right)$ and $\Omega^{\prime}=\left\{x\in R:\left\|x\right\|_{1}\leq\tau^{\star}\right\}$. Finally, we deduce that:

[TABLE]

where $\mathbf{prox}_{\lambda\max}\left(v\right)$ is the proximal operator given by Equation (41) on page 41.

We consider the example of multi-asset allocation on page 8. Let us assume that the current allocation is a 50/50 asset mix policy, where the weight of each asset class is $12.5\%$. In Table 11, we have reported the solution for different turnover limits $\tau^{\star}$. If the turnover limit is very low, the optimized RB portfolio is close to the current allocation. We verify that the optimized portfolio tends to the ERC portfolio when we increase the turnover limit.

5 Discussion and limitations of constrained risk budgeting portfolios

In this section, we discuss the concept of constrained risk budgeting allocation, and shows that it is not natural and has some limitations.

5.1 Coherent risk measures and the homogeneity property

Many people believe that the risk budgeting allocation is only related to the Euler decomposition:

[TABLE]

implying that the risk measure is “convex”. However, this concept is not always well-defined, and is often confused with the concept of coherent risk measure. Following Artzner et al. (1999), a risk measure $\mathcal{R}\left(x\right)$ is said to be coherent if it satisfies the following properties:

Subadditivity

[TABLE]

The risk of two portfolios should be less than adding the risk of the two separate portfolios. 2. 2.

Homogeneity

[TABLE]

Leveraging or deleveraging of the portfolio increases or decreases the risk measure in the same magnitude. 3. 3.

Monotonicity

[TABLE]

If Portfolio $x_{2}$ has a better return than Portfolio $x_{1}$ under all scenarios, risk measure $\mathcal{R}\left(x_{1}\right)$ should be higher than risk measure $\mathcal{R}\left(x_{2}\right)$. 4. 4.

Translation invariance

[TABLE]

Adding a cash position of amount $m$ to the portfolio reduces the risk by $m$.

Föllmer and Schied (2002) propose replacing the homogeneity and subadditivity conditions by a weaker condition called the convexity property:

[TABLE]

This condition means that diversification should not increase the risk. Saying that the risk measure is convex is ambiguous. For some people, this means that $\mathcal{R}\left(x\right)$ satisfies the convexity property of Föllmer and Schied (2002), while for other people, this means that $\mathcal{R}\left(x\right)$ satisfies the Euler decomposition. It is true that these two concepts are related (Tasche, 2008), but they recover two different things (Kalkbrener, 2005). First, there is a confusion between the Euler decomposition and the Euler allocation principle (Tasche, 2008). Second, the Euler allocation principle does make sense only if the risk measure is subadditive (Kalkbrener, 2005). But another property is very important. Roncalli (2015) shows that risk budgeting is valid only if the homogeneity property is satisfied — $\mathcal{R}\left(\lambda x\right)=\lambda\mathcal{R}\left(x\right)$ — because this property ensures that there is a solution and the solution is unique. By definition, this property is related to the scaling property of the RB portfolio when there are no constraints. When we impose some constraints, it is obvious that the homogeneity property is valid only if the constraints $\Omega$ are compatible with the scaling property. This is not generally true except for some special cases. However, we don’t need the homogeneity property to be satisfied for all values $\lambda\geq 0$. We need the homogeneity property to be met for a range of $\lambda$ around the unconstrained risk budgeting portfolio. This is why imposing tight constraints can lead to a numerical solution without knowing if it corresponds to the true constrained risk budgeting portfolio.

5.2 The scaling puzzle

We consider the example given on page 8. If we consider the optimized portfolio subject to the constraint $\sum_{i=5}^{8}x_{i}\geq 30\%$, we may think that it is equivalent to the optimized portfolio subject to the constraint $\sum_{i=1}^{4}x_{i}\leq 70\%$. Results are given in Table 12, when we use the equally-weighted portfolio as the starting value $x^{\left(0\right)}$ in the ADMM algorithms. It is surprising to obtain two different solutions. Nevertheless, they are very close. If we do the same exercise with a minimum allocation of 40% in the equity asset class, the two solutions are very different (see columns 6 and 8 in Table 12). The problem is that we assume that the constraint $\sum_{i=1}^{n}x_{i}=1$ is managed by the set $\Omega$. This is not the case, because the constraint $\sum_{i=1}^{n}x_{i}=1$ is managed by the Lagrange multiplier $\lambda$ associated to the logarithmic barrier. Here, we face an important issue called the scaling compatibility problem. This means that a solution is acceptable if and only if the constraints $\Omega$ are “compatible” with the homogeneity property151515For example, $\Omega=\left\{x\in\mathbb{R}^{n}:x_{2}\geq 2x_{1}\right\}$ is compatible with the scaling property..

How do we explain these results? The first reason is the choice of the starting value for initializing the algorithm. It is obvious that the bisection algorithm takes a road that depends on the initialization step. However, this reason is not the primary answer, because we observe that we converge to the same solution whatever the starting value or the algorithm when we consider box or pointwise constraints. In fact, the main reason is the scaling property. In this case, the convergence path is crucial, because we can obtain local minima.

When the constraints are incompatible with the homogeneity property of the risk measure $\mathcal{R}\left(x\right)$, we may wonder if one solution is better than the others. In order to answer this question, we recall that the numerical algorithms solve the minimization problem $x^{\star}\left(\lambda\right)=\arg\min\mathcal{L}\left(x;\lambda\right)$ where:

[TABLE]

The idea is then to calculate $\mathcal{L}\left(x^{\star}\left(\lambda^{\star}\right);\lambda^{\star}\right)$ for the different solutions and to take the solution that gives the lowest value. In Table 12, the first and third portfolios are the best solutions. They have the lowest volatility and Lagrange function.

Remark 10

Let us assume that the set of constraints includes the standard simplex: $\mathcal{S}\subset\Omega$. By construction, the sum of weights is always equal to $1$ whatever the value of the Lagrange multiplier $\lambda$:

[TABLE]

It gives the impression that there are two solutions. However, there is only one solution which corresponds to the portfolio with the lowest risk measure.

6 Conclusion

In this paper161616This aim of this paper is also to help to popularize large-scale optimization algorithms that are very popular in machine learning, but are not well-known in finance. However, we think that they are relevant for many financial applications, in particular for portfolio optimization. Our paper which is focused on risk-budgeting optimization can then be seen as a companion work of Bourgeron et al. (2018) which focuses on mean-variance optimization., we propose an approach to find the risk budgeting portfolio when we impose some constraints. The underlying idea is to consider the logarithmic barrier problem and to use recent optimization algorithms in order to find the numerical solution. In particular, we use cyclical coordinate descent (CCD), alternative direction method of multipliers (ADMM), proximal operators and Dykstra’s algorithm.

We provide different examples that are focused on the ERC portfolio. This portfolio is very interesting since it imposes the same risk contribution between the assets of the investment universe. We may then wonder what does an ERC portfolio mean when we impose some constraints. Most of the times, we observe that the ERC property continues to be satisfied for the assets that are not impacted by the constraints. This type of approach is very appealing when we would like to manage the liquidity of a portfolio, the small cap bias of risk-based indices or the turnover of a risk parity fund.

This study also highlights the importance of the homogeneity property of RB portfolios. Roncalli (2015) has already pointed out that the risk measure must be coherent, which can be incompatible when imposing constraints. The homogeneity property really asks the question of the compatibility between risk budgeting allocation and imposing some constraints. However, even if the calibration of the Lagrange multiplier $\lambda$ remains an issue in some cases, our approach extends the optimization framework defined by Richard and Roncalli (2015), and reinforces that idea that risk budgeting and risk minimization are highly connected.

Appendix

Appendix A Optimization algorithms

A.1 ADMM algorithm

The alternating direction method of multipliers (ADMM) is an algorithm introduced by Gabay and Mercier (1976) to solve problems which can be expressed as171717We follow the standard presentation of Boyd et al. (2011) on ADMM.:

[TABLE]

where $A\in\mathbb{R}^{p\times n}$, $B\in\mathbb{R}^{p\times m}$, $c\in\mathbb{R}^{p}$, and the functions $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{+\infty\}$ and $g:\mathbb{R}^{m}\rightarrow\mathbb{R}\cup\{+\infty\}$ are proper closed convex functions. Boyd et al. (2011) show that the ADMM algorithm consists of three steps:

The $x$-update is:

[TABLE] 2. 2.

The $z$-update is:

[TABLE] 3. 3.

The $u$-update is:

[TABLE]

In this approach, $u^{\left(k\right)}$ is the dual variable of the primal residual $r=Ax+Bz-c$ and $\varphi$ is the $\boldsymbol{\ell}_{2}$ penalty variable. In the paper, we use the notations $f^{\left(k\right)}\left(x\right)$ and $g^{\left(k\right)}\left(z\right)$ when referring to the objective functions that are defined in the $x$- and $z$-steps.

A.2 Proximal operator

In what follows, we give the main results that are summarized in Bourgeron et al. (2017). Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\left\{+\infty\right\}$ be a proper closed convex function. The proximal operator $\mathbf{prox}_{f}\left(v\right):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is defined by:

[TABLE]

Since the function $f_{v}\left(x\right)=f\left(x\right)+\dfrac{1}{2}\left\|x-v\right\|_{2}^{2}$ is strongly convex, it has a unique minimum for every $v\in\mathbb{R}^{n}$ (Parikh and Boyd, 2014). For example, if we consider the logarithmic barrier function $f\left(x\right)=-\ln x$, we have:

[TABLE]

The first-order condition is $-x^{-1}+2x-v=0$. We obtain two roots with opposite signs. Since the logarithmic function is defined for $x>0$, we deduce that the proximal operator is:

[TABLE]

More generally, if we consider $f\left(x\right)=-\lambda\sum_{i=1}^{n}\ln x_{i}$, we have:

[TABLE]

Let us now consider some special cases. If we assume that $f\left(x\right)=\mathds{1}_{\Omega}\left(x\right)$ where $\Omega$ is a convex set, we have:

[TABLE]

where $\mathcal{P}_{\Omega}\left(v\right)$ is the standard projection. Here, we give the results of Parikh and Boyd (2014) for some simple polyhedra:

[TABLE]

where $A^{\dagger}$ is the Moore-Penrose pseudo-inverse of $A$, and $\mathcal{T}\left(v;x^{-},x^{+}\right)$ is the truncation operator:

[TABLE]

In the case of complex polyhedra, the reader can find analytical formulas and numerical algorithms in Parikh and Boyd (2014), and Combettes and Pesquet (2011).

We also have:

[TABLE]

If we define $\Omega$ as follows:

[TABLE]

we obtain:

[TABLE]

Imposing linear constraints is then equivalent to solving a standard QP problem.

We now consider the case of norm functions. For that, we need a preliminary result. In the case of the pointwise maximum function $f\left(x\right)=\max x$, we have:

[TABLE]

where $s^{\star}$ is the solution of the following equation:

[TABLE]

If we assume that $f\left(x\right)=\left\|x\right\|_{p}$, we obtain181818The proximal operator $S_{\lambda}\left(v\right)$ is known as the soft thresholding operator. In particular, it is used for solving lasso regression problems (Friedman et al., 2010).:

[TABLE]

An important property of the proximal operator is the Moreau decomposition theorem:

[TABLE]

where $f^{\ast}$ is the convex conjugate of $f$. If $f\left(x\right)$ is a $\boldsymbol{\ell}_{p}$-norm function, then $f^{\ast}\left(x\right)=\mathds{1}_{\mathcal{B}_{p}}\left(x\right)$ where $\mathcal{B}_{p}$ is the $\boldsymbol{\ell}_{p}$ unit ball. Since we have $\mathbf{prox}_{f^{\ast}}\left(v\right)=\mathcal{P}_{\mathcal{B}_{p}}\left(v\right)$, we deduce that:

[TABLE]

More generally, we have:

[TABLE]

It follows that the projection on $\boldsymbol{\ell}_{p}$ ball can be deduced from the proximal operator of the $\boldsymbol{\ell}_{p}$-norm function. Let $\mathcal{B}_{p}\left(c,\lambda\right)=\left\{x\in\mathbb{R}^{n}:\left\|x-c\right\|_{p}\leq\lambda\right\}$ be the $\boldsymbol{\ell}_{p}$ ball with center $c$ and radius $\lambda$. We obtain:

[TABLE]

In the case where the center $c$ is not equal to $\mathbf{0}$, we consider the translation property:

[TABLE]

where $g\left(x\right)=f\left(x+c\right)$. Since we have the equivalence $\mathcal{B}_{p}\left(\mathbf{0},\lambda\right)=\left\{x\in\mathbb{R}^{n}:f\left(x\right)\leq\lambda\right\}$ where $f\left(x\right)=\left\|x\right\|_{p}$, we deduce that:

[TABLE]

A.3 Dykstra’s algorithm

We consider the following proximal problem:

[TABLE]

where $f\left(x\right)=\mathds{1}_{\Omega}\left(x\right)$ and:

[TABLE]

The solution can be found thanks to Dykstra’s algorithm (Dykstra, 1983; Bauschke and Borwein, 1994), which consists in the following two steps until convergence:

The $x$-update is:

[TABLE] 2. 2.

The $z$-update is:

[TABLE]

where $x^{\left(0\right)}=v$, $z^{\left(k\right)}=\mathbf{0}$ for $k<0$ and $\mathop{\rm mod}\left(k,m\right)$ denotes the modulo operator taking values in $\left\{1,\ldots,m\right\}$.

Let us consider the case $\Omega=\left\{x\in\mathbb{R}^{n}:Cx\leq D\right\}$ where the number of inequality constraints is equal to $m$. We can write:

[TABLE]

where $\Omega_{j}=\left\{x\in\mathbb{R}^{n}:c_{\left(j\right)}^{\top}x\leq d_{\left(j\right)}\right\}$, $c_{\left(j\right)}^{\top}$ corresponds to the $j^{\mathrm{th}}$ row of $C$ and $d_{\left(j\right)}$ is the $j^{\mathrm{th}}$ element of $D$. We follow Tibshirani (2017) to define the corresponding algorithm. In particular, we introduce two iteration indices $j$ and $k$. The index $j$ refers to the constraint number, whereas the index $k$ refers to the main loop. Algorithm 5 describes the Dykstra’s approach for solving this proximal problem.

If we define $\Omega$ as follows:

[TABLE]

we decompose $\Omega$ as the intersection of three basic convex sets:

[TABLE]

where $\Omega_{1}=\left\{x\in\mathbb{R}^{n}:Ax=B\right\}$, $\Omega_{2}=\left\{x\in\mathbb{R}^{n}:Cx\leq D\right\}$ and $\Omega_{3}=\left\{x\in\mathbb{R}^{n}:x^{-}\leq x\leq x^{+}\right\}$. Using Dykstra’s algorithm is equivalent to formulating Algorithm 6.

Remark 11

*An alternative approach is to write the constraints in the following way191919We use the following properties:

$Ax=B\Leftrightarrow Ax\leq B\text{ and }Ax\geq B$

and:

$x^{-}\leq x\leq x^{+}\Leftrightarrow-x\leq-x^{-}\text{ and }x\leq x^{+}$

:*

[TABLE]

where $C^{\star}=\left(A,-A,C,-I_{n},I_{n}\right)$ and $D^{\star}=\left(B,-B,D,-x^{-},x^{+}\right)$. Therefore, we can use Algorithm 5 to find the solution.

A.4 Proximal operator of the risk budgeting logarithmic barrier

We have:

[TABLE]

where:

[TABLE]

The first-order condition is:

[TABLE]

We deduce that $z_{i}^{\left(k\right)}$ is the solution of the quadratic equation:

[TABLE]

We have:

[TABLE]

Since $\Delta>0$ and $-\lambda\varphi b_{i}<0$, we have two roots with opposite signs. Therefore, the solution is equal to:

[TABLE]

A.5 CCD algorithm with separable constraints

The coordinate update proposed by Nesterov (2012) and Wright (2015) is:

[TABLE]

where $\xi$ is a positive scalar, $\eta>0$ is the stepsize of the gradient descent and $g_{i}$ is the first-derivative of the function with respect to $x_{i}$:

[TABLE]

The objective function is equivalent to:

[TABLE]

By taking $\xi=\eta^{-1}$, we deduce that:

[TABLE]

where $\psi\left(x\right)=\mathds{1}_{\Omega_{i}}\left(x\right)$.

A.6 Accelerated bisection algorithm

The classic bisection algorithm consists in calculating the solution $x^{\star}\left(\lambda\right)=\arg\min\mathcal{L}\left(x;\lambda\right)$ and updating the bounds of the interval $\left[a_{\lambda},b_{\lambda}\right]$ that contains the solution $\lambda^{\star}$ such that $\sum\limits_{i=1}^{n}x^{\star}\left(\lambda\right)=1$. One of the issues is that $x^{\star}\left(\lambda\right)$ is obtained by an optimization algorithm and is not an analytical formula. This means that we can write:

[TABLE]

where $x^{\left(0\right)}$ is the starting value of the ADMM or CCD algorithm. Therefore, we can update the starting value $x^{\left(0\right)}$ of the algorithm at each iteration of the bisection method. This helps to accelerate the convergence of the $x$-update. In practice, we can replace Algorithm 1 by Algorithm 7.

A.7 Adaptive penalization parameter

The convergence of the ADMM algorithm holds regardless of the value of the penalization parameter $\varphi>0$. But the choice of $\varphi$ affects the speed of convergence. In practice, the penalization parameter $\varphi$ may be changed at each iteration, implying that $\varphi$ is replaced by $\varphi^{\left(k\right)}$. This may improve the convergence and make the performance of the ADMM algorithm less dependent of the initial choice $\varphi^{\left(0\right)}$. To update $\varphi^{\left(k\right)}$ in practice, He et al. (2000) and Wang and Liao (2001) provide a simple and efficient scheme. Let $r^{\left(k\right)}=Ax^{\left(k\right)}+Bz^{\left(k\right)}-c$ and $s^{\left(k\right)}=\varphi A^{\top}B\left(z^{\left(k\right)}-z^{\left(k-1\right)}\right)$ be the primal and dual residual variables (Boyd et al., 2011). On the one hand, the $x$ and $z$-updates essentially comes from placing a penalty on $\left\|r^{\left(k\right)}\right\|_{2}^{2}$. As a consequence, if $\varphi^{\left(k\right)}$ is large, $\left\|r^{\left(k\right)}\right\|_{2}^{2}$ tends to be small. On the other hand, $s^{\left(k\right)}$ depends linearly on $\varphi$. As a consequence, if $\varphi^{\left(k\right)}$ is small, $\left\|s^{\left(k\right)}\right\|_{2}^{2}$ is small and $\left\|r^{\left(k\right)}\right\|_{2}^{2}$ may be large. To keep $\left\|r^{\left(k\right)}\right\|_{2}^{2}$ and $\left\|s^{\left(k\right)}\right\|_{2}^{2}$ within a factor $\mu$, one may consider the following scheme:

If $\left\|r^{\left(k\right)}\right\|_{2}^{2}>\mu\left\|s^{\left(k\right)}\right\|_{2}^{2}$, the $\varphi$-update is:

[TABLE] 2. 2.

If $\left\|s^{\left(k\right)}\right\|_{2}^{2}>\mu\left\|r^{\left(k\right)}\right\|_{2}^{2}$, we have:

[TABLE] 3. 3.

Otherwise, the penalization parameter remains the same $\varphi^{\left(k+1\right)}\leftarrow\varphi^{\left(k\right)}$, implying that we do not rescale the dual variable $u^{\left(k+1\right)}$.

The previous scheme corresponds to the $\varphi$-update and must be placed after the $u$-update of the ADMM algorithm. In practice, we use the following default values: $\varphi^{\left(0\right)}=1$, $u^{\left(0\right)}=0$, $\mu=10^{6}$ and $\tau=\tau^{\prime}=2$.

Appendix B Implementation

A Python implementation is available on the following webpage:

https://github.com/jcrichard