英文原文(Writing Functions Taking Eigen Types as Parameters)
Eigen使用表达式模板的方式导致每个表达式的类型可能都不同。如果将这样的表达式传递给一个需要Matrix类型参数的函数,则表达式将隐式地被评估为一个临时Matrix,然后再传递给函数。这意味着失去了表达式模板的好处。具体而言,这有两个缺点:
- 将表达式评估为临时变量可能是无用且低效的;
- 这会限制函数只能读取表达式,无法对表达式进行修改。
幸运的是,所有这些表达式类型都有一些共同的基类和模板。通过让函数接受这些基类的模板参数,可以让它们与Eigen的表达式模板很好地协作。
一些开始的示例
本节将为 Eigen 提供的不同类型的对象提供简单的示例。在开始实际示例之前,我们需要概括一下可以使用哪些基础对象(另请参阅类层次结构)。
- MatrixBase:所有稠密矩阵表达式的公共基类(与数组表达式相对,与稀疏和特殊矩阵类相对)。在仅适用于稠密矩阵的函数中使用它。
- ArrayBase:所有稠密数组表达式(与矩阵表达式相对)的公共基类。在仅适用于数组的函数中使用它。
- DenseBase:所有稠密矩阵表达式的公共基类,即
MatrixBase
和ArrayBase
的基类。它可以用在同时适用于矩阵和数组的函数中。 - EigenBase:该基类统一了可以被计算为稠密矩阵或数组的所有对象类型,例如对角矩阵、置换矩阵等特殊的矩阵类。它可以用于处理任何此类通用类型的函数中。
EigenBase示例
打印 Eigen 中存在的最通用对象的尺寸。它可以是任何矩阵表达式、任何稠密或稀疏矩阵以及任何数组。
#include <iostream>
#include <Eigen/Core>
template <typename Derived>
void print_size(const Eigen::EigenBase<Derived>& b)
{
std::cout << "size (rows, cols): " << b.size() << " (" << b.rows()
<< ", " << b.cols() << ")" << std::endl;
}
int main()
{
Eigen::Vector3f v;
print_size(v);
// v.asDiagonal() returns a 3x3 diagonal matrix pseudo-expression
print_size(v.asDiagonal());
}
输出:
size (rows, cols): 3 (3, 1)
size (rows, cols): 9 (3, 3
DenseBase示例
打印稠密表达式的子块。接受任何稠密矩阵或数组表达式,但不接受稀疏对象,也不接受特殊矩阵类(例如 DiagonalMatrix
)。
template <typename Derived>
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
{
std::cout << "block: " << b.block(x,y,r,c) << std::endl;
}
ArrayBase示例
打印数组或数组表达式的最大系数。
template <typename Derived>
void print_max_coeff(const ArrayBase<Derived> &a)
{
std::cout << "max: " << a.maxCoeff() << std::endl;
}
MatrixBase示例
打印给定矩阵或矩阵表达式的逆条件数。
template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType& sing_vals = a.jacobiSvd().singularValues();
std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}
多个模板化参数示例
计算两点之间的欧几里德距离。
template <typename DerivedA,typename DerivedB>
typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
{
return (p1-p2).squaredNorm();
}
请注意,我们使用了两个模板参数,每个参数一个。这允许函数处理不同类型的输入,例如:
squaredist(v1,2*v2)
其中第一个参数 v1
是向量,第二个参数 2*v2
是表达式。
这些示例只是为了让读者对如何编写接受普通和常量Matrix或Array参数的函数有第一印象。它们还旨在为读者提供有关最常见基类是作为函数的最佳候选的想法。在下一节中,将更详细地说明一个示例以及可以实现它的不同方式,同时讨论每个实现的问题和优点。在下面的讨论中,Matrix
和Array
以及MatrixBase
和ArrayBase
可以互换使用,所有参数仍然保持不变。
如何编写通用但非模板化的函数?
在所有以前的示例中,函数都必须是模板函数。这种方法允许编写非常通用的代码,但通常希望编写非模板函数并仍然保持一定程度的通用性,以避免对参数进行愚蠢的复制。典型的例子是编写接受MatrixXf
或MatrixXf
块的函数。这正是Ref类的目的。这里是一个简单的例子:
#include <iostream>
#include <Eigen/SVD>
float inv_cond(const Eigen::Ref<const Eigen::MatrixXf>& a)
{
const Eigen::VectorXf sing_vals = a.jacobiSvd().singularValues();
return sing_vals(sing_vals.size()-1) / sing_vals(0);
}
int main()
{
Eigen::MatrixXf m = Eigen::MatrixXf::Random(4, 4);
std::cout << "matrix m:\n" << m << "\n\n";
std::cout << "inv_cond(m): " << inv_cond(m) << "\n";
std::cout << "inv_cond(m(1:3,1:3)): " << inv_cond(m.topLeftCorner(3,3)) << "\n";
std::cout << "inv_cond(m+I): " << inv_cond(m+Eigen::MatrixXf::Identity(4, 4)) << "\n";
}
输出:
matrix m:
0.68 0.823 -0.444 -0.27
-0.211 -0.605 0.108 0.0268
0.566 -0.33 -0.0452 0.904
0.597 0.536 0.258 0.832
inv_cond(m): 0.0562343
inv_cond(m(1:3,1:3)): 0.0836819
inv_cond(m+I): 0.160204
在对inv_cond
的前两个调用中,不会发生复制,因为参数的内存布局与Ref <MatrixXf>
接受的内存布局匹配。然而,在最后一次调用中,我们有一个通用表达式,将通过Ref <>
对象自动评估为一个临时的MatrixXf
。
Ref
对象也可以是可写的。这是一个计算两个输入矩阵的协方差矩阵的函数示例,其中每行都是一个观测值:
void cov(const Ref<const MatrixXf> x, const Ref<const MatrixXf> y, Ref<MatrixXf> C)
{
const float num_observations = static_cast<float>(x.rows());
const RowVectorXf x_mean = x.colwise().sum() / num_observations;
const RowVectorXf y_mean = y.colwise().sum() / num_observations;
C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}
下面是两个不带任何副本调用 cov
的示例:
MatrixXf m1, m2, m3
cov(m1, m2, m3);
cov(m1.leftCols<3>(), m2.leftCols<3>(), m3.topLeftCorner<3,3>());
Ref<>
类还有另外两个可选模板参数,允许控制无需任何副本即可接受的内存布局类型。有关详细信息,请参阅类 Ref 文档。