Hi,

I am new here, so let me start by saying thanks for producing and providing Eigen. It's an amazing library.

To my issue: I am trying to write a fairly generic implementation of an algorithm. At a few points in my algorithm I need to solve linear problems. To keep things DRY, I decided to wrap the corresponding series of Eigen-calls into a function. Furthermore, to also avoid need-less copies I thought it would be a good idea to return the corresponding expression template instead of a plain matrix object. Unfortunately, this fails due to a violated assertion inside Eigen:

When compiled with Clang 3.4:

Code: Select all

fails: /usr/include/eigen3/Eigen/src/QR/ColPivHouseholderQR.h:364: Index Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >::nonzeroPivots() const [MatrixType = Eigen::Matrix<double, 2, 2, 0, 2, 2>]: Assertion `m_isInitialized && "ColPivHouseholderQR is not initialized."' failed.
Aborted

And when compiled with GCC 4.8.1:

Code: Select all

fails: /usr/include/eigen3/Eigen/src/Core/MapBase.h:148: Eigen::MapBase<Derived, 0>::MapBase(Eigen::MapBase<Derived, 0>::PointerType, Eigen::MapBase<Derived, 0>::Index, Eigen::MapBase<Derived, 0>::Index) [with Derived = Eigen::Block<const Eigen::Matrix<double, 2, 2>, -1, 1, false>; Eigen::MapBase<Derived, 0>::PointerType = const double*; Eigen::MapBase<Derived, 0>::Index = long int]: Assertion `(dataPtr == 0) || ( nbRows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == nbRows) && nbCols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == nbCols))' failed.

The above assertion is generated by the following example code:

Code: Select all

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

using namespace Eigen;

template <class DerivedA, class DerivedB>
auto solve(const MatrixBase<DerivedA> &A, const MatrixBase<DerivedB> &B)
    -> decltype(A.colPivHouseholderQr().solve(B)) {
    return A.colPivHouseholderQr().solve(B);
}

int main() {
    Matrix2d A;
    Vector2d b, c;

    A << 1, 2,
         4, 3;
    b << 5, 10;

    c = solve(A, b);

    std::cout << c << "\n";
}

It seems to me, that the result of `solve()` has some internal state, and does not quite like losing its scope before you assign it to something.

A GDB back-trace of the GCC compiled version follows in the end.

A quite hack-ish way to remedy this problem is to transpose the result twice. I.e. simply replace the `solve` function above by the following piece of code:

Code: Select all

template <class DerivedA, class DerivedB>
auto solve(const MatrixBase<DerivedA> &A, const MatrixBase<DerivedB> &B)
    -> decltype(A.colPivHouseholderQr().solve(B).transpose().transpose()) {
    return A.colPivHouseholderQr().solve(B).transpose().transpose();
}

Now, my question is: Am I doing something wrong, or is this an issue within Eigen? And why does the double `transpose` fix it? Does it contain an implicit `eval` or is there something else to it?

Btw, I've tried it with `fullPivLu` as well. The outcome is the same.

Cheers

Andreas

EDIT: This is with Eigen version 3.2.2

---

The promised back-trace

Code: Select all

(gdb) bt
#0  0x00007ffff7241849 in __GI_raise (sig=sig@entry=6)
    at ../nptl/sysdeps/unix/sysv/linux/raise.c:56
#1  0x00007ffff7242cd8 in __GI_abort () at abort.c:89
#2  0x00007ffff723a616 in __assert_fail_base (
    fmt=0x7ffff7374f38 "%s%s%s:%u: %s%sAssertion `%s' failed.\n%n",
    assertion=assertion@entry=0x41f630 "(dataPtr == 0) || ( nbRows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == nbRows) && nbCols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == nbCols))",
    file=file@entry=0x41f3d8 "/usr/include/eigen3/Eigen/src/Core/MapBase.h", line=line@entry=148,
    function=function@entry=0x427000 <Eigen::MapBase<Eigen::Block<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false>, 0>::MapBase(double const*, long, long)::__PRETTY_FUNCTION__> "Eigen::MapBase<Derived, 0>::MapBase(Eigen::MapBase<Derived, 0>::PointerType, Eigen::MapBase<Derived, 0>::Index, Eigen::MapBase<Derived, 0>::Index) [with Derived = Eigen::Block<const Eigen::Matrix<doub"...) at assert.c:92
#3  0x00007ffff723a6c2 in __GI___assert_fail (
    assertion=0x41f630 "(dataPtr == 0) || ( nbRows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == nbRows) && nbCols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == nbCols))",
    file=0x41f3d8 "/usr/include/eigen3/Eigen/src/Core/MapBase.h",
    line=148,
    function=0x427000 <Eigen::MapBase<Eigen::Block<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false>, 0>::MapBase(double const*, long, long)::__PRETTY_FUNCTION__> "Eigen::MapBase<Derived, 0>::MapBase(Eigen::MapBase<Derived, 0>::PointerType, Eigen::MapBase<Derived, 0>::Index, Eigen::MapBase<Derived, 0>::Index) [with Derived = Eigen::Block<const Eigen::Matrix<doub"...) at assert.c:101
#4  0x000000000041a0b6 in Eigen::MapBase<Eigen::Block<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false>, 0>::MapBase (this=0x7fffffffd190,
    dataPtr=0x7fffffffd498, nbRows=-1, nbCols=1)
    at /usr/include/eigen3/Eigen/src/Core/MapBase.h:146
#5  0x0000000000418d40 in Eigen::internal::BlockImpl_dense<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false, true>::BlockImpl_dense (
    this=0x7fffffffd190, xpr=..., startRow=3, startCol=2, blockRows=-1,
    blockCols=1) at /usr/include/eigen3/Eigen/src/Core/Block.h:350
#6  0x00000000004176d8 in Eigen::BlockImpl<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false, Eigen::Dense>::BlockImpl (this=0x7fffffffd190,
    xpr=..., a_startRow=3, a_startCol=2, blockRows=-1, blockCols=1)
    at /usr/include/eigen3/Eigen/src/Core/Block.h:159
#7  0x0000000000416036 in Eigen::Block<Eigen::Matrix<double, 2, 2, 0, 2, 2> const, -1, 1, false>::Block (this=0x7fffffffd190, xpr=..., a_startRow=3,
    a_startCol=2, blockRows=-1, blockCols=1)
    at /usr/include/eigen3/Eigen/src/Core/Block.h:136
#8  0x0000000000412452 in Eigen::internal::hseq_side_dependent_impl<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, 1>::essentialVector (h=..., k=2)
    at /usr/include/eigen3/Eigen/src/Householder/HouseholderSequence.h:85
#9  0x000000000041068b in Eigen::HouseholderSequence<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, 1>::essentialVector
    (this=0x7fffffffd2a0, k=2)
    at /usr/include/eigen3/Eigen/src/Householder/HouseholderSequence.h:199
#10 0x000000000040e4f5 in Eigen::HouseholderSequence<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, 1>::applyThisOnTheLeft<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 1, 1, 1, 1, 1> > (this=0x7fffffffd2a0, dst=..., workspace=...)
    at /usr/include/eigen3/Eigen/src/Householder/HouseholderSequence.h:314
#11 0x000000000040c4b3 in Eigen::HouseholderSequence<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, 1>::applyThisOnTheLeft<Eigen::Matrix<double, 2, 1, 0, 2, 1> > (this=0x7fffffffd2a0, dst=...)
    at /usr/include/eigen3/Eigen/src/Householder/HouseholderSequence.h:303
#12 0x000000000040a941 in Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >::applyOnTheLeft<Eigen::HouseholderSequence<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, 1> > (
    this=0x7fffffffd290, other=...)
    at /usr/include/eigen3/Eigen/src/Core/EigenBase.h:156
#13 0x0000000000408a13 in Eigen::internal::solve_retval<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> >::evalTo<Eigen::Matrix<double, 2, 1, 0, 2, 1> > (
    this=0x7fffffffd5d0, dst=...)
    at /usr/include/eigen3/Eigen/src/QR/ColPivHouseholderQR.h:540
#14 0x0000000000406d6f in Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> >::evalTo<Eigen::Matrix<double, 2, 1, 0, 2, 1> > (
    this=0x7fffffffd5d0, dst=...)
    at /usr/include/eigen3/Eigen/src/misc/Solve.h:51
#15 0x00000000004053c3 in Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> > >::evalTo<Eigen::Matrix<double, 2, 1, 0, 2, 1> > (this=0x7fffffffd5d0, dst=...)
    at /usr/include/eigen3/Eigen/src/Core/ReturnByValue.h:61
#16 0x0000000000403ed9 in Eigen::internal::assign_selector<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> >, false, false>::evalTo<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::ReturnByValue<Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> > > > (dst=..., other=...)
    at /usr/include/eigen3/Eigen/src/Core/Assign.h:522
#17 0x0000000000402f65 in Eigen::MatrixBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >::operator=<Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> > > (this=0x7fffffffd530, other=...)
    at /usr/include/eigen3/Eigen/src/Core/Assign.h:578
#18 0x000000000040266e in Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >::operator=<Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> > > (this=0x7fffffffd530, func=...)
    at /usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:418
#19 0x0000000000402023 in Eigen::Matrix<double, 2, 1, 0, 2, 1>::operator=<Eigen::internal::solve_retval_base<Eigen::ColPivHouseholderQR<Eigen::Matrix<double, 2, 2, 0, 2, 2> >, Eigen::Matrix<double, 2, 1, 0, 2, 1> > > (
    this=0x7fffffffd530, func=...)
    at /usr/include/eigen3/Eigen/src/Core/Matrix.h:190
#20 0x00000000004013d6 in main () at fails.cpp:20

This is because the ColPivHouseholderQR object is deleted at the end of the function and therefore the object returned by the solve() method stores dead references. In other words, this is because A.colPivHouseholderQr() does not return an abstract expression, but an object containing the result of the decomposition.

@ggael Thanks for your response.

Okay, I suspected something like this. But then, why does the double transpose fix this? Does it internally do an `eval`, or is it something else?

In other words, what is the correct way to implement my function `solve`? And is it possible to return an expression to avoid an extra copy?

The first transpose evaluate the result of the solve operation which is put on the stack. The cleanest way would be to either make the result an argument of your function: solve(A,b,x) (C-style) or make your solve function returns a pseudo expression inheriting EigenBase, storing references to A and b, and implementing an evalTo(x) member function doing the actual factorisation and solving.

@ggael Thank you for these suggestions.

Unfortunately, for technical reasons, the C-type signature is not really an option in my case.

So, I gave the pseudo expression a shot. Is there some documentation available on how to do this? I've read this page about customizing/extending Eigen, but couldn't find any instructions on how to write your own expressions. The documentation on expression templates also seems rather sparse.

Here is what I've got so far:

Code: Select all

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

namespace Eigen {

template <class Decomposition, class DerivedB> class SolveOp;

namespace internal {
template <class DerivedA, class DerivedB>
struct traits<SolveOp<
    DerivedA, DerivedB>> : public solve_retval<ColPivHouseholderQR<DerivedA>,
                                               DerivedB> {};
}

template <class DerivedA, class DerivedB>
class SolveOp : public EigenBase<SolveOp<DerivedA, DerivedB>> {
    DerivedA A_;
    DerivedB B_;

    typedef ColPivHouseholderQR<DerivedA> Decomposition;
    typedef internal::solve_retval<Decomposition, DerivedB> Sol;

  public:
    enum {
        RowsAtCompileTime = Sol::RowsAtCompileTime,
        ColsAtCompileTime = Sol::ColsAtCompileTime,
        MaxRowsAtCompileTime = Sol::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = Sol::MaxColsAtCompileTime,
        Flags = Sol::Flags
    };

  public:
    SolveOp() = default;
    SolveOp(const DerivedA &A, const DerivedB &B) : A_(A), B_(B) {}

    template <class Dest>
    void evalTo(Dest &dest) const {
        dest = A_.colPivHouseholderQr().solve(B_);
    }
};

}

using namespace Eigen;

template <class DerivedA, class DerivedB>
SolveOp<DerivedA, DerivedB>
solve(const MatrixBase<DerivedA> &A, const MatrixBase<DerivedB> &B) {
    return SolveOp<DerivedA, DerivedB>(A, B);
}

int main() {
    Matrix2d A;
    Vector2d b, c;

    A << 1, 2,
         4, 3;
    b << 5, 10;

    c = solve(A, b);

    std::cout << c << "\n";
}

It compiles, but when run, it fails with a segmentation fault. According to GDB that happens in this line:

Code: Select all

SolveOp(const DerivedA &A, const DerivedB &B) : A_(A), B_(B) {}

And the stack-trace consists of a very, very, very long list of frames like the following

Code: Select all

#7  0x000000000042309d in Eigen::EigenBase<Eigen::SolveOp<Eigen::Matrix<double, 2, 2, 0, 2, 2>, Eigen::Matrix<double, 2, 1, 0, 2, 1> > >::rows (
    this=0x7fffffffd4f0)
    at /usr/include/eigen3/Eigen/src/Core/EigenBase.h:44

.

What am I doing wrong?

Unfortunately, I still haven't been able to solve this problem.

May I ask for some pointers on how to properly write my own expressions in Eigen?

ggael wrote:This is because the ColPivHouseholderQR object is deleted at the end of the function and therefore the object returned by the solve() method stores dead references. In other words, this is because A.colPivHouseholderQr() does not return an abstract expression, but an object containing the result of the decomposition.

ggael wrote:The first transpose evaluate the result of the solve operation which is put on the stack. The cleanest way would be to either make the result an argument of your function: solve(A,b,x) (C-style) or make your solve function returns a pseudo expression inheriting EigenBase, storing references to A and b, and implementing an evalTo(x) member function doing the actual factorisation and solving.

Thanks a lot!! Actually you pointed the miss things for me.

ggael wrote:This is because the ColPivHouseholderQR object is deleted at the end of the function and therefore the object returned by the solve() method stores dead references. In other words, this is because A.colPivHouseholderQr() does not return an abstract expression, but an object containing the result of the decomposition.

It makes sense. The reason for this issue could be because the A.colPivHouseholderQr() method returns an object that contains the result of a matrix decomposition operation, rather than an abstract expression.