nmwsharp / geometry-central

Applied 3D geometry in C++, with a focus on surface meshes.
https://geometry-central.net
MIT License
1.01k stars 141 forks source link

Fix SuiteSparse call for underdetermined systems #178

Closed zoemarschner closed 2 weeks ago

zoemarschner commented 3 weeks ago

Currently, solving underdetermined systems using SuiteSparse by calling solve does not work. For example,

#include "geometrycentral/numerical/linear_solvers.h"

Eigen::Matrix<double, 2, 3> A;
A << 1, 0, 1, 0, 1, 0 ;
Vector<double> b(2);
b << 1, 2;
SparseMatrix<double> sparseA = A.sparseView();

Vector<double> sol = solve(sparseA, b);

std::cout << sol << std::endl;

will print out 1 1 1, but the correct answer is 0.5 2 0.5. This is because SuiteSparse's QR decomposition includes a permutation matrix (i.e., E such that QR=AE), but geometry central does not properly account for this in the call to SuiteSparseQR_solve. (Thank you @MarkGillespie for helping to find the bug!)

nmwsharp commented 2 weeks ago

Looks great, thank you so much for this! I haven't tested this code, but the problem makes sense and the fix looks reasonable so I'm going to merge it.

Also, I love that this removes a // TODO what is this E doing here? comment 😂