Solving underdetermined systems with SparseQR Factorization

[romainmesnil]

Hello, First, thanks a lot for this library, I'm glad that there is finally a good native C# open-source sparse matrix library! I have a beginner's question regarding the solution of underdetermined systems.

I am currently doint a full C# implementation of the Heat Method by Keenan Crane et al., which requires to solve a linear system Ax=y, where A is a square matrix, which is not full-rank (A is actually the Laplacian matrix of a graph, so the vector (1,1,1,1,1,1,1...,1) is in the null space of A). I would like to solve the system in a least square sense, and I have read that QR factorization is well-suited for this purpose.

I applied the same code as the sparseQR example, but with a laplacian matrix. The result of Solve is a vector proportional (1,1,1,1,1,1,1...,1) with a very large norm (each component of the vector has a value superior to 1 million, while the yvector has components inferior to 10.

Is there a specific trick to solve this in the least square sense with the SparseQR Factorization?

Best, Romain

[wo80]

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.

[romainmesnil]

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

[wo80]

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.

[romainmesnil]

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.

[wo80]

FYI, I've pushed MINRES to https://github.com/wo80/csparse-extensions/commit/185dd0d351cc4a4970bf551e3c0b184699e77911 . Not tested in practice, though.

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