Minimizing a quadratic function with nonlinear constraints

Question

what would be good methods (and/or software packages) to try for solving a problem minimizing a quadratic function $f(x) = \sum_{i=1}^N{(x_i - y_i)^2}$, s.t. $0 \leq x_i \leq 1$, and there are more constraints some of which are non-linear (and non-differentiable), e.g. $\sum_i x_i \mathbf{1}_{x_i>a} < b$ ?

I am thinking about $N \approx 100$. FWIW, Matlab is apparently using an "active set method, similar to that of Gill et al.", which has somewhat uneven performance.

Answer 1

If you have nonsmooth constraints, it is no help that you have a convex quadratic objective. You need a constrained nonsmooth solver.

See my web page
http://www.mat.univie.ac.at/~neum/glopt/software_l.html#nonsm
for suitable software.

Answer 2

You say in the comment that you cannot get it to work as it is not quadratic enough. I don't see any reason for that. The problem is easily coded as a mixed-integer quadratic program.

If I understand your problem definition, you want to constraint the sum of the varables larger than a threshold. Introduce a binary variable indicating if x is larger than a, and introduce another variable z which should be equal to x when this holds, and zero otherwise, and use the sum of the new variables.

Using the MATLAB toolbox YALMIP to interface CPLEX (or Gurobi or any other MIQP solver), the problem is trivially solved in fractions of a second. Here, a random example, implemented using both a manually derived model, and a model which exploits the high-level modelling capabilities in YALMIP

% Create random data
N = 100;
y = rand(N,1);
a = rand(N,1);
b = rand(1);

% Decision variables
x = sdpvar(N,1);
z = sdpvar(N,1);
d = binvar(N,1);

% Define objective
Objective = (x-y)'*(x-y)

% High-level model
Con1 = [0 <= x <= 1,0 <= z <=1];
Con2 = [implies(x-a>=0,d), implies(d,z==x), implies(1-d,z==0)]
Con3 = [sum(z) <= b];

% Solve problem
solvesdp([Con1,Con2,Con3],Objective)

% display solution
[double(x) a double(z)]

% Manually derived model
Con1 = [0 <= x <= 1,0 <= z <=1];
Con2 = [x-a <= d, -(1-d) <= x-z <= 1-d, z <= d];
Con3 = [sum(z) <= b];
Objective = (x-y)'*(x-y)
solvesdp([Con1,Con2,Con3],Objective)
[double(x) a double(z)]