8
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Update: the link in Arnold's answer below is useful. For this particular problem though, it was helpful to rewrite it as a mixed integer problem (defining indicators as new variables) and use a mixed integer solver (I couldn't get CPLEX to work though, as the constraint was apparently not "quadratic enough" for it). $\endgroup$ – laxxy May 26 '12 at 15:06
6
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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)]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.