fix: integrate fmin internally

__tests__/unit/data/venn/fmin/fmin.spec.ts

@@ -0,0 +1,132 @@
+import {
+  conjugateGradient,
+  conjugateGradientSolve,
+  gradientDescent,
+  gradientDescentLineSearch,
+  nelderMead,
+} from '../../../../../src/data/utils/venn/fmin';
+
+const SMALL = 1e-5;
+
+function nearlyEqual(
+  left,
+  right,
+  tolerance = SMALL,
+  message = 'assertNearlyEqual',
+) {
+  expect(Math.abs(left - right)).toBeLessThan(tolerance);
+  console.log(`${message}: ${left} ~== ${right}`);
+}
+
+function lessThan(test, left, right, message) {
+  message = message || 'lessThan';
+  test.ok(left < right, `${message}: ${left} < ${right}`);
+}
+
+const optimizers = [
+  nelderMead,
+  gradientDescent,
+  gradientDescentLineSearch,
+  conjugateGradient,
+];
+
+const optimizerNames = [
+  'Nelder Mead',
+  'Gradient Descent',
+  'Gradient Descent w/ Line Search',
+  'Conjugate Gradient',
+];
+
+describe('fmin', () => {
+  test('himmelblau', () => {
+    // due to a bug, this used to not converge to the minimum
+    const x = 4.9515014216303825;
+    const y = 0.07301421370357275;
+
+    const params = { learnRate: 0.1 };
+
+    const himmelblau = (X, fxprime = [0, 0]) => {
+      const [x, y] = X;
+      fxprime[0] = 2 * (x + 2 * y - 7) + 4 * (2 * x + y - 5);
+      fxprime[1] = 4 * (x + 2 * y - 7) + 2 * (2 * x + y - 5);
+      // biome-ignore lint/style/useExponentiationOperator: TODO: use **
+      return Math.pow(x + 2 * y - 7, 2) + Math.pow(2 * x + y - 5, 2);
+    };
+
+    optimizers.forEach((optimizer, index) => {
+      const solution = optimizer(himmelblau, [x, y], params);
+      nearlyEqual(solution.fx, 0, SMALL, `himmelblau:${optimizerNames[index]}`);
+    });
+  });
+
+  test('banana', () => {
+    const x = 1.6084564160555601;
+    const y = -1.5980748860165477;
+
+    const banana = (X, fxprime) => {
+      fxprime = fxprime || [0, 0];
+      const x = X[0];
+      const y = X[1];
+      fxprime[0] = 400 * x * x * x - 400 * y * x + 2 * x - 2;
+      fxprime[1] = 200 * y - 200 * x * x;
+      return (1 - x) * (1 - x) + 100 * (y - x * x) * (y - x * x);
+    };
+
+    const params = { learnRate: 0.0003, maxIterations: 50000 };
+    for (let i = 0; i < optimizers.length; ++i) {
+      const solution = optimizers[i](banana, [x, y], params);
+      nearlyEqual(solution.fx, 0, 1e-3, `banana:${optimizerNames[i]}`);
+    }
+  });
+
+  test('quadratic1D', () => {
+    const loss = (x, xprime) => {
+      xprime = xprime || [0, 0];
+      xprime[0] = 2 * (x[0] - 10);
+      return (x[0] - 10) * (x[0] - 10);
+    };
+
+    const params = { learnRate: 0.5 };
+
+    for (let i = 0; i < optimizers.length; ++i) {
+      const solution = optimizers[i](loss, [0], params);
+      nearlyEqual(solution.fx, 0, SMALL, `quadratic_1d:${optimizerNames[i]}`);
+    }
+  });
+
+  test('nelderMead', () => {
+    const loss = (X) => {
+      const x = X[0];
+      const y = X[1];
+      return Math.sin(y) * x + Math.sin(x) * y + x * x + y * y;
+    };
+
+    const solution = nelderMead(loss, [-3.5, 3.5]);
+    nearlyEqual(solution.fx, 0, SMALL, 'nelderMead');
+  });
+
+  test('conjugateGradientSolve', () => {
+    // matyas function
+    let A = [
+      [0.52, -0.48],
+      [-0.48, 0.52],
+    ];
+    let b = [0, 0];
+    const initial = [-9.08, -7.83];
+    let x = conjugateGradientSolve(A, b, initial);
+    nearlyEqual(x[0], 0, SMALL, 'matyas.x');
+    nearlyEqual(x[1], 0, SMALL, 'matyas.y');
+
+    // booth's function
+    const history = [];
+    A = [
+      [10, 8],
+      [8, 10],
+    ];
+    b = [34, 38];
+    x = conjugateGradientSolve(A, b, initial, history);
+    nearlyEqual(x[0], 1, SMALL, 'booth.x');
+    nearlyEqual(x[1], 3, SMALL, 'booth.y');
+    console.log(history);
+  });
+});

package.json

@@ -78,7 +78,6 @@
     "@antv/util": "^3.3.10",
     "@antv/vendor": "^1.0.8",
     "flru": "^1.0.2",
-    "fmin": "0.0.2",
     "pdfast": "^0.2.0"
   },
   "devDependencies": {

src/data/utils/venn/fmin/bisect.ts

@@ -0,0 +1,32 @@
+/** finds the zeros of a function, given two starting points (which must
+ * have opposite signs */
+export function bisect(f, a, b, parameters?: any) {
+  parameters = parameters || {};
+  const maxIterations = parameters.maxIterations || 100;
+  const tolerance = parameters.tolerance || 1e-10;
+  const fA = f(a);
+  const fB = f(b);
+  let delta = b - a;
+
+  if (fA * fB > 0) {
+    throw 'Initial bisect points must have opposite signs';
+  }
+
+  if (fA === 0) return a;
+  if (fB === 0) return b;
+
+  for (let i = 0; i < maxIterations; ++i) {
+    delta /= 2;
+    const mid = a + delta;
+    const fMid = f(mid);
+
+    if (fMid * fA >= 0) {
+      a = mid;
+    }
+
+    if (Math.abs(delta) < tolerance || fMid === 0) {
+      return mid;
+    }
+  }
+  return a + delta;
+}

src/data/utils/venn/fmin/blas1.ts

@@ -0,0 +1,42 @@
+// need some basic operations on vectors, rather than adding a dependency,
+// just define here
+export function zeros(x) {
+  const r = new Array(x);
+  for (let i = 0; i < x; ++i) {
+    r[i] = 0;
+  }
+  return r;
+}
+export function zerosM(x, y) {
+  return zeros(x).map(() => zeros(y));
+}
+
+export function dot(a, b) {
+  let ret = 0;
+  for (let i = 0; i < a.length; ++i) {
+    ret += a[i] * b[i];
+  }
+  return ret;
+}
+
+export function norm2(a) {
+  return Math.sqrt(dot(a, a));
+}
+
+export function scale(ret, value, c?: any) {
+  for (let i = 0; i < value.length; ++i) {
+    ret[i] = value[i] * c;
+  }
+}
+
+export function weightedSum(ret, w1, v1, w2, v2) {
+  for (let j = 0; j < ret.length; ++j) {
+    ret[j] = w1 * v1[j] + w2 * v2[j];
+  }
+}
+
+export function gemv(output, A, x) {
+  for (let i = 0; i < output.length; ++i) {
+    output[i] = dot(A[i], x);
+  }
+}

src/data/utils/venn/fmin/conjugateGradient.ts

@@ -0,0 +1,106 @@
+import { dot, gemv, norm2, scale, weightedSum } from './blas1';
+import { wolfeLineSearch } from './linesearch';
+
+export function conjugateGradient(f, initial, params) {
+  // allocate all memory up front here, keep out of the loop for perfomance
+  // reasons
+  let current = { x: initial.slice(), fx: 0, fxprime: initial.slice() };
+  let next = { x: initial.slice(), fx: 0, fxprime: initial.slice() };
+  const yk = initial.slice();
+  let temp;
+  let a = 1;
+
+  params = params || {};
+  const maxIterations = params.maxIterations || initial.length * 20;
+
+  current.fx = f(current.x, current.fxprime);
+  const pk = current.fxprime.slice();
+  scale(pk, current.fxprime, -1);
+
+  for (let i = 0; i < maxIterations; ++i) {
+    a = wolfeLineSearch(f, pk, current, next, a);
+
+    // todo: history in wrong spot?
+    if (params.history) {
+      params.history.push({
+        x: current.x.slice(),
+        fx: current.fx,
+        fxprime: current.fxprime.slice(),
+        alpha: a,
+      });
+    }
+
+    if (!a) {
+      // faiiled to find point that satifies wolfe conditions.
+      // reset direction for next iteration
+      scale(pk, current.fxprime, -1);
+    } else {
+      // update direction using Polak–Ribiere CG method
+      weightedSum(yk, 1, next.fxprime, -1, current.fxprime);
+
+      const delta_k = dot(current.fxprime, current.fxprime);
+      const beta_k = Math.max(0, dot(yk, next.fxprime) / delta_k);
+
+      weightedSum(pk, beta_k, pk, -1, next.fxprime);
+
+      temp = current;
+      current = next;
+      next = temp;
+    }
+
+    if (norm2(current.fxprime) <= 1e-5) {
+      break;
+    }
+  }
+
+  if (params.history) {
+    params.history.push({
+      x: current.x.slice(),
+      fx: current.fx,
+      fxprime: current.fxprime.slice(),
+      alpha: a,
+    });
+  }
+
+  return current;
+}
+
+/// Solves a system of lienar equations Ax =b for x
+/// using the conjugate gradient method.
+export function conjugateGradientSolve(A, b, x, history?: any) {
+  const r = x.slice();
+  const Ap = x.slice();
+  let rsold;
+  let rsnew;
+  let alpha;
+
+  // r = b - A*x
+  gemv(Ap, A, x);
+  weightedSum(r, 1, b, -1, Ap);
+  const p = r.slice();
+  rsold = dot(r, r);
+
+  for (let i = 0; i < b.length; ++i) {
+    gemv(Ap, A, p);
+    alpha = rsold / dot(p, Ap);
+    if (history) {
+      history.push({ x: x.slice(), p: p.slice(), alpha: alpha });
+    }
+
+    //x=x+alpha*p;
+    weightedSum(x, 1, x, alpha, p);
+
+    // r=r-alpha*Ap;
+    weightedSum(r, 1, r, -alpha, Ap);
+    rsnew = dot(r, r);
+    if (Math.sqrt(rsnew) <= 1e-10) break;
+
+    // p=r+(rsnew/rsold)*p;
+    weightedSum(p, 1, r, rsnew / rsold, p);
+    rsold = rsnew;
+  }
+  if (history) {
+    history.push({ x: x.slice(), p: p.slice(), alpha: alpha });
+  }
+  return x;
+}

src/data/utils/venn/fmin/gradientDescent.ts

@@ -0,0 +1,75 @@
+import { dot, norm2, scale, weightedSum, zeros } from './blas1';
+import { wolfeLineSearch } from './linesearch';
+
+export function gradientDescent(f, initial, params) {
+  params = params || {};
+  const maxIterations = params.maxIterations || initial.length * 100;
+  const learnRate = params.learnRate || 0.001;
+  const current = { x: initial.slice(), fx: 0, fxprime: initial.slice() };
+
+  for (let i = 0; i < maxIterations; ++i) {
+    current.fx = f(current.x, current.fxprime);
+    if (params.history) {
+      params.history.push({
+        x: current.x.slice(),
+        fx: current.fx,
+        fxprime: current.fxprime.slice(),
+      });
+    }
+
+    weightedSum(current.x, 1, current.x, -learnRate, current.fxprime);
+    if (norm2(current.fxprime) <= 1e-5) {
+      break;
+    }
+  }
+
+  return current;
+}
+
+export function gradientDescentLineSearch(f, initial, params) {
+  params = params || {};
+  let current = { x: initial.slice(), fx: 0, fxprime: initial.slice() };
+  let next = { x: initial.slice(), fx: 0, fxprime: initial.slice() };
+  const maxIterations = params.maxIterations || initial.length * 100;
+  let learnRate = params.learnRate || 1;
+  const pk = initial.slice();
+  const c1 = params.c1 || 1e-3;
+  const c2 = params.c2 || 0.1;
+  let temp;
+  let functionCalls = [];
+
+  if (params.history) {
+    // wrap the function call to track linesearch samples
+    const inner = f;
+    f = (x, fxprime) => {
+      functionCalls.push(x.slice());
+      return inner(x, fxprime);
+    };
+  }
+
+  current.fx = f(current.x, current.fxprime);
+  for (let i = 0; i < maxIterations; ++i) {
+    scale(pk, current.fxprime, -1);
+    learnRate = wolfeLineSearch(f, pk, current, next, learnRate, c1, c2);
+
+    if (params.history) {
+      params.history.push({
+        x: current.x.slice(),
+        fx: current.fx,
+        fxprime: current.fxprime.slice(),
+        functionCalls: functionCalls,
+        learnRate: learnRate,
+        alpha: learnRate,
+      });
+      functionCalls = [];
+    }
+
+    temp = current;
+    current = next;
+    next = temp;
+
+    if (learnRate === 0 || norm2(current.fxprime) < 1e-5) break;
+  }
+
+  return current;
+}