Eigen库类类型作为函数参数

最新推荐文章于 2024-05-30 12:30:34 发布

-名字不重要-

最新推荐文章于 2024-05-30 12:30:34 发布

阅读量3.5k

点赞数 1

分类专栏： C/C++ 文章标签： Eigen C++ 数值计算库函数参数矩阵库

C/C++ 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

英文原文衔接：http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html

Eigen库类类型作为函数参数

由于Eigen使用了表达式模板（expression templates），潜在的导致了每一个表达式都可能是一个不同的类型。如果传递这样的表达式给一个矩阵(Martix type)类型参数的函数，表达式将被隐含地提升（evaluated）为一个临时矩阵（temporary Matrix），然后传递给函数。这意味着不能有效利用表达式模板抑制临时变量的好处。具体地说，这有两个缺点：

*转换为临时变量可能是无用且不高效的；

*只允许函数读取表达式，不能修改表达式。

幸运的是,所有的这些表达式的类型的都有一个共同点,他们都继承至一些共同的,模板化的基类。通过让函数的参数为这些基本类型的模板参数,可让函数很好地处理Eigen的表达式模板。

开头的一些例子

这部分将针对Eigen库提供的不同类型的对象，展示一些例子。在实际例子开始之前，我们需要概括哪些基本对象我们可以利用（也可查看类继承部分）。

MatrixBase: 所有稠密矩阵表达式的共有基类(不是数组表达式,也不是稀疏和特殊矩阵类)。在函数中使用它，意味着只能使用稠密矩阵。

ArrayBase: 所有密集数组表达式的共有基类(不是矩阵表达式,等等)。在函数中使用它，意味着只能使用数组。

DenseBase: 所有稠密矩阵、数组表达式的共有基类,即MatrixBase和ArrayBase的基类。它可以用于使用矩阵和数组的函数。

EigenBase: 所有类型的对象的统一基类，可以被提升为稠密矩阵或数组,例如特殊矩阵类,如对角矩阵（diagonal matrices）,置换矩阵（permutation matrices）等。它可用于使用任何这样的一般类型的函数。

EigenBase 例子

打印在Eigen库中最常用的对象的维数. 它可以是任何矩阵表达式,任何稠密或稀疏矩阵和数组。

Example:
#include <iostream>
#include <Eigen/Core>
using namespace Eigen;
template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
  std::cout << "size (rows, cols): " << b.size() << " (" << b.rows()
            << ", " << b.cols() << ")" << std::endl;
}
int main()
{
    Vector3f v;
    print_size(v);
    // v.asDiagonal() returns a 3x3 diagonal matrix pseudo-expression
    print_size(v.asDiagonal());
}

Output:
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是一个表达式。

这些例子只是为了给读者第一印象：如何写一个函数以普通的常数矩阵或数组为参数。他们也旨在让读者了解关于函数的最优候选的基类。在下一节中我们将更详细地看一个例子,它可以实现的不同的方式,在讨论每个实现的问题和优势。下面讨论,矩阵和数组以及MatrixBase ArrayBase仍可以交换和所有参数仍不变。

如何写一个通用而非模板的函数?

在前面的例子中,所有函数都被定义为模板函数。这种方法允许编写通用代码,但通常我们渴望写非模板函数，而且仍然保持一定程度的通用性（genericity）避免繁琐的参数复制（相当于使用引用）。典型的例子是写函数接受一个MatrixXf或MatrixXf的子区块。这正是Ref类的目的。

下面是一个简单的例子:

Example:

#include <iostream>
#include <Eigen/SVD>
using namespace Eigen;
using namespace std;
float inv_cond(const Ref<const MatrixXf>& a)
{
  const VectorXf sing_vals = a.jacobiSvd().singularValues();
  return sing_vals(sing_vals.size()-1) / sing_vals(0);
}
int main()
{
  Matrix4f m = Matrix4f::Random();
  cout << "matrix m:" << endl << m << endl << endl;
  cout << "inv_cond(m):          " << inv_cond(m)                      << endl;
  cout << "inv_cond(m(1:3,1:3)): " << inv_cond(m.topLeftCorner(3,3))   << endl;
  cout << "inv_cond(m+I):        " << inv_cond(m+Matrix4f::Identity()) << endl;
}

Output:

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.0562344
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类参考文档。

什么情况下函数能有效的接受纯 Matrix或Array作为实参数？

不使用模板函数，不使用Ref类，前一个cov函数可以简单的表示为：

MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
  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;
  return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

起初,你可能会认为这个实现很好,但事实上这恰恰相反，除非你需要一个通用的实现与双矩阵,除非你不关心临时对象。这是为什么?涉及到临时变量在哪里?下面给出的代码如何编译?

MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);

在这种特殊情况下,这个例子很好且有效，因为两个参数都声明为常量引用。编译器创建一个临时变量并且将表达式y + z赋值给这个临时变量。一旦函数处理,临时变量被释放,结果被赋给C。

注：接受常引用Matrix (or Array)的函数可以以临时变量的代价处理表达式。

什么情况下函数接受纯Matrix或者Array实参失败？

在这里,我们考虑上面给出的函数略微的修改版本。这一次,我们不返回结果而是通过传递非常量参数来储存结果。一个简单的实现可能看起来如下。

// Note: 此代码有缺陷!
void cov(const MatrixXf& x, const MatrixXf& y, 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;
}

当尝试执行下面的代码：

MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));

编译器将会报错，因为在返回值时，将MatrixXf::block() 转换为非常量MatrixXf&是不可行的。出现这种情况,因为编译器想保护你，避免将结果写入一个临时对象。在这种特殊情况下这种保护不是目的,我们想写一个临时对象。那么我们如何解决这个问题呢?

目前首选的解决方案是基于一个小技巧（ a little hack）。需要通过一个常量引用矩阵和内部常量性需要抛弃。C98兼容的编译器正确的实现如下：

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;
  const Scalar num_observations = static_cast<Scalar>(x.rows());
  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;
  const_cast< MatrixBase<OtherDerived>& >(C) =
    (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

上面的实现现在不仅可以使用临时表达式,它还允许使用任意浮点标量类型矩阵的函数。

注意:const_cast技巧只能使用模板函数。它不会与MatrixXf实现,因为它是不可能把一区块给一个矩阵参考!

怎样通用实现改变矩阵大小功能？

有人可能会认为我们已经做了,对吧? 这并不完全正确,因为为了我们的协方差函数一般适用,我们希望下面的代码工作

MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);

当我们使用MatrixBase作为参数,情况就不是这样了。一般来说,Eigen支持自动调整,但在表达式上这样做是不可能的。为什么要允许调整矩阵的块?这是一个引用子矩阵,我们绝对不想调整。如果我们不能调整MatrixBase，那么我们如何能完成调整? 解决方法是调整的派生对象实现。

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C_)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;
  const Scalar num_observations = static_cast<Scalar>(x.rows());
  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;
  MatrixBase<OtherDerived>& C = const_cast< MatrixBase<OtherDerived>& >(C_);
 
  C.derived().resize(x.cols(),x.cols()); // resize the derived object
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

当参数是表达式或者参数是矩阵，但矩阵大小不正确时，上面的实现也能正常工作。这意味着,传递一个错误尺寸的表达式会导致运行时错误(仅在调试模式),而传递正确大小的表达式就能工作的很好。

注意:在上面的讨论中， Matrix 、 Array、MatrixBase和 ArrayBase仍可以交换和所有参数保持。

总结

总而言之,函数的实现采取常量引用(non-writable)对象不是一个大问题,在编译和运行您的程序时不会导致问题。然而,一个天真的实现是可能在代码中引入不必要的临时对象。为了避免产生临时变量, 将它们作为(常量)引用传递给MatrixBase或ArrayBase(所以模板化你的函数)。

函数接受(non-const)实参数，必须采取常量引用，然后在函数体内抛弃常量性（constness）。

函数接收MatrixBase(或ArrayBase)对象的参数，并可能潜在性需要调整大小(如果它们是可调整大小的)，必须在衍生类(derived class)上调用resize()方法，返回时用derived()方法。

Eigen库还是比较好的一个C++数值计算库，在官网上可以看到其在许多项目中都有使用，使用时只包含头文件即可，使用可谓非常方便。

-名字不重要-

关注

1
点赞
踩
16

收藏

觉得还不错? 一键收藏
0
评论
Eigen库类类型作为函数参数

英文原文衔接：http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.htmlEigen库类类型作为函数参数由于Eigen使用了表达式模板（expression templates），潜在的导致了每一个表达式都可能是一个不同的类型。如果传递这样的表达式给一个矩阵(Martix type)类型参数的函数，表达式将被...
复制链接

扫一扫

专栏目录