SparseLU regression, not working with eigen 3.4.X for some matrices (but working with 3.3.7 and older)

Summary

I am working with the SparseLU solver. I noticed that it consistently fails to factorize a matrix for which older versions (and other solvers) worked fine.

It is a square matrix (n=1323) with 6 non-zero elements per column. I can confirm that it is a bug and a regression, since:

I have confirmed that the SparseLU solver from the version 3.3.7 is factorizing properly the same matrix.
I have checked that the normal LU solver (FullPivLU) of both versions (3.3.7 + 3.4.X) works on the dense representation of this matrix
I have also validated that the LU solver from the scipy library (in Python) solves the problem.

Even though this problem does not occur with all the big sparse matrices that I am working with, it seems like it is not so uncommon. Indeed, out of roughly 20 tests that I have (each one with different matrices), it happened 3 times.

Environment

Operating System : I first encounter the problem with Windows (MSVC or clang) and I reproduced it with Linux Ubuntu 20. The example that I will provide reproduces the minimal problem for Windows or Linux.
Architecture : x64
Eigen Version : 3.4.90
Compiler Version : g++ (Ubuntu 9.4.0-1ubuntu1~20.04.1) 9.4.0
Compile Flags : does not influence
Vector Extension : I don't know

Minimal Example

To reproduce the problem, the first step was to dump the problematic matrix to reuse it. I have done this using this code. I have copy pasted it here, at the last paragraph.

In my code, I use these functions to save two binary matrices: the sparse one ( that causes the problem) and its dense representation.

// Sparse binary
EigenIO::write_binary_sparse("A_sparse.dat", A);

// Dense binary
const MatrixXd ADense = MatrixXd(A);
EigenIO::write_binary("A_dense.dat", ADense);

Here are the matrices : A_sparse.dat, A_dense.dat

The code that reproduces the problem is the following. This code performs the following tasks

it reads the dense matrix and the sparse matrix.
it performs LU using the dense matrix (takes about 1min) and checks that LU worked to provided the correct inv(P)LUinv(Q) matrix
it performs the SparseLU factorization and checks that it worked using the output of the info() method.

#include <iostream>
#include <fstream>
#include <random>

#include <Eigen/Dense>
#include <Eigen/Sparse>

using std::cout;
using std::endl;
using namespace Eigen;
typedef Eigen::SparseMatrix<double, Eigen::ColMajor> SparseMatCol;

int main() {
        // Read the dense matrix
        MatrixXd A_dense ;
        EigenIO::read_binary("/mnt/c/data/A_dense.dat", A_dense);
        cout << "Dense matrix's shape: " << A_dense.rows() << ", " << A_dense.cols() << endl;

        // LU of the dense matrix
        // This block checks that the LU on the dense matrix works
        // Takes about 20 seconds on my machine
        const int row = 1323;
        const int col = 1323;
        Eigen::FullPivLU<MatrixXd> lu(A_dense);
        MatrixXd l = Matrix<double, row, col>::Identity();
        l.block<row,col>(0,0).triangularView<StrictlyLower>() = lu.matrixLU();
        MatrixXd u = lu.matrixLU().triangularView<Upper>();
        MatrixXd A_LU = lu.permutationP().inverse() * l * u * lu.permutationQ().inverse();
        cout << "[Dense LU] LU worked ? " << A_dense.isApprox(A_LU) << endl;


        // Read the sparse matrix
        SparseMatCol A_sparse ;
        EigenIO::read_binary_sparse("/mnt/c/data/A_sparse.dat", A_sparse);
        cout << "Sparse matrix's shape: " << A_sparse.rows() << ", " << A_sparse.cols() << endl;

        // LU of the sparse matrix
        // This block runs the SparseLU solver
        // This will return the info = 1 (NumericalIssue) with eigen 3.4.90
        Eigen::SparseLU<SparseMatCol, Eigen::COLAMDOrdering<SparseMatCol::StorageIndex>> solver;
        solver.analyzePattern(A_sparse);
        solver.factorize(A_sparse);
        cout << "[Sparse LU] Solver info: " << solver.info() << endl;
        if (solver.info())
                cout << solver.lastErrorMessage() << endl;

        return 0;
}

Steps to reproduce

On my machine, I have two versions of Eigen:

Eigen 3.3.7 under /usr/include/eigen3
Eigen 3.4.0 under /home/abricq/dev/eigen

To reproduce the problem, compile with g++ eigen_matrix_io.cpp -I/home/abricq/dev/eigen

To show that it is a regression, compile with g++ eigen_matrix_io.cpp -I/usr/include/eigen3. The problem will not occur.

What is the current bug behavior?

After running the code while linking Eigen 3.4.90, the output of the code is:

Dense matrix's shape: 1323, 1323
[Dense LU] LU worked ? 1
Sparse matrix's shape: 1323, 1323
[Sparse LU] Solver info: 1
THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT 1323

I was able to link this error at this line: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/SparseLU/SparseLU.h#L757

I have very little understanding of your code base, it seems to me that it is unable to find a good pivot, at this point: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/SparseLU/SparseLU_pivotL.h#L96.

I looked if these files have changed in the last versions and I found that they did not. Therefore I have no idea of what could be causing this issue.

What is the expected correct behavior?

The expected behavior is the one obtained using the previous version. The output looks now like this when linking Eigen 3.3.7

Dense matrix's shape: 1323, 1323
[Dense LU] LU worked ? 1
Sparse matrix's shape: 1323, 1323
[Sparse LU] Solver info: 1

Please, let me know if I can provide anymore information. Thanks in advance for the help.

IO code

Here is the code that I use to read and write to file.

// found at https://stackoverflow.com/a/25389481/11927397
namespace EigenIO {
        template<class Matrix>
        inline void write_binary(const std::string& filename, const Matrix& matrix){
                std::ofstream out(filename, std::ios::out | std::ios::binary | std::ios::trunc);
                if(out.is_open()) {
                        typename Matrix::Index rows=matrix.rows(), cols=matrix.cols();
                        out.write(reinterpret_cast<char*>(&rows), sizeof(typename Matrix::Index));
                        out.write(reinterpret_cast<char*>(&cols), sizeof(typename Matrix::Index));
                        out.write(reinterpret_cast<const char*>(matrix.data()),
                                        rows*cols*static_cast<typename Matrix::Index>(sizeof(typename Matrix::Scalar)) );
                        out.close();
                }
                else {
                        std::cout << "Can not write to file: " << filename << std::endl;
                }
        }

        template<class Matrix>
        inline void read_binary(const std::string& filename, Matrix& matrix){
                std::ifstream in(filename, std::ios::in | std::ios::binary);
                if (in.is_open()) {
                        typename Matrix::Index rows=0, cols=0;
                        in.read(reinterpret_cast<char*>(&rows),sizeof(typename Matrix::Index));
                        in.read(reinterpret_cast<char*>(&cols),sizeof(typename Matrix::Index));
                        matrix.resize(rows, cols);
                        in.read(reinterpret_cast<char*>(matrix.data()),
                                        rows*cols*static_cast<typename Matrix::Index>(sizeof(typename Matrix::Scalar)) );
                        in.close();
                }
                else {
                        std::cout << "Can not open binary matrix file: " << filename << std::endl;
                }
        }

        // https://scicomp.stackexchange.com/a/21438
        template <class SparseMatrix>
        inline void write_binary_sparse(const std::string& filename, const SparseMatrix& matrix) {
                assert(matrix.isCompressed() == true);
                std::ofstream out(filename, std::ios::binary | std::ios::out | std::ios::trunc);
                if(out.is_open())
                {
                        typename SparseMatrix::Index rows, cols, nnzs, outS, innS;
                        rows = matrix.rows()     ;
                        cols = matrix.cols()     ;
                        nnzs = matrix.nonZeros() ;
                        outS = matrix.outerSize();
                        innS = matrix.innerSize();

                        out.write(reinterpret_cast<char*>(&rows), sizeof(typename SparseMatrix::Index));
                        out.write(reinterpret_cast<char*>(&cols), sizeof(typename SparseMatrix::Index));
                        out.write(reinterpret_cast<char*>(&nnzs), sizeof(typename SparseMatrix::Index));
                        out.write(reinterpret_cast<char*>(&outS), sizeof(typename SparseMatrix::Index));
                        out.write(reinterpret_cast<char*>(&innS), sizeof(typename SparseMatrix::Index));

                        typename SparseMatrix::Index sizeIndexS =
                                static_cast<typename SparseMatrix::Index>(sizeof(typename SparseMatrix::StorageIndex));
                        typename SparseMatrix::Index sizeScalar =
                                static_cast<typename SparseMatrix::Index>(sizeof(typename SparseMatrix::Scalar));
                        out.write(reinterpret_cast<const char*>(matrix.valuePtr()),       sizeScalar * nnzs);
                        out.write(reinterpret_cast<const char*>(matrix.outerIndexPtr()),  sizeIndexS  * outS);
                        out.write(reinterpret_cast<const char*>(matrix.innerIndexPtr()),  sizeIndexS  * nnzs);
                        out.close();
                }
                else {
                        std::cout << "Can not write to file: " << filename << std::endl;
                }
        }

        template <class SparseMatrix>
        inline void read_binary_sparse(const std::string& filename, SparseMatrix& matrix) {
                std::ifstream in(filename, std::ios::binary | std::ios::in);
                if(in.is_open()) {
                        typename SparseMatrix::Index rows, cols, nnz, inSz, outSz;
                        typename SparseMatrix::Index sizeScalar =
                                static_cast<typename SparseMatrix::Index>(sizeof(typename SparseMatrix::Scalar      ));
                        typename SparseMatrix::Index sizeIndex  =
                                static_cast<typename SparseMatrix::Index>(sizeof(typename SparseMatrix::Index       ));
                        typename SparseMatrix::Index sizeIndexS =
                                static_cast<typename SparseMatrix::Index>(sizeof(typename SparseMatrix::StorageIndex));
                        in.read(reinterpret_cast<char*>(&rows ), sizeIndex);
                        in.read(reinterpret_cast<char*>(&cols ), sizeIndex);
                        in.read(reinterpret_cast<char*>(&nnz  ), sizeIndex);
                        in.read(reinterpret_cast<char*>(&outSz), sizeIndex);
                        in.read(reinterpret_cast<char*>(&inSz ), sizeIndex);

                        matrix.resize(rows, cols);
                        matrix.makeCompressed();
                        matrix.resizeNonZeros(nnz);

                        in.read(reinterpret_cast<char*>(matrix.valuePtr())     , sizeScalar * nnz  );
                        in.read(reinterpret_cast<char*>(matrix.outerIndexPtr()), sizeIndexS * outSz);
                        in.read(reinterpret_cast<char*>(matrix.innerIndexPtr()), sizeIndexS * nnz );

                        matrix.finalize();
                        in.close();
                } // file is open
                else {
                        std::cout << "Can not open binary sparse matrix file: " << filename << std::endl;
                }
        }
}