Whilst I agree with the general consensus of responders that this is not the best way to solve the minimisation problem in the question, I have now resolved the challenge and can answer my own question to share the way one might overcome similar issues in using penalty methods to resolve optimisation problems in Python.
The key mathematical issue is indeed the non-differentiability of the penalty functions; it seems that best practice is to use a polynomial of the same order as the objective function; in this way you ensure the behaviour of the penalty is complementary to your objective function. Considering <PERSON> series, it seems conceptually clear that you should be able to form a sum of these polynomials, of different coefficients, to well approximate almost any penalty you might desire.
The key programming issue (which will cause the SciPy error in the question) comes from providing incorrect Jacobian (and/or Hessian) functions to scipy.optimize.minimize(). The addition of the penalty function makes the calculation of the gradient vector and Hessian matrix considerably more difficult, and I had to calculate these by hand. 768
Fortunately SciPy provides a function to test your gradient function: check_grad(F, dF, x_k)), which compares the norm of your gradient function at x_k against an inbuilt finite difference approximation over a small region around x_k. Typically if this returns something $<10^{-4}$ then your function is likely correct (well, correct enough). This doesn't hold for the Hessian matrix, so more careful calculation is required.
So, for this question, I swapped the absolute, linear terms in the both penalty functions above with quadratic approximations:
def objective(x):
Q = DF.cov() # Covariance matrix
var = np.matmul(W.transpose(),np.matmul(Q,x)) # Variance vector
# Penalty Function method
penalty1 = (np.sum(x)-1)**2 # Large for sum(x) <> 1
penalty2 = 100*(R_min - np.dot(Mus, x))**2 # Large for returns <> minR
return var + penalty1 + penalty2
These are therefore twice differentiable and, after some painstaking matrix and vector derivative calculations, yielded gradient and Hessian functions as below:
def der_objective(x):
Q = DF.cov()
der = 2 * np.matmul(Q,x)
penalty1_der = np.array([2*np.sum(x)-2 for i,xi in enumerate(x)])
penalty2_der = 100*np.array([ 2*Mus[i]*(np.dot(Mus,x) - R_min) for i,xi in enumerate(x)])
return der + penalty1_der + penalty2_der
def hess_objective(x):
Q = DF.cov()
hess = 2 * Q.values
# Assemble Hessian terms for penalty1 function
penalty1_hess = 2*np.ones_like(hess)
# Assemble Hessian terms for penalty2 function
bcast = np.broadcast(Mus.reshape(len(Mus),1),Mus.reshape(1,len(Mus)))
penalty2_hess = np.empty(bcast.shape)
penalty2_hess.flat = 100*np.array([2*a*b for (a,b) in bcast])
return hess + penalty1_hess + penalty2_hess
These are then typically passed to the simple interface provided by SciPy to select a minimisation routine:
result = scipy.optimize.minimize(fun=objective, x0=x_k, method='BFGS',
jac=der_objective, hess=hess_objective,
options={'tol':1e-6,'maxiter':1e3}) 509