Unnecessary allocation in TriDiagonalize

Summary

While computing eigenvalues on a SelfAdjoint matrix through the SelfAdjointEigenSolver::compute() method, there should be no memory allocation as everything should be pre-allocated by the constructor for dynamic-size matrices. However, there are some allocations in tridiagonalization_inplace_selector method (line 444 of Tridiagonalization.h).

This has a performance impact as my program should run thousands of eigen value computations on a matrix of constant size (constant at runtime).

Environment

  • Operating System : Windows
  • Architecture : x64
  • Eigen Version : 3.3.9
  • Compiler Version : Gcc9.3.0
  • Compile Flags : -O3
  • Vector Extension : SSE/AVX/NEON

Minimal Example

#define EIGEN_RUNTIME_NO_MALLOC // Define this symbol to enable runtime tests for allocations
#include <Eigen/Dense>

int main(int argc, char *argv[])
{
    //At the beginning, allocation can be done
    Eigen::internal::set_is_malloc_allowed(true);
    Eigen::MatrixXd A(4, 4);
    A.setIdentity();
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(A);

    //After this line, there should be no allocation
    Eigen::internal::set_is_malloc_allowed(false);
    for (int i = 0; i < 10; i++)
    {
        es.compute(A);
    }
    return 0;
}

Steps to reproduce

Run the above code sample.

What is the current bug behavior?

Throws an exception:

assertion failed: is_malloc_allowed() && "heap allocation is forbidden (EIGEN_RUNTIME_NO_MALLOC is defined and g_is_malloc_allowed is false)", file .../eigen3/Eigen/src/Core/util/Memory.h, line 143

What is the expected correct behavior?

According to the compute() method documentation:

This method reuses the memory in the SelfAdjointEigenSolver object that was allocated when the object was constructed, if the size of the matrix does not change.

So, there should be no allocation and therefore no exception.

Relevant logs

Stack trace:

Eigen::internal::tridiagonalization_inplace_selector<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, false>::run<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> & mat, Eigen::Matrix<double, -1, 1, 0, -1, 1> & diag, Eigen::Matrix<double, -1, 1, 0, -1, 1> & subdiag, bool extractQ) (...\lib\include\eigen3\Eigen\src\Eigenvalues\Tridiagonalization.h:444)
Eigen::internal::tridiagonalization_inplace<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> & mat, Eigen::Matrix<double, -1, 1, 0, -1, 1> & diag, Eigen::Matrix<double, -1, 1, 0, -1, 1> & subdiag, bool extractQ) (...\lib\include\eigen3\Eigen\src\Eigenvalues\Tridiagonalization.h:430)
Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> >::compute<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, -1, -1, 0, -1, -1> > * const this, const Eigen::EigenBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > & a_matrix, int options) (...\lib\include\eigen3\Eigen\src\Eigenvalues\SelfAdjointEigenSolver.h:437)
main(int argc, char ** argv) (...\main.cpp:18)
  • Have a plan to fix this issue. The fix is too complex for me to propose:
  1. add additionnal member CoeffVectorType hCoeffs to SelfAdjointEigenSolver class and pass it as an additional argument to tridiagonalization_inplace_selector::run method => this I can do
  2. Fix the allocation that comes from the following line in tridiagonalization_inplace_selector::run => this I don't know how to fix HouseholderSequenceType(mat, hcoeffs.conjugate()).setLength(mat.rows() - 1).setShift(1);
Edited by Marc Teyssier