BertrandBev / eigen-js

⚡ Eigen-js is a port of the Eigen C++ linear algebra library
https://bertrandbev.github.io/eigen-js/#/
92 stars 15 forks source link

LU decomposition returns wrong values #22

Closed idiotWu closed 3 years ago

idiotWu commented 3 years ago

Problem

eigen.Decompositions.lu() returns wrong values which don't satisfy LU = A.

Steps to Reproduce

Run the following code in demo page...

const M = new eig.Matrix([[2, 1], [1, 4]]);
const lu = eig.Decompositions.lu(M);
return {...lu, "L*U": lu.L.matMul(lu.U)}

...and you will get:

image

While the expected values are:

L = [[1, 0], [0.5, 1]],
U = [[2, 1], [0, 3.5]]
idiotWu commented 3 years ago

Just checked the Eigen documentation for FullPivLU and found this method decomposes matrix A as A = P^-1 * L * U * Q^-1.

Maybe it's better to change to following comment to match the original docs?

https://github.com/BertrandBev/eigen-js/blob/1aa37c8d1d61c83132091f8a4230ba8487b98839/src/classes/Decompositions.js#L20


Edit: tried P^-1 * L * U * Q^-1 but got [[0.25, 1], [0, 0]] which didn't equal to the given matrix [[2, 1], [1, 4]], am I missing something else?

BertrandBev commented 3 years ago

Good point. In practice, L is strictly lower triangular, which means that we need to add the identity matrix to complete the reconstruction

I will update the documentation to

LU decomposition consists in decomposing a square matrix A as a product A = P^-1 (L + I) U * Q^-1 Where U is upper triangular L lower triangular, I is the identity matrix, and P & Q are two permutation matrices

BertrandBev commented 3 years ago

Here we go! https://bertrandbev.github.io/eigen-js/#/decompositions

idiotWu commented 3 years ago

Thanks!