Solving underdetermined systems with SparseQR Factorization about csparse.net HOT 5 CLOSED

Scanning over the document I couldn't find a reference to QR, neither to the graph Laplacian.

Regarding rank deficient systems in general: the pseudo-inverse may be used to get the best solution (in a least squares sense) and the QR factorization can be used to find it, though you have to be careful when dealing with sparse matrices. Here's some Matlab code [improved qrginv]:

function ginv = IMqrginv(A)
    [Q,R,P] = qr(A);
    r=sum(any(abs(R)>1e-5,2));
    Q = Q(:, 1:r);
    R = R(1:r, :);
    ginv = P*(R'/(R*R')*Q');
end

EDIT Some remarks:

• Of course, you don't want to compute the last matrix product, since it will be dense. You'd have to follow the solving procedure of the SparseQR class.
• SparseQR doesn't provide the full Q matrix (and you should never compute it, see here), but you don't need it here. Just apply the first r Househoulder vectors (r = rank of the matrix, in case of the Laplacian of a connected graph = n-1).
• Use the code from Factorization members wiki page to get the QR factorization data.

Alternatively, if you're interested in only one particular solution, you could try the MINRES iterative solver, which can handle singular systems.

from csparse.net.

Thanks for the reply! The equation involving the Laplacian matrix in the document is Poisson's equation Lc.phi=d, at the fourth page, left column.

We try to retrieve a function whose gradient is equal to a given value (the eikonal equation), this is thus natural that the solution is given up to a constant. Since the matrix is semi-definite positive, I tried a Tikhonov regularization by solving (Lc+epsilon*Identity)*phi=d together with a Cholesky Factorization. It works very well for small graphs, but not for large ones, I am still investigating why. I wonder if you have ever tried this kind of regularization scheme with CSparse.NET?

Best,
Romain

from csparse.net.

This seems like a valid approach, though I never had to use this kind of regularization. Obviously, depending on your choice of epsilon, the system gets ill-conditioned and errors in the solution process might amplify in higher dimensions. Additionally, mesh quality always plays a role, so if the solution fails, you might check for sharp angles or other defects.

from csparse.net.

Hello,
thanks for your reply. I performed some tests with 10k vertices (which is enough for my application), and with epsilon=1e-6. The LU factorization gives satisfying results with the Tikhonov regularization technique. I guess that some numerical errors mesh processing and laplacian matrix construction make my matrix slightly non-symmetric and make the Cholesky factorization algorithm unstable (for the record, the mesh is a Delaunay triangulation and the Laplacian matrix is supposed to by semi-definite negative, so a small shift and multiplication by -1 should make the problem sdp). The performance with LU factorization is enough for what I want to achieve, I will try to implement a solution inspired by your first reply for the sake of completeness.

from csparse.net.

FYI, I've pushed MINRES to wo80/csparse-extensions@185dd0d . Not tested in practice, though.

Here's an example how to use it: MinResTest.cs

from csparse.net.