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?

Try this http://www.numericjs.com/workshop.php?link=bf9eb195e289cf58fe226b11ff9e70d3360865f1aa1565d23990ad12673fd7de — Axalix, May 16 '15 at 00:05
or [sylvester.js](http://sylvester.jcoglan.com/) both libs include implementations of EVD but one can just take the generic algorithm and implement it in js as is — Nikos M., May 16 '15 at 00:26
@NikosM. I don't see any eigenvalue decomposition (let alone generalized eigen decomposition) in sylvester.js –where is it? I'm not sure numeric.js supports generalized eigen decomposition either (i.e., `eig(A,B)` rather than just `eig(A)`. — horchler, May 16 '15 at 04:19
@horchler, yeap my mistake, it has routines for singular and determinant (had used it some time back, did not remember correctly) — Nikos M., May 16 '15 at 04:31
@NikosM. yeah, could be useful primitives should I try to implement it myself, but my linear algebra knowledge would need a lot of holes filled first. I'm hoping there's a solution that's, if anything, heavier on the code and lighter on the math. — squidpickles, May 17 '15 at 05:24
@slugchewer, i c ant find the js lib i had used, which implemented eigenvalues, so i'll add a lib from a past project mine in php [Simulacra](https://github.com/foo123/Simulacra) which implements matrix/linear algebra algorithms like SVD, LU decomposition and more, it is easy to port it to javascript and modify as needed or add new features — Nikos M., May 17 '15 at 10:37
Computing generalised eigenvalues/eigenvectors in the general case is rather tricky. There is a very important class of generalised eigenvector problems where A is symmetric (hermitian) and B is symmetric (hermitian) positive definite. This latter case is easier to solve and can be reformulated in terms of standard eigenvalue/eigenvector problems. Please explain what type of problems (or the origin of them) you are interested in solving. — Stefano M, May 17 '15 at 22:10
@StefanoM this is part of an ellipse fitting algorithm - http://www.mathworks.com/matlabcentral/fileexchange/22683-ellipse-fit--taubin-method- - that I'm porting to JavaScript. — squidpickles, May 17 '15 at 22:51
It seems that you are in the special case, confirmed also from the fact that you are using `Eigen::GeneralizedSelfAdjointEigenSolver`. This means that you can reduce the generalized eigenvalue problem to a regular one by computing the Cholesky factorisation of B. — Stefano M, May 18 '15 at 15:22
@StefanoM the truth is, I just learned what an eigenvector was in context of porting an ellipse fit algorithm from MATLAB to JavaScript. It's possible my use of `Eigen::GeneralizedSelfAdjointEigenSolver` isn't correct; I'll keep your comments in mind if I run into problems with this implementation. In any case, the approach seems to be working. Thanks for the help! — squidpickles, May 18 '15 at 17:52
Don't worry: using `Eigen::GeneralizedSelfAdjointEigenSolver` for a non self-adjoint problem would completely blow up your solution (and you would for sure notice it.) Let me more explicit: your choice of the self-adjoint solver looks appropriate. Sometimes beginner's luck outperforms Murphy's law! — Stefano M, May 18 '15 at 18:25

squidpickles · Accepted Answer · 2015-05-18T17:39:22.583

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!

The details provided in this answer are invaluable, thanks! – mucaho Dec 07 '15 at 03:17 — mucaho, Dec 07 '15 at 03:17

score 0 · Answer 2 · answered Jun 28 '20 at 14:10

I made a javascript module to calculate eigenvalues, eigenvectors, or RREF(Reduced Row Echelon Form). You can install it using npm i @ahmaddynugroho/eig

Example:

const e = require('@ahmaddynugroho/eig')

const A = [[7,3],[3,-1]]
const ans = e.eig(A)
console.log(ans)

That will log the ans object as follow:

{
  eigval: [ 8.000000000000009, -2.0000000000000004 ],
  eigvec: [
    [ -0.9486832980505129, 0.31622776601684044 ],
    [  0.3162277660168379, 0.9486832980505138  ]
  ]
}

That's it, you just need to pass a matrix as a parameter for eig(). Btw, here's the documentation

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

2 Answers2