Minimize quadratic function subject to linear equality constraints with SciPy

This is a partial answer to the question that I'm putting here to keep the question from getting even bigger, but I'd still love to see a more comprehensive and explanatory answer. These answers are based on comments from two others, but neither of them fully wrote out the code, and I thought it would make sense to make that explicit so here it is:

Fixing 2a (trust-constr with jacobian)

It's seems that the key here with regard to the Jacobian and Hessian is to specify neither or both (but not the jacobian only). @SubhaneilLahiri commented to this effect and there was also an error message to this effect that I initially failed to notice:

UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.

So I fixed it by defining the hessian function:

def constr_hess(x,v):
    return np.zeros([n,n])

and adding it to the constraint

nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac, constr_hess )

Fixing 3a & 3b (SLSQP)

This just seemed to be a matter of making the tolerance smaller as suggested by @user545424. So I just added options={'ftol':1e-15} to the minimization:

resSQjac = minimize( obj_func, x0, method='SLSQP',
                     options={'ftol':1e-15},
                     jac = obj_jac, constraints = eq_cons )

Here is how this problem could be solved using nlopt which is a library for nonlinear optimization which I've been pretty impressed with.

First, the objective function and gradient are both defined using the same function:

def obj_func(x, grad):
    if grad.size > 0:
        grad[:] = obj_jac(x)
    return ( ( ( x/x0 - 1 )) ** 2 ).sum()

def obj_jac(x):
    return 2. * ( x - x0 ) / x0 ** 2

def constr_func(x, grad):
    if grad.size > 0:
        grad[:] = constr_jac(x)
    return x.sum() - target

def constr_jac(x):
    return np.ones(n)

Then, to run the minimization using Nelder-Mead and SLSQP:

opt = nlopt.opt(nlopt.LN_NELDERMEAD,len(x0)-1)
opt.set_min_objective(unconstr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0[1:].copy())
xopt = np.hstack([target - xopt.sum(), xopt])
fval = opt.last_optimum_value()
print_res(xopt,fval,"Nelder-Mead");

opt = nlopt.opt(nlopt.LD_SLSQP,len(x0))
opt.set_min_objective(obj_func)
opt.add_equality_constraint(constr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0.copy())
fval = opt.last_optimum_value()
print_res(xopt,fval,"SLSQP w/ jacobian");

And here are the results:

 *****  Nelder-Mead  ***** 

obj func value at solution 0.00454545454546
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0


 *****  SLSQP w/ jacobian  ***** 

obj func value at solution 0.00454545454545
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0

When testing this out, I think I discovered what the issue with the original attempt was. If I set the absolute tolerance on the function to 1e-8 which is what the scipy functions default to I get:

 *****  Nelder-Mead  ***** 

obj func value at solution 0.0045454580693
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30816 41454 52274]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0


 *****  SLSQP w/ jacobian  ***** 

obj func value at solution 0.0146361108503
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10999 21000 31000 41000 51000]
% diff [9 5 3 2 2]
target achieved? 155000.0 155000.0

which is exactly what you were seeing. So my guess is that the minimizer ends up somewhere in the likelihood space during SLSQP where the next jump is less than 1e-8 from the last place.