Tail call optimization with eigen

Question

I'm trying to make the compiler tail-call optimize eigen-expressions that are passed recursively. The following will optimize using GCC -O2:


inline double recursive(const double &A,
                        const double &b,
                        const int &N)
{
    if (N <= 0)
        return b;
    return recursive(A, b/A, N - 1);
}

int main()
{
    double A = 1;
    double b = 2;
    double res2 = recursive(A, b, 10);
}

The following does not:

#include <iostream>
#include <Eigen/Dense>
using Mat = Eigen::MatrixXd;
using Vec = Eigen::VectorXd;
using Eigen::MatrixBase;


inline Vec recursive(const Mat &A,
                        const Vec &b,
                        const int &N)
{
    if (N <= 0)
        return b;
    return recursive(A, A.ldlt().solve(b), N - 1);
}

int main()
{
    Mat A(2, 2);
    A(0, 0) = 10.0;
    A(1, 1) = 10.0;

    Vec b(2);
    b(0) = 1.0;
    b(1) = 5.0;

    Vec res(2);

    Vec res2 = recursive(A, b, 10);
}

Even when I pass the computation as sol = A.ldlt().solve(b).Eval(); it is unable to perform optimization. Why?

Think about where temporaries are created and destructed. Try your first example, but replace the type of `b` with a class with a non-trivial destructor. — chtz, Nov 01 '21 at 09:52
It's not a tail call. Destructors invisibly turn apparent tail calls into non-tail calls, since they need to execute afterwards. — molbdnilo, Nov 01 '21 at 09:57
To elaborate: `A.ldlt().solve(b)` creates two temporaries that need to be destroyed after the recursion returns. Thus, the recursion is not in tail position. Storing the result in a local variable doesn't help, since that also needs to be destroyed after recursing. — molbdnilo, Nov 01 '21 at 10:07
I see, so then I assume there is no way to make this work with the Eigen-datatypes? — Anon232, Nov 01 '21 at 10:09
You could just write it as a loop yourself (this would also avoid calculating `A.ldlt()` each step). Also, even for trivial types, the C++ standard does not require compilers to optimize away tail-recursion (as far as I know ...), so just avoid it if possible. — chtz, Nov 01 '21 at 10:27
The Matrix-class inherits a protected destructor from MatrixBase in Eigen. This destructor is default, so then I assume that the protected declaration is the root cause of the problem? — Anon232, Nov 01 '21 at 11:09

score 1 · Answer 1 · answered Nov 03 '21 at 13:21

As others have commented, the code as it stands cannot be tail-recursive because destructors need to run after the recursion returns. You can rewrite the code to achieve this, though:

static void recursive_impl(const Mat &A,
                          Vec &b,
                          int N)
{
    if (N <= 0)
        return b;
    b = A.ldlt().solve(b);
    recursive_impl(A, b, N - 1);
}
inline Vec recursive(const Mat &A,
                        const Vec &b,
                        const int N)
{
  Vec rtrn = b;
  recursive_impl(a, rtrn, N);
  return rtrn;
}

Now the recursive_impl function has no temporary objects that exist when it does its tail call. ldlt().solve() may allocate memory temporarily, but it is released before the tail call.

Note that it is pointless and often detrimental to pass arguments to builtin types like int and double as const &. The reference is a pointer in disguise. The object to which the pointer points has to live somewhere, presumably on the stack. This prevents tail call optimization. Even your first example with the double recursive only works because the compiler is clever enough to figure out that the reference is unnecessary.

Tail call optimization with eigen

1 Answers1