Thursday, June 19, 2008

Median DNLP

The median (http://en.wikipedia.org/wiki/Median) can be found by minimizing sum(i, abs(x-a(i))). This is not a well-defined optimization model, and many solvers have problems with it.

An AMPL representation can be found here: http://orfe.princeton.edu/~rvdb/ampl/nlmodels/median/median.mod. A GAMS translation is:

$set m 19

set i /i1*i%m%/;

parameter a(i);
a(i) = uniform(0,1);

variable
x 'median'
z 'objective'
;

equation sumdist;

sumdist.. z =e= sum(i, abs(x-a(i)));

x.l = 0.5;

model median /sumdist/;
solve median minimizing z using dnlp;

Some results:
solverxzNote
Minos0.50024.9424The current point cannot be improved.
Snopt0.50024.9424The current point cannot be improved.
Conopt0.50024.9424Convergence too slow
Mosek1.864341E+123.542248E+13TERMINATED BY SOLVER
Ipopt0.50024.9424Restoration Failed!
Knitro0.5002-0.0002TERMINATED BY SOLVER
Baron,optcr=0,reslim=250.50024.9424Listing file has contradictory info
OQNLP0.50024.9424NORMAL COMPLETION


Baron mentions on the log that bounds are converging to 0.494241D+01, but in the listing file it says Best possible = -1E50.

Note: This model can be reformulated as an LP:

$set m 19

set i /i1*i%m%/;

parameter a(i);
a(i) = uniform(0,1);

variable
x 'median'
e(i) 'abs deviation'
z 'objective'
;

equation
sumdist
e1(i)
e2(i)
;

sumdist.. z =e= sum(i, e(i));
e1(i).. -e(i) =l= x-a(i);
e2(i).. x-a(i) =l= e(i);


model median /all/;
solve median minimizing z using lp;

3 comments:

  1. This is really interesting! Do you happen to know if this can be extended to other quantile calculations?

    ReplyDelete
  2. I don't know if anybody reads this after such a long time, but I think it should be noted, that the linear model doesn't work für an even number of Points. The model gives 8 instead of 7.5 for x with the following parameters:

    set i /i1*i6/;
    parameter a(i)
    /i1 1
    i2 2
    i3 7
    i4 8
    i5 1000
    i6 1000
    /;

    Apart from this, many thanks for your interesting blog.
    Best regards
    Joerg Schlutter

    ReplyDelete
  3. Depends on the definition. I would say median = 50% percentile is not uniquely defined for N even.

    ReplyDelete