Software 10.5.1 (Solving the KKT conditions). Example 10.2.2: for each of the \(2^3\) choices of active constraints, solve stationarity together with \(g_i = 0\) (active) or \(\lambda _i = 0\) (inactive), and keep the solutions that are feasible with \(\lambda \geq 0\):

>>> import sympy as sp
>>> from itertools import product
>>> x, y = sp.symbols('x y', real=True)
>>> f = (x - 2)**2 + (y - 1)**2            # minimise f subject to g_i <= 0
>>> g = [x + y - 2, -x, -y]
>>> lam = sp.symbols('l1:4', real=True)
>>> L = f + sum(l * gi for l, gi in zip(lam, g))
>>> for active in product([0, 1], repeat=3):        # guess the active set
...     eqs = [L.diff(x), L.diff(y)]
...     eqs += [gi if a else l for a, gi, l in zip(active, g, lam)]
...     for s in sp.solve(eqs, [x, y, *lam], dict=True):
...         ok = (all(gi.subs(s) <= 0 for gi in g)
...               and all(s[l] >= 0 for l in lam))
...         if ok: print(active, s, f.subs(s))
...
(1, 0, 0) {l1: 1, l2: 0, l3: 0, x: 3/2, y: 1/2} 1/2

Reading the output. Only the active set \(\{x + y \leq 2\}\) gives a KKT point: \((\frac 32, \frac 12)\) with \(\lambda _1 = 1\) and value \(\frac 12\). Enumerating active sets is exponential in general, but it shows exactly what the conditions say.

MATLAB/Octave: solve with the Symbolic Math Toolbox on the same systems.

Software 10.5.2 (The efficient frontier). The portfolio QP solved by SLSQP for five required returns, giving the weights and the standard deviation; then the minimum-variance portfolio from the formula \(w = \Sigma ^{-1}\mathbf 1/\mathbf 1^{\mathsf T}\Sigma ^{-1}\mathbf 1\), which ignores \(w \geq 0\):

>>> import numpy as np
>>> from scipy.optimize import minimize
>>> mu = np.array([0.08, 0.12, 0.15])      # expected returns of three assets
>>> S = np.array([[0.04, 0.006, 0.010],    # covariance matrix of the returns
...               [0.006, 0.09, 0.030],
...               [0.010, 0.030, 0.16]])
>>> def frontier(target):  # min w S w, mu w >= target, sum w = 1
...     cons = [{'type': 'eq', 'fun': lambda w: w.sum() - 1},
...             {'type': 'ineq', 'fun': lambda w: mu @ w - target}]
...     r = minimize(lambda w: w @ S @ w, np.ones(3) / 3,
...                  jac=lambda w: 2 * S @ w, constraints=cons,
...                  bounds=[(0, 1)] * 3, method='SLSQP')
...     return r.x, np.sqrt(r.x @ S @ r.x)
...
>>> for t in (0.08, 0.10, 0.12, 0.14, 0.15):
...     w, sd = frontier(t)
...     print(t, np.round(w, 3), round(sd, 4))
...
0.08 [0.662 0.245 0.093] 0.1699
0.1 [0.6   0.268 0.133] 0.1712
0.12 [0.262 0.389 0.349] 0.2169
0.14 [0.    0.333 0.667] 0.3073
0.15 [0. 0. 1.] 0.4
>>> one = np.ones(3); Si = np.linalg.inv(S)  # unconstrained minimum variance
>>> wmv = Si @ one / (one @ Si @ one); np.round(wmv, 4), round(mu @ wmv, 4)
(array([0.6614, 0.2456, 0.093 ]), 0.0963)

Reading the output. For required returns up to \(0.096\) the constraint is not active and the answer is the minimum-variance portfolio \((0.661, 0.246, 0.093)\) with \(\sigma = 0.170\); the formula gives the same weights, all positive, so the bounds were not needed. For higher returns the weight shifts to assets \(2\) and \(3\): at \(0.14\), asset \(1\) is dropped entirely, and a return of \(0.15\) requires holding only asset \(3\), at \(\sigma = 0.4\).

MATLAB/Octave: quadprog(2*S, [], -mu’, -r, ones(1,3), 1, zeros(3,1)).

Software 10.5.3 (A barrier method). Minimise \(x_1^2 + x_2^2\) subject to \(x_1x_2 \geq 1\), by minimising \(x_1^2 + x_2^2 - t\log (x_1x_2 - 1)\) for decreasing \(t\); each line shows \(t\), the minimiser, the objective and the multiplier estimate \(t/(x_1x_2 - 1)\); then SLSQP:

>>> def barrier(t):  # min x1^2 + x2^2 - t log(x1 x2 - 1)
...     def phi(z):
...         gap = z[0] * z[1] - 1
...         return z @ z - t * np.log(gap) if gap > 0 else np.inf
...     opts = {'xatol': 1e-10, 'fatol': 1e-12, 'maxiter': 4000}
...     return minimize(phi, [2.0, 2.0], method='Nelder-Mead', options=opts).x
...
>>> for t in (1.0, 0.1, 0.01, 0.001):
...     z = barrier(t); lam_est = t / (z[0] * z[1] - 1)  # multiplier
...     print(t, np.round(z, 5), round(z @ z, 5), round(lam_est, 4))
...
1.0 [1.22474 1.22474] 3.0 2.0
0.1 [1.0247 1.0247] 2.1 2.0
0.01 [1.0025 1.0025] 2.01 2.0
0.001 [1.00025 1.00025] 2.001 2.0
>>> cons = [{'type': 'ineq', 'fun': lambda z: z[0] * z[1] - 1}]
>>> r = minimize(lambda z: z @ z, [2.0, 0.6], constraints=cons,
...              method='SLSQP')
>>> np.round(r.x, 6), round(r.fun, 6)
(array([0.999999, 1.000001]), 2.0)

Reading the output. The minimisers \((1.2247, 1.2247)\), \((1.0247, 1.0247)\), …approach the optimum \((1, 1)\) from inside the feasible region, and the objective \(2 + t\) exceeds the optimum \(2\) by exactly \(mt\) with \(m = 1\). The multiplier estimate is \(2\) throughout, the exact KKT multiplier: at \((1, 1)\), \(\nabla f = (2, 2) = 2\nabla (x_1x_2)\). SLSQP finds \((1, 1)\) directly.

MATLAB/Octave: fmincon(@(z) z*z’, [2 0.6], [], [], [], [], [], [], @(z) deal(1 - z(1)*z(2), [])).