Eigen 中 maxCoeff 的奇怪错误
Strange bug with maxCoeff in Eigen
我有:
VectorXd my_vector;
MatrixXd my_matrix;
和维度在下面的表达式中是兼容的。
我记得在 Eigen 文档中读到,作为一个整体评估表达式总是比评估较小的块然后重新组合结果更好,因为它为 Eigen 提供了更多优化机会。
考虑到这一点,我会使用:
// (1)
int index;
(my_vector.transpose() * my_matrix.leftCols(n_cols+1)).maxCoeff(&index);
而不是:
// (2)
int index;
RowVectorXd temp = my_vector.transpose() * my_matrix.leftCols(n_cols+1);
temp.maxCoeff(&index);
然而,当 n_cols
= 1 时 (1) 崩溃,而在同一个地方使用 (2) 而不是 (1) 似乎总是有效。 (注意 my_matrix
有超过 1 列...)
问题:
为什么 (1) 崩溃而 (2) 不崩溃???
如何防止 (1) 在将其保持为单个表达式时崩溃?
错误来自 Eigen 的内部机制断言,在我看来,这是因为在 (1) 中 maxCoeff
似乎没有意识到它是在行向量而不是一般矩阵上调用的- 如果我错了请纠正我。
作为解决方法,您可以使用 .eval()
。这实际上与 (2) 相同,只是没有额外的代码行:
(my_vector.transpose() * my_matrix.leftCols(n_cols + 1)).eval().maxCoeff(&index);
关于 Eigen 的内部工作原理,当您在情况 (1) 中调用 maxCoeff
时,被调用的对象是(正如您明确希望的那样)惰性求值结构(在本例中 GeneralProduct<Transpose<..>,Block<...> >
) 由调用 mat.coeff(0,0)
触发,其中 mat
是 GeneralProduct
。在内部,mat
有 (n_cols + 1)
列,而不是断言后的行所需的 1 (Matrix<Scalar,1,1> result = *this;
)
这就是导致评估的原因,因此我们每次只想评估一个元素。如果没有断言,尝试将 *this
分配给 result
会触发调整大小,导致程序崩溃。
当使用临时对象时,对象是 Matrix<double,1,-1,1,1,-1>
所以 coeff(i,j)
只是对 m_storage.data()[colId + rowId * m_storage.cols()]
.
的简单调用
所有这些都成为 Eigen 3.3 中的一个有争议的点,因为整个 coeff(i,j)
现在是 return m_evaluator.coeff(row, col);
TLDR
升级到 Eigen 3.3(alpha),问题就消失了。
为了完成 Avi 的回答,让我们首先补充一点,以便有效地评估,无论如何,产品都必须评估为临时的,所以无论 Eigen 的版本如何,调用 .eval()
都完全没问题。
您正在观察的当前断言:
Assertion failed: (this->rows() == 1 && this->cols() == 1), function coeff, file ../eigen3.2/Eigen/src/Core/ProductBase.h, line 148
是为了确保乘积表达式不会意外地按系数求值(内积除外)。可以在断言消息中添加此解释以使其更清楚。
然而,如果您没有明确调用 coeff()
或 operator(i,j)
自己,因此这是 Eigen 中的错误,更准确地说是 Eigen::Visitor
中的错误,应该评估为您嵌套产品表达式。这将在 3.2.8 中修复。
最后,如果您知道 n_cols+1
和 my_vector
都非常小,您可以使用惰性产品来避免这种临时情况:
my_vector.transpose().lazyProduct(my_matrix.leftCols(n_cols+1)).maxCoeff(&index);
编辑: 这里记录的是 fix.
我有:
VectorXd my_vector;
MatrixXd my_matrix;
和维度在下面的表达式中是兼容的。
我记得在 Eigen 文档中读到,作为一个整体评估表达式总是比评估较小的块然后重新组合结果更好,因为它为 Eigen 提供了更多优化机会。
考虑到这一点,我会使用:
// (1)
int index;
(my_vector.transpose() * my_matrix.leftCols(n_cols+1)).maxCoeff(&index);
而不是:
// (2)
int index;
RowVectorXd temp = my_vector.transpose() * my_matrix.leftCols(n_cols+1);
temp.maxCoeff(&index);
然而,当 n_cols
= 1 时 (1) 崩溃,而在同一个地方使用 (2) 而不是 (1) 似乎总是有效。 (注意 my_matrix
有超过 1 列...)
问题:
为什么 (1) 崩溃而 (2) 不崩溃???
如何防止 (1) 在将其保持为单个表达式时崩溃?
错误来自 Eigen 的内部机制断言,在我看来,这是因为在 (1) 中 maxCoeff
似乎没有意识到它是在行向量而不是一般矩阵上调用的- 如果我错了请纠正我。
作为解决方法,您可以使用 .eval()
。这实际上与 (2) 相同,只是没有额外的代码行:
(my_vector.transpose() * my_matrix.leftCols(n_cols + 1)).eval().maxCoeff(&index);
关于 Eigen 的内部工作原理,当您在情况 (1) 中调用 maxCoeff
时,被调用的对象是(正如您明确希望的那样)惰性求值结构(在本例中 GeneralProduct<Transpose<..>,Block<...> >
) 由调用 mat.coeff(0,0)
触发,其中 mat
是 GeneralProduct
。在内部,mat
有 (n_cols + 1)
列,而不是断言后的行所需的 1 (Matrix<Scalar,1,1> result = *this;
)
这就是导致评估的原因,因此我们每次只想评估一个元素。如果没有断言,尝试将 *this
分配给 result
会触发调整大小,导致程序崩溃。
当使用临时对象时,对象是 Matrix<double,1,-1,1,1,-1>
所以 coeff(i,j)
只是对 m_storage.data()[colId + rowId * m_storage.cols()]
.
所有这些都成为 Eigen 3.3 中的一个有争议的点,因为整个 coeff(i,j)
现在是 return m_evaluator.coeff(row, col);
TLDR
升级到 Eigen 3.3(alpha),问题就消失了。
为了完成 Avi 的回答,让我们首先补充一点,以便有效地评估,无论如何,产品都必须评估为临时的,所以无论 Eigen 的版本如何,调用 .eval()
都完全没问题。
您正在观察的当前断言:
Assertion failed: (this->rows() == 1 && this->cols() == 1), function coeff, file ../eigen3.2/Eigen/src/Core/ProductBase.h, line 148
是为了确保乘积表达式不会意外地按系数求值(内积除外)。可以在断言消息中添加此解释以使其更清楚。
然而,如果您没有明确调用 coeff()
或 operator(i,j)
自己,因此这是 Eigen 中的错误,更准确地说是 Eigen::Visitor
中的错误,应该评估为您嵌套产品表达式。这将在 3.2.8 中修复。
最后,如果您知道 n_cols+1
和 my_vector
都非常小,您可以使用惰性产品来避免这种临时情况:
my_vector.transpose().lazyProduct(my_matrix.leftCols(n_cols+1)).maxCoeff(&index);
编辑: 这里记录的是 fix.