How would I implement MATLAB's "eig(A, B)" function in JavaScript

Question

According to MATLAB's documentation:

[V,D] = eig(A,B) returns diagonal matrix D of generalized eigenvalues and full matrix V whose columns are the corresponding right eigenvectors, so that A*V = B*V*D.

As I read available source code (It seems all implementations I've looked at Octave, R, Scipy) lead to LAPACK's DGGEV routine, which isn't available from JavaScript.

Rabbit holes include Use Emscripten with Fortran: LAPACK binding and learning enough Fortran and linear algebra to do it myself.

Anyone know of something more accessible?

Answer 1

I ended up finding success using the Eigen library, combined with Emscripten.

Right now, my test code is hard-coded to 5x5 matrices, but that's just a matter of template arguments.

I'm passing data to and from the function by using row major 1D arrays.

The code looks something like:

#include <Eigen/Eigenvalues>
#typedef double ArrayMat5d[25];
#typedef double ArrayVec5d[5];
#typedef Eigen::Matrix<double, 5, 5, Eigen::RowMajor> Matrix5dR;
#typedef Eigen::Matrix<double, 5, 1> Vector5d;

extern "C" void eig(const ArrayMat5d A, const ArrayMat5d B,
        ArrayMat5d V, ArrayVec5d D) {
    Eigen::Map<const Matrix5dR> a(A);
    Eigen::Map<const Matrix5dR> b(B);
    const Eigen::GeneralizedSelfAdjointEigenSolver<Matrix5dR> solver(a, b);
    Eigen::Map<Matrix5dR> v(V);
    Eigen::Map<Vector5d> d(D);
    v = solver.eigenvectors();
    d = solver.eigenvalues();
}

And I'm compiling the code using:

emcc -I /usr/include/eigen3 -O2 -o eig.js -s "DISABLE_EXCEPTION_CATCHING = 1" \ 
-s "NO_FILESYSTEM = 1" -s "NO_BROWSER = 1" -s "EXPORTED_FUNCTIONS = ['_eig']" \
-s "NO_EXIT_RUNTIME = 1" eig.cpp

From the JavaScript side:

// builds reference to eig function with argument type checking
var eig = Module.cwrap('eig', null, ['number', 'number', 'number', 'number']);

// sets up the two matrices
var P = new Float64Array([ 92.31360, 11.75040, -15.84640, -21.88800, -0.83200, 11.75040, 15.76960, -4.37760, -0.83200, 2.11200, -15.84640, -4.37760, 4.24960, 2.11200, -1.15200, -21.88800, -0.83200, 2.11200, 15.04000, -1.44000, -0.83200, 2.11200, -1.15200, -1.44000, -2.24000 ]);
var Q = new Float64Array([ 60.16, -2.88, 0.0, 0.0, 0.0, -2.88, 17.28, -2.88, 0.0, 0.0, 0.0, -2.88, 8.96, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 ]);

// allocates memory for input and output matrices
var matLength = 25;
var vecLength = 5;
var matSize = matLength * P.BYTES_PER_ELEMENT;
var vecSize = vecLength * P.BYTES_PER_ELEMENT;

var Pptr = Module._malloc(matSize);
var Qptr = Module._malloc(matSize);
var Vptr = Module._malloc(matSize);
var Dptr = Module._malloc(vecSize);

// gets references to Emscripten heap
var Pheap = new Uint8Array(Module.HEAPU8.buffer, Pptr, matSize);
var Qheap = new Uint8Array(Module.HEAPU8.buffer, Qptr, matSize);
var Vheap = new Uint8Array(Module.HEAPU8.buffer, Vptr, matSize);
var Dheap = new Uint8Array(Module.HEAPU8.buffer, Dptr, vecSize);

// copies input matrices into Emscripten heap
Pheap.set(new Uint8Array(P.buffer));
Qheap.set(new Uint8Array(Q.buffer));

// calls the function (finally!)
eig(Pheap.byteOffset, Qheap.byteOffset, Vheap.byteOffset, Dheap.byteOffset);

// Gets double views into Emscripten heap containing results
var Vresult = new Float64Array(Vheap.buffer, Vheap.byteOffset, P.length);
var Dresult = new Float64Array(Dheap.buffer, Dheap.byteOffset, vecLength);

console.log(Vresult);
console.log(Dresult);

// Frees up allocated memory
Module._free(Pheap.byteOffset);
Module._free(Qheap.byteOffset);
Module._free(Vheap.byteOffset);
Module._free(Dheap.byteOffset);

The whole thing works quite well. At level -O2, I'm getting times of about 800 ms to run 10000 iterations, and the results exactly match my original C++ test code. (It's just about exactly 10x slower at -O0.)

Now to finish that ellipse fit!