# Solving a system of two equations with non-convex matrix multiplication for MV optimization

Scenario: I am trying to do a variation of the MV optimization for a portfolio. In this instance, I already have a vector of mean returns ($\mu$), a vector of ones, a covariance matrix ($\Sigma$), and $\phi^{-1}$ which is the inverse of the standard normal cumulative distribution function (function of probability alpha).

Problem: From the two equations stated below, I am going to input: $\mu$, $\Sigma$, 1 (vector of ones), $H$, and $\alpha$, and try to get a value for $\gamma$ (also, the ' represents the transpose of the vector).

Equations:

$$H = w(\gamma)'\mu + \Phi^{−1}(α)[w(\gamma)' \Sigma w(γ)]^{1/2}$$ $$w(\gamma) = \frac{1}{\gamma}\Sigma^{-1}\Big[\mu − \Big(\frac{\mathbf{1}' \Sigma^{-1} \mu-\gamma}{\mathbf{1}' \Sigma^{-1} \mathbf{1}}\Big) \mathbf{1}\Big]$$

What I already tried: I am trying to solve this system using a python solver. But I am not able to do it since the matrix multiplication apparently becomes non-convex.

Question: How to solve this problem?

Obs: This question comes from the article Portfolio Optimization with mental accounts from Sanjiv Das, Harry Markowitz, Jonathan Scheid, and Meir Statman, Equations (7) and (8). The article is found at: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.410.8747&rep=rep1&type=pdf

• Starting with a value of $\gamma$, you can compute $w(\gamma)$ from the second equation and plug it into the RHS of the first equation and see what you get. If it does not match the desired $H$ you can try again with a different $\gamma$. This is what the article seems to suggest. Are you saying that this process does not work/fails to converge? Sep 20, 2017 at 15:31
• @noob2 Thanks for the comment. As the article states, specifying gamma with precision is not intuitive, so what I am trying to do is set up a optimization procedure that allows me to solve for it (given H and alpha). What you said is correct, I can just try values of gamma until I reach the desired result, but I don't know to to optimize that with a convex only solver. Sep 21, 2017 at 13:58

I will use a case in which I randomly generate both $\Sigma$ and $\mu$ randomly in 4 dimensions

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import bisect

n = 4
np.random.seed(0)

S = np.random.uniform(0, 1, [n, n])
S = 0.5 * (S + S.transpose())
mu = np.random.uniform(-1, 1, n)


The trick with the transpose is just there to ensure that S is symmetric. I will choose $\alpha = 0.90$, and define the function $H$ as

Sinv = np.linalg.inv(S)
ones = np.ones_like(mu)

def H(gamma):

# w
w = (ones.dot(Sinv).dot(mu) - gamma) / (ones.dot(Sinv).dot(ones))
w = Sinv.dot(mu - w) / gamma

# H
H = w.dot(mu) + norm.ppf(alpha) * np.sqrt(w.dot(S).dot(w))

return H


This is a plot of $H = H(\gamma)$.

 g = np.linspace(4, 50, num = 500)
h = np.array(map(H, g))
i = h > 0
plt.plot(g[i], h[i])
plt.show() The function is indeed sort of problematic for large values of $\gamma$ in the sense that any optimization algorithm (pretty much all of them) that uses derivatives will fail. But still you can use bisection to find the root

As an example I will use $H = 1.04$ (up in the flat region)

f = lambda x: H(np.exp(x)) - 1.04
r = bisect(f, 1.5, 6, xtol = 1e-12)
g0 = np.exp(r)
print g0, H(g0)


which generates the output

231.206991393 1.04


So the root is $\gamma = 231.20699$

• Many thanks for the answer, this is exactly what I wanted. Just a couple of questions: how to you plot the graph you used in the answer? And for that graph, would it be possible to plot it only for H values above zero (since that is the return I want from the portfolio)? Sep 27, 2017 at 15:49
• @DGMS89 I modified the answer to take include the commands for plotting $H>0$ Sep 27, 2017 at 16:02