sloisel / numeric

Numerical analysis in Javascript
http://www.numericjs.com/
Other
1.42k stars 176 forks source link

uncmin doesn't work #69

Open Meekohi opened 8 years ago

Meekohi commented 8 years ago

Seems like the most simple test possible... am I missing something?

numeric.uncmin(function(x){
  return x[0];
},[10])
f: 9.000000000037858
solution: 9.000000000037858
message: "Line search step size smaller than tol"
angeris commented 8 years ago

You're definitely right. Doing a vanilla pull of numeric and running the unit tests (under tools/build.sh) fails miserably in uncmin.

In particular, it seems that the inverse Hessian is incorrectly calculated for your constant f(x). Trying it out with my code on node.js (v8) as the engine I get:

> sq = function(x) {return x*x;}
[Function]
> numeric.uncmin(sq, [10])
{ solution: [ 2.3305801732931286e-11 ],
  f: 5.431603944147029e-22,
  gradient: [ 4.6611603465862566e-11 ],
  invHessian: [ [ 0.5000000000011653 ] ],
  iterations: 2,
  message: 'Newton step smaller than tol' }

which gets the right result. Trying it out on a simple linear function, screws up hilariously:

> cons = function(x) {return x;}
[Function]
> numeric.uncmin(cons, [10])
{ solution: [ 10 ],
  f: [ 10 ],
  gradient: [ 0 ],
  invHessian: [ [ 1 ] ],
  iterations: 0,
  message: 'Newton step smaller than tol' }

Note in particular that the hessian is zero, thus the inverse of the hessian is not defined. There should be an explicit check for this.

In terms of the unit test, it fails when you do test 179:

> 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;
}
[Function]
> f(numeric.uncmin(f,[1,2,1,1,1,1],1e-14).solution) < 1e-20
false

on the current head. I'm unsure as to why this is, but I feel it has to do with the computation of the current Hessian inverse.

Meekohi commented 8 years ago

In the two simple examples you gave I think you're supposed to reference x[0] instead of x inside those functions -- the function to minimized is passed an array, not a value. Maybe it is converting it by lucky chance to a number in that case.

angeris commented 8 years ago

You're right, my mistake. It evaluates to the same result on both cases, though, as

> cons = function(x) {return x[0];}
[Function]
> numeric.uncmin(cons, [10])
{ solution: [ 9.000000000037858 ],
  f: 9.000000000037858,
  gradient: [ 1.0000000000368852 ],
  invHessian: [ [ -13379220933.331156 ] ],
  iterations: 62,
  message: 'Line search step size smaller than tol' }
> sq = function(x) { return x[0]*x[0]; }
[Function]
> numeric.uncmin(sq, [10])
{ solution: [ 2.3305801732931286e-11 ],
  f: 5.431603944147029e-22,
  gradient: [ 4.6611603465862566e-11 ],
  invHessian: [ [ 0.5000000000011653 ] ],
  iterations: 2,
  message: 'Newton step smaller than tol' }

It still fails to correctly compute the inverse hessian (or check that it is completely ill-conditioned for invertibility). I fear this may be a result of the way the hessian is computed, but I'm unsure as it seems to be stuck on the line-search step for enough iterations that it exits there. In particular,

//[...] Some previous code defining things
if(typeof maxit === "undefined") maxit = 1000;
//[...] code that begins the newton-raphson method on the derivative

//Begin the line search
while(it < maxit) {
    if(t*nstep < tol) { break; }
    s = mul(step,t);
    x1 = add(x0,s);
    f1 = f(x1);
    if(f1-f0 >= 0.1*t*df0 || isNaN(f1)) {
        t *= 0.5;
        ++it;
        continue;
    }
    break;
}

This code begins in line 2782 of numeric.js from the repository.

carstenschwede commented 7 years ago

What is the current state of this issue?

angeris commented 7 years ago

Afaik, I don’t think there’s been any pull request accepted for a few years… unsure if there’s also been any active branches from this project.

On Mar 29, 2017, at 12:09 PM, Carsten Schwede notifications@github.com wrote:

What is the current state of this issue?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/sloisel/numeric/issues/69#issuecomment-290194353, or mute the thread https://github.com/notifications/unsubscribe-auth/ACZjr8rVdJvT3UXA_opOTxpgQh7vyEicks5rqqxkgaJpZM4GhSKy.