How to Solve Optimisation Problems using Penalty Functions in Python

Question

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?

Answer 1

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