# sloisel/numeric

New feature: linear programming (numeric.solveLP)

1 parent 9eb4e48 commit fd3831f8bd45fc3bde6c2affe4e96909e1350935 Sébastien Loisel committed Sep 1, 2012
File renamed without changes.
 @@ -137,6 +137,7 @@
setBlockSet a block of a matrix
sinPointwise Math.sin(x)
solveSolve Ax=b +
solveLPSolve a linear programming problem
splineCreate a Spline object
Spline.atEvaluate the Spline at a point @@ -398,6 +399,17 @@

Utility functions

OUT> [1,1.5,2,2.5,3] + +

Arithmetic operations

The standard arithmetic operations have been vectorized: @@ -1051,31 +1063,84 @@

Unconstrained optimization

We also include a simple unconstrained optimization routine. Here are some demos from from Moré et al., 1981:
-IN> sqr = function(x) { return x*x; }; numeric.uncmin(function(x) { return sqr(10*(x[1]-x[0]*x[0])) + sqr(1-x[0]); },[-1.2,1]).solution
+IN> sqr = function(x) { return x*x; };
+    numeric.uncmin(function(x) { return sqr(10*(x[1]-x[0]*x[0])) + sqr(1-x[0]); },[-1.2,1]).solution
OUT> [1,1]
-IN> f = function(x) { return sqr(-13+x[0]+((5-x[1])*x[1]-2)*x[1])+sqr(-29+x[0]+((x[1]+1)*x[1]-14)*x[1]); }; x0 = numeric.uncmin(f,[0.5,-2]).solution
+IN> f = function(x) { return sqr(-13+x[0]+((5-x[1])*x[1]-2)*x[1])+sqr(-29+x[0]+((x[1]+1)*x[1]-14)*x[1]); };
+    x0 = numeric.uncmin(f,[0.5,-2]).solution
OUT> [11.41,-0.8968]
-IN> f = function(x) { return sqr(1e4*x[0]*x[1]-1)+sqr(Math.exp(-x[0])+Math.exp(-x[1])-1.0001); }; x0 = numeric.uncmin(f,[0,1]).solution
+IN> f = function(x) { return sqr(1e4*x[0]*x[1]-1)+sqr(Math.exp(-x[0])+Math.exp(-x[1])-1.0001); };
+    x0 = numeric.uncmin(f,[0,1]).solution
OUT> [1.098e-5,9.106]
-IN> f = function(x) { return sqr(x[0]-1e6)+sqr(x[1]-2e-6)+sqr(x[0]*x[1]-2)}; x0 = numeric.uncmin(f,[0,1]).solution
+IN> f = function(x) { return sqr(x[0]-1e6)+sqr(x[1]-2e-6)+sqr(x[0]*x[1]-2)};
+    x0 = numeric.uncmin(f,[0,1]).solution
OUT> [1e6,2e-6]
-IN> f = function(x) { return sqr(1.5-x[0]*(1-x[1]))+sqr(2.25-x[0]*(1-x[1]*x[1]))+sqr(2.625-x[0]*(1-x[1]*x[1]*x[1])); }; x0 = numeric.uncmin(f,[1,1]).solution
+IN> f = function(x) { return sqr(1.5-x[0]*(1-x[1]))+sqr(2.25-x[0]*(1-x[1]*x[1]))+sqr(2.625-x[0]*(1-x[1]*x[1]*x[1])); };
+    x0 = numeric.uncmin(f,[1,1]).solution
OUT> [3,0.5]
-IN> f = function(x) { var ret = 0,i; for(i=1;i<=10;i++) ret+=sqr(2+2*i-Math.exp(i*x[0])-Math.exp(i*x[1])); return ret; }; x0 = numeric.uncmin(f,[0.3,0.4]).solution
+IN> f = function(x) {
+        var ret = 0,i;
+        for(i=1;i<=10;i++) ret+=sqr(2+2*i-Math.exp(i*x[0])-Math.exp(i*x[1]));
+         return ret;
+    };
+    x0 = numeric.uncmin(f,[0.3,0.4]).solution
OUT> [0.2578,0.2578]
-IN> y = [0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,0.37,0.58,0.73,0.96,1.34,2.10,4.39]; f = function(x) { var ret = 0,i; for(i=1;i<=15;i++) ret+=sqr(y[i-1]-(x[0]+i/((16-i)*x[1]+Math.min(i,16-i)*x[2]))); return ret; }; x0 = numeric.uncmin(f,[1,1,1]).solution
+IN> y = [0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,0.37,0.58,0.73,0.96,1.34,2.10,4.39];
+    f = function(x) {
+        var ret = 0,i;
+        for(i=1;i<=15;i++) ret+=sqr(y[i-1]-(x[0]+i/((16-i)*x[1]+Math.min(i,16-i)*x[2])));
+        return ret;
+    };
+    x0 = numeric.uncmin(f,[1,1,1]).solution
OUT> [0.08241,1.133,2.344]
-IN> y = [0.0009,0.0044,0.0175,0.0540,0.1295,0.2420,0.3521,0.3989,0.3521,0.2420,0.1295,0.0540,0.0175,0.0044,0.0009]; f = function(x) { var ret = 0,i; for(i=1;i<=15;i++) ret+=sqr(x[0]*Math.exp(-x[1]*sqr((8-i)/2-x[2])/2)-y[i-1]); return ret; }; x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[1,1,1]).solution,1000)),1000)
+IN> y = [0.0009,0.0044,0.0175,0.0540,0.1295,0.2420,0.3521,0.3989,0.3521,0.2420,0.1295,0.0540,0.0175,0.0044,0.0009];
+    f = function(x) {
+        var ret = 0,i;
+        for(i=1;i<=15;i++)
+        ret+=sqr(x[0]*Math.exp(-x[1]*sqr((8-i)/2-x[2])/2)-y[i-1]);
+        return ret;
+    };
+    x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[1,1,1]).solution,1000)),1000)
OUT> [0.399,1,0]
-IN> f = function(x) { return sqr(x[0]+10*x[1])+5*sqr(x[2]-x[3])+sqr(sqr(x[1]-2*x[2]))+10*sqr(x[0]-x[3]); }; x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[3,-1,0,1]).solution,1e5)),1e5)
+IN> f = function(x) { return sqr(x[0]+10*x[1])+5*sqr(x[2]-x[3])+sqr(sqr(x[1]-2*x[2]))+10*sqr(x[0]-x[3]); };
+    x0 = numeric.div(numeric.round(numeric.mul(numeric.uncmin(f,[3,-1,0,1]).solution,1e5)),1e5)
OUT> [0,0,0,0]
-IN> f = function(x) { return sqr(10*(x[1]-x[0]*x[0]))+sqr(1-x[0])+90*sqr(x[3]-x[2]*x[2])+sqr(1-x[2])+10*sqr(x[1]+x[3]-2)+0.1*sqr(x[1]-x[3]); }; x0 = numeric.uncmin(f,[-3,-1,-3,-1]).solution
+IN> f = function(x) {
+        return (sqr(10*(x[1]-x[0]*x[0]))+sqr(1-x[0])+
+                90*sqr(x[3]-x[2]*x[2])+sqr(1-x[2])+
+                10*sqr(x[1]+x[3]-2)+0.1*sqr(x[1]-x[3])); };
+    x0 = numeric.uncmin(f,[-3,-1,-3,-1]).solution
OUT> [1,1,1,1]
-IN> y = [0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246]; u = [4,2,1,0.5,0.25,0.167,0.125,0.1,0.0833,0.0714,0.0625]; f = function(x) { var ret=0, i; for(i=0;i<11;++i) ret += sqr(y[i]-x[0]*(u[i]*u[i]+u[i]*x[1])/(u[i]*u[i]+u[i]*x[2]+x[3])); return ret; }; x0 = numeric.uncmin(f,[0.25,0.39,0.415,0.39]).solution
+IN> y = [0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246];
+    u = [4,2,1,0.5,0.25,0.167,0.125,0.1,0.0833,0.0714,0.0625];
+    f = function(x) {
+        var ret=0, i;
+        for(i=0;i<11;++i) ret += sqr(y[i]-x[0]*(u[i]*u[i]+u[i]*x[1])/(u[i]*u[i]+u[i]*x[2]+x[3]));
+        return ret;
+    };
+    x0 = numeric.uncmin(f,[0.25,0.39,0.415,0.39]).solution
OUT> [     0.1928,     0.1913,     0.1231,     0.1361]
-IN> y = [0.844,0.908,0.932,0.936,0.925,0.908,0.881,0.850,0.818,0.784,0.751,0.718,0.685,0.658,0.628,0.603,0.580,0.558,0.538,0.522,0.506,0.490,0.478,0.467,0.457,0.448,0.438,0.431,0.424,0.420,0.414,0.411,0.406]; f = function(x) { var ret=0, i; for(i=0;i<33;++i) ret += sqr(y[i]-(x[0]+x[1]*Math.exp(-10*i*x[3])+x[2]*Math.exp(-10*i*x[4]))); return ret; }; x0 = numeric.uncmin(f,[0.5,1.5,-1,0.01,0.02]).solution
+IN> y = [0.844,0.908,0.932,0.936,0.925,0.908,0.881,0.850,0.818,0.784,0.751,0.718,
+         0.685,0.658,0.628,0.603,0.580,0.558,0.538,0.522,0.506,0.490,0.478,0.467,
+         0.457,0.448,0.438,0.431,0.424,0.420,0.414,0.411,0.406];
+    f = function(x) {
+        var ret=0, i;
+        for(i=0;i<33;++i) ret += sqr(y[i]-(x[0]+x[1]*Math.exp(-10*i*x[3])+x[2]*Math.exp(-10*i*x[4])));
+        return ret;
+    };
+    x0 = numeric.uncmin(f,[0.5,1.5,-1,0.01,0.02]).solution
OUT> [     0.3754,      1.936,     -1.465,    0.01287,    0.02212]
-IN> f = function(x) { var ret=0, i,ti,yi,exp=Math.exp; for(i=1;i<=13;++i) { ti = 0.1*i; yi = exp(-ti)-5*exp(-10*ti)+3*exp(-4*ti); ret += sqr(x[2]*exp(-ti*x[0])-x[3]*exp(-ti*x[1])+x[5]*exp(-ti*x[4])-yi); } return ret; }; x0 = numeric.uncmin(f,[1,2,1,1,1,1],1e-14).solution; f(x0)<1e-20;
+IN> f = function(x) {
+        var ret=0, i,ti,yi,exp=Math.exp;
+        for(i=1;i<=13;++i) {
+            ti = 0.1*i;
+            yi = exp(-ti)-5*exp(-10*ti)+3*exp(-4*ti);
+            ret += sqr(x[2]*exp(-ti*x[0])-x[3]*exp(-ti*x[1])+x[5]*exp(-ti*x[4])-yi);
+        }
+        return ret;
+    };
+    x0 = numeric.uncmin(f,[1,2,1,1,1,1],1e-14).solution;
+    f(x0)<1e-20;
OUT> true
There are optional parameters to numeric.uncmin(f,x0,tol,gradient,maxit,callback). The iteration stops when @@ -1099,6 +1164,70 @@

Unconstrained optimization

{i:6, x:[1.571], f:-1 , g: [ 2.242e-11] , H:[[0.25 ]] }] +

Linear programming

+ +Linear programming is available: + +
+IN> x = numeric.solveLP([1,1],                   /* minimize [1,1]*x                */
+                        [[-1,0],[0,-1],[-1,-2]], /* matrix of inequalities          */
+                        [0,0,-3]                 /* right-hand-side of inequalities */
+                        );
+    numeric.trunc(x.solution,1e-12);
+OUT> [0,1.5]
+
+ +The function numeric.solveLP(c,A,b) minimizes dot(c,x) subject to dot(A,x) <= b. +The algorithm used is very basic. For alpha>0, define the ``barrier function'' +f0 = dot(c,x) - alpha*sum(log(b-A*x)). The function numeric.solveLP calls numeric.uncmin +on f0 for smaller and smaller values of alpha. This is a basic ``interior point method''. + +We also handle infeasible problems: +
+IN> numeric.solveLP([1,1],[[1,0],[0,1],[-1,-1]],[-1,-1,-1])
+OUT> { solution: NaN, message: "Infeasible" }
+
+ +Unbounded problems: +
+IN> numeric.solveLP([1,1],[[1,0],[0,1]],[0,0]);
+OUT> { solution:[-1.009,-1.009], message:"Unbounded" }
+
+ +With an equality constraint: +
+IN> numeric.solveLP([1,2,3],                      /* minimize [1,2,3]*x                */
+                    [[-1,0,0],[0,-1,0],[0,0,-1]], /* matrix A of inequality constraint */
+                    [0,0,0],                      /* RHS b of inequality constraint    */
+                    [[1,1,1]],                    /* matrix Aeq of equality constraint */
+                    [3]                           /* vector beq of equality constraint */
+                    );
+OUT> { solution: [3,2.355e-17,9.699e-17], message:"" }
+
+ + +

Solving ODEs

The function numeric.dopri() is an implementation of the Dormand-Prince-Runge-Kutta integrator with
 @@ -1 +1,2 @@ *.js +*.tar.gz