CODEWITHSUNDEEP
c++linear-algebraeigen
dari
dari

Reputation: 2455

Fastest way to solve small homogeneous linear systems with eigen

I need to solve a lot of small (n=4) homogeneous linear systems of the form Ax=0 with A being a singular matrix. I'm currently using the following code:

void solve(const matrix_t& A, vector_t& x){
    auto svd = A.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV);
    auto V = svd.matrixV();
    x = V.col( A.rows() - 1 );
    x.normalize();
}

Are there any faster ways to do that?

Upvotes: 0

Views: 1643

Answers (1)

dari
dari

Reputation: 2455

The fastest way to obtain the nullspace of a general matrix with eigen is by using its LU decomposition. In praxis I'm using the Householder QR decomposition instead of LU, because it appears to be more stable when the input matrix isn't perfectly singular. QR is still alot faster than the SVD proposed in the question and gives very similiar results for my problems. A benchmark of the different eigen decompositions can be found here: https://eigen.tuxfamily.org/dox/group__DenseDecompositionBenchmark.html

Code for computing the nullspace with LU, QR and SVD (note: x.normalize() is not required but helpful to compare the solutions):

template<typename matrix_t, typename vector_t>
void solveNullspaceLU(const matrix_t& A, vector_t& x){
    x = A.fullPivLu().kernel();
    x.normalize();
}

template<typename matrix_t, typename vector_t>
void solveNullspaceQR(const matrix_t& A, vector_t& x){
    auto qr = A.transpose().colPivHouseholderQr();
    matrix_t Q = qr.householderQ();
    x = Q.col(A.rows() - 1);
    x.normalize();
}

template<typename matrix_t, typename vector_t>
void solveNullspaceSVD(const matrix_t& A, vector_t& x){
    x = A.jacobiSvd(Eigen::ComputeFullV).matrixV().col( A.rows() - 1 );
    x.normalize();
}

Upvotes: 7

Related Questions