# How to Solve Optimisation Problems using Penalty Functions in Python

I am working on a implementing a simple quadratic optimisation problem:

$$\min _x \; {\underline{x}}^T Q {\underline{x}}$$ $$s.t. \,\quad {\underline{\mu}}^T{\underline{x}} = R^*$$ $$\quad \quad \underline{1}^T\underline{x} = 1$$ I also expect to go on to include the inequality constraint, as an additional complexity, once the above is working.

$$x_i \geq 0$$ The method I think is simplest, and which I understand best for implementing these constraints, is the penalty function method, where we modify the objective function to 'steer' the optimisation away from forbidden regions. By carefully parameterising the size of the penalties, I have achieved good results using SciPy's built-in Nelder-Mead Simplex algorithm, using the objective function below.

def objective(x):
Q = DF.cov()     # Covariance matrix

# Penalty Function method

penalty1 =  0.0005 * abs(np.sum(x)-1)                             # Large for sum(x) <> 1
penalty2 =  0.05 * abs(R_min - np.matmul(Mus.transpose(), x))     # Large for returns <> R_min

return np.matmul(x.transpose(),np.matmul(Q,x)) + penalty1 + penalty2


Now, I wish to use other optimisation algorithms (in particular BFGS and Newton-CG), which require the gradient and Hessian of the objective function. I have implemented the derivative functions in the unconstrained case, but by adding the penalty terms to the objective (and the derivatives of the penalties to the gradient function) optimisation fails with the following error:

Warning: Desired error not necessarily achieved due to precision loss.
Current function value: 0.000056
Iterations: 0
Function evaluations: 780
Gradient evaluations: 96


(Previously Iterations would be a few hundred). This strictly occurs with penalty1, but not penalty2 by itself, so either my derivative is wrong for penalty1:

penalty1_der = np.sign(x)


Or can I not use the L1 norm in this way? I also tried replacing the constraint with a smoother, quadratic approximation:

penalty1 = np.matmul((x - vector_ones).transpose(), (x - vector_ones))


but unfortunately, although this prevents the error, Minimize() seems to completely ignore my penalty functions (even with vastly increased parameters).

How can I implement my constraints such that I can solve the problem using BFGS/Newton-CG?

• The L1 penalty isn’t differential enough and thus this method isn’t applicable. – Brian Borchers Mar 31 '18 at 0:33
• Thanks, Brian - are you saying it would need to be twice differentiable? Is that true also for BFGS (since this doesn't require the Hessian)? Finally, I'm assuming a sigmoid or other smooth approximation could still work instead? – Zac Mar 31 '18 at 0:39
• Could you write out the actual constraints that you're trying to impose? It's likely that we can help to suggest either a more effective penalization or another way to solve the problem. It should be noted that if you have only equality constraints like $\sum_i x_i = 1$, the optimization problem has a closed-form solution, and you need not go through the hassle of an iterative procedure. – Tyler Olsen Mar 31 '18 at 0:39
• Please write out your problem in mathematical notation and explain what you are trying to accomplish. – Brian Borchers Mar 31 '18 at 0:49
• Just use a Quadratic Programming (QP) solver, and don't use penalty terms. – Mark L. Stone Mar 31 '18 at 1:29

## 1 Answer

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 Taylor 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. 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})