How to understand this convex optimization method to find risk budget portfolio

Question

Both the short course material coded by the CVXPY developers and an answer on Quant SE suggest that given a desired risk budget $b$, we can find the full-investment portfolio with weights $w$ that has the risk budget (as defined in those materials) equal to $b$ by performing the following convex optimization:

$$\begin{align}\text{Minimize}\;&\frac{1}{2}w'\Sigma w - \sum_i b_i\,\text{log}w_i \\ \text{subject to} \;& 1'w=1\end{align}$$

However, having done this myself in Python with CVXPY, I found the resulting risk budgets were not the same as the desired $b$.

I then tried to calculate this minimization by hand, and I found that the solution $w$ has

$$(\Sigma w)_i - \frac{b_i}{w_i} = \lambda\;\;\;\text{for all }i$$

where $\lambda$ is the Lagrangian multiplier.

In other words, since we can show that the $i$-th risk budget ${b_w}_i$ by definition is equal to $\frac{w_i(\Sigma w)_i}{w'\Sigma w}$, the solution to that optimization problem has:

$${b_w}_i = \dfrac{b_i+\lambda w_i}{w'\Sigma w} = \dfrac{b_i+\lambda w_i}{1 + \lambda}$$

and this in general is not equal to $b_i$ (otherwise, $w_i=b_i$). This has been verified by the optimization that I ran in Python.

But I'm sure this method is not wrong--there's an academic paper written to explain it by Spinu (2013), which is beyond my capabilities. So, I'd really appreciate anyone who can explain this formulation!

Update:

Here's the Python code that I wrote. It's an exercise that's part of the CVXPY short course.

import numpy as np
import cvxpy as cp


#input data

Sigma = np.array([[6.1, 2.9, -0.8, 0.1], 
                  [2.9, 4.3, -0.3, 0.9], 
                  [-0.8, -0.3, 1.2, -0.7],
                  [0.1, 0.9, -0.7, 2.3]])
b = np.ones(4)/4  #risk parity


# optimization

w = cp.Variable(4)  #portfolio weight
obj = 0.5 * cp.quad_form(w, Sigma) - cp.sum(cp.multiply(b, cp.log(w)))  #objective
constr = [cp.sum(w) == 1, w >= 0] # constraint
prob = cp.Problem(cp.Minimize(obj), constr)
prob.solve()


# print the solution weight and solution risk budget

b_w = cp.multiply(w, Sigma @ w) / cp.quad_form(w, Sigma)  #solution risk budget
print("The solution weight is", w.value)
print("The solution risk budget is", b_w.value)

And the printed outputs are:

The solution weight is [0.16073365 0.14918463 0.42056612 0.2695156 ]

The solution risk budget is [0.32355772 0.33307394 0.10944985 0.23391849]

Answer 1

I believe the problem should be solved in two steps: the optimization without the condition that the sum equals 1, followed by a normalization step which divides w by sum(w) to produce the desired solution.

Here is the modified code:

import numpy as np
import cvxpy as cp


#input data

Sigma = np.array([[6.1, 2.9, -0.8, 0.1], 
                  [2.9, 4.3, -0.3, 0.9], 
                  [-0.8, -0.3, 1.2, -0.7],
                  [0.1, 0.9, -0.7, 2.3]])
b = np.ones(4)/4  #risk parity


# optimization

w = cp.Variable(4)  #portfolio weight
obj = 0.5 * cp.quad_form(w, Sigma) - cp.sum(cp.multiply(b, cp.log(w)))  
#objective
constr = [w >= 0] # constraint
prob = cp.Problem(cp.Minimize(obj), constr)
prob.solve()

# normalize

w = w/cp.sum(w)


# print the solution weight and solution risk budget

b_w = cp.multiply(w, Sigma @ w) / cp.quad_form(w, Sigma)  #solution risk budget
print("The solution weight is", w.value)
print("The solution risk budget is", b_w.value)

The result

The solution weight is [0.13765302 0.11336252 0.4758825  0.27310195]
The solution risk budget is [0.25000012 0.25000013 0.24999967 0.25000007]

The marginal values are equal (within numerical error).