Eigen: Obtain the kernel of a sparse matrix

Question

Given a sparse matrix A and a vector b, I would like to obtain a solution x to the equation A * x = b as well as the kernel of A.

One possibility is to convert A to a dense representation.

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/SparseQR>

int main()
{
    // This is a toy problem. My actual matrix
    // is of course bigger and sparser.
    Eigen::SparseMatrix<double> A(2,2);
    A.insert(0,0) = 1;
    A.insert(0,1) = 2;
    A.insert(1,0) = 4;
    A.insert(1,1) = 8;
    A.makeCompressed();

    Eigen::Vector2d b;
    b << 3, 12;

    Eigen::SparseQR<Eigen::SparseMatrix<double>,
                    Eigen::COLAMDOrdering<int> > solver;
    solver.compute(A);
    std::cout << "Solution:\n" << solver.solve(b) << std::endl;

    Eigen::Matrix2d A_dense(A);
    std::cout << "Kernel:\n" << A_dense.fullPivLu().kernel() << std::endl;
    return 0;
}

Is it possible to do the same directly in the sparse representation? I could not find a function kernel() anywhere except in FullPivLu.

To lazy to implement it myself, but what should work is computing the SparseQR of `A.transpose()` and obtain the last `qr.rank()` columns of `Q`, by multiplying `qr.matrixQ()` by the last columns of a sparse (or dense) identity matrix. — chtz, Feb 19 '19 at 18:41
@chtz Wouldn't you like to turn your comment into an answer? If there would be a little explanation as to why it works, I'd be happy to accept it. — Dominik Mokriš, Feb 27 '19 at 20:10

Hongkai Dai · Accepted Answer · 2019-05-22T23:19:40.717

I think @chtz's answer is almost correct, except we need to take the last A.cols() - qr.rank() columns. Here is a mathematical derivation.

Say we do a QR decomposition of your matrix Aᵀ as

Aᵀ * P = [Q₁ Q₂] * [R; 0] = Q₁ * R

where P is the permutation matrix, thus

Aᵀ = Q₁ * R * P⁻¹.

We can see that Range(Aᵀ) = Range(Q₁ * R * P⁻¹) = Range(Q₁) (because both P and R are invertible).

Since Aᵀ and Q₁ have the same range space, this implies that A and Q₁ᵀ will also have the same null space, namely Null(A) = Null(Q₁ᵀ). (Here we use the property that Range(M) and Null(Mᵀ) are complements to each other for any matrix M, hence Null(A) = complement(Range(Aᵀ)) = complement(Range(Q₁)) = Null(Q₁ᵀ)).

On the other hand, since the matrix [Q₁ Q₂] is orthonormal, Null(Q₁ᵀ) = Range(Q₂), thus Null(A) = Range(Q₂), i.e., kernal(A) = Q₂.

Since Q₂ is the right A.cols() - qr.rank() columns, you could call rightCols(A.cols() - qr.rank()) to retrieve the kernal of A.

For more information on kernal space, you could refer to https://en.wikipedia.org/wiki/Kernel_(linear_algebra)

Eigen: Obtain the kernel of a sparse matrix

1 Answers1