Eigen教程(3)

整理下Eigen库的教程,参考:http://eigen.tuxfamily.org/dox/index.html

矩阵和向量的运算

提供一些概述和细节:关于矩阵、向量以及标量的运算。

介绍

Eigen提供了matrix/vector的运算操作,既包括重载了c++的算术运算符+/-/*,也引入了一些特殊的运算比如点乘dot、叉乘cross等。

对于Matrix类(matrix和vectors)这些操作只支持线性代数运算,比如:matrix1*matrix2表示矩阵的乘机,vetor+scalar是不允许的。如果你想执行非线性代数操作,请看下一篇(暂时放下)。

加减

左右两侧变量具有相同的尺寸(行和列),并且元素类型相同(Eigen不自动转化类型)操作包括:

  • 二元运算 + 如a+b
  • 二元运算 - 如a-b
  • 一元运算 - 如-a
  • 复合运算 += 如a+=b
  • 复合运算 -= 如a-=b
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  Matrix2d a;
  a << 1, 2,
       3, 4;
  MatrixXd b(2,2);
  b << 2, 3,
       1, 4;
  std::cout << "a + b =n" << a + b << std::endl;
  std::cout << "a - b =n" << a - b << std::endl;
  std::cout << "Doing a += b;" << std::endl;
  a += b;
  std::cout << "Now a =n" << a << std::endl;
  Vector3d v(1,2,3);
  Vector3d w(1,0,0);
  std::cout << "-v + w - v =n" << -v + w - v << std::endl;
}

输出:

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

标量乘法和除法

乘/除标量是非常简单的,如下:

  • 二元运算 * 如matrix*scalar
  • 二元运算 * 如scalar*matrix
  • 二元运算 / 如matrix/scalar
  • 复合运算 *= 如matrix*=scalar
  • 复合运算 /= 如matrix/=scalar
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  Matrix2d a;
  a << 1, 2,
       3, 4;
  Vector3d v(1,2,3);
  std::cout << "a * 2.5 =n" << a * 2.5 << std::endl;
  std::cout << "0.1 * v =n" << 0.1 * v << std::endl;
  std::cout << "Doing v *= 2;" << std::endl;
  v *= 2;
  std::cout << "Now v =n" << v << std::endl;
}

结果

a * 2.5 =
2.5   5
7.5  10
0.1 * v =
0.1
0.2
0.3
Doing v *= 2;
Now v =
2
4
6

表达式模板

这里简单介绍,在高级主题中会详细解释。在Eigen中,线性运算比如+不会对变量自身做任何操作,会返回一个“表达式对象”来描述被执行的计算。当整个表达式被评估完(一般是遇到=号),实际的操作才执行。

这样做主要是为了优化,比如

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;

Eigen会编译这段代码最终遍历一次即可运算完成。

for(int i = 0; i < 50; ++i)
  a[i] = 3*b[i] + 4*c[i] + 5*d[i];

因此,我们不必要担心大的线性表达式的运算效率。

转置和共轭

img 表示transpose转置

img 表示conjugate共轭

img 表示adjoint(共轭转置) 伴随矩阵

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix an" << a << endl;
cout << "Here is the matrix a^Tn" << a.transpose() << endl;
cout << "Here is the conjugate of an" << a.conjugate() << endl;
cout << "Here is the matrix a^*n" << a.adjoint() << endl;

输出

Here is the matrix a
 (-0.211,0.68) (-0.605,0.823)
 (0.597,0.566)  (0.536,-0.33)
Here is the matrix a^T
 (-0.211,0.68)  (0.597,0.566)
(-0.605,0.823)  (0.536,-0.33)
Here is the conjugate of a
 (-0.211,-0.68) (-0.605,-0.823)
 (0.597,-0.566)    (0.536,0.33)
Here is the matrix a^*
 (-0.211,-0.68)  (0.597,-0.566)
(-0.605,-0.823)    (0.536,0.33)

对于实数矩阵,conjugate不执行任何操作,adjoint等价于transpose。

transpose和adjoint会简单的返回一个代理对象并不对本省做转置。如果执行 b=a.transpose() ,a不变,转置结果被赋值给b。如果执行 a=a.transpose() Eigen在转置结束之前结果会开始写入a,所以a的最终结果不一定等于a的转置。

Matrix2i a; a << 1, 2, 3, 4;
cout << "Here is the matrix a:n" << a << endl;
a = a.transpose(); // !!! do NOT do this !!!
cout << "and the result of the aliasing effect:n" << a << endl;

Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4

这被称为“别名问题”。在debug模式,当assertions打开的情况加,这种常见陷阱可以被自动检测到。

a=a.transpose() 这种操作,可以执行in-palce转置。类似还有adjointInPlace。

MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
cout << "Here is the initial matrix a:n" << a << endl;
a.transposeInPlace();
cout << "and after being transposed:n" << a << endl;

Here is the initial matrix a:
1 2 3
4 5 6
and after being transposed:
1 4
2 5
3 6

矩阵-矩阵的乘法和矩阵-向量的乘法

向量也是一种矩阵,实质都是矩阵-矩阵的乘法。

  • 二元运算 *如a*b
  • 复合运算 *=如a*=b
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
  Matrix2d mat;
  mat << 1, 2,
         3, 4;
  Vector2d u(-1,1), v(2,0);
  std::cout << "Here is mat*mat:n" << mat*mat << std::endl;
  std::cout << "Here is mat*u:n" << mat*u << std::endl;
  std::cout << "Here is u^T*mat:n" << u.transpose()*mat << std::endl;
  std::cout << "Here is u^T*v:n" << u.transpose()*v << std::endl;
  std::cout << "Here is u*v^T:n" << u*v.transpose() << std::endl;
  std::cout << "Let's multiply mat by itself" << std::endl;
  mat = mat*mat;
  std::cout << "Now mat is mat:n" << mat << std::endl;
}

输出

Here is mat*mat:
 7 10
15 22
Here is mat*u:
1
1
Here is u^T*mat:
2 2
Here is u^T*v:
-2
Here is u*v^T:
-2 -0
 2  0
Let's multiply mat by itself
Now mat is mat:
 7 10
15 22

m=m*m并不会导致别名问题,Eigen在这里做了特殊处理,引入了临时变量。实质将编译为:

tmp = m*m
m = tmp

如果你确定矩阵乘法是安全的(并没有别名问题),你可以使用noalias()函数来避免临时变量 c.noalias() += a*b

点运算和叉运算

dot()执行点积,cross()执行叉积,点运算得到1*1的矩阵。当然,点运算也可以用u.adjoint()*v来代替。

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
  Vector3d v(1,2,3);
  Vector3d w(0,1,2);
  cout << "Dot product: " << v.dot(w) << endl;
  double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
  cout << "Dot product via a matrix product: " << dp << endl;
  cout << "Cross product:n" << v.cross(w) << endl;
}

输出

Dot product: 8
Dot product via a matrix product: 8
Cross product:
 1
-2
 1

注意:点积只对三维vector有效。对于复数,Eigen的点积是第一个变量共轭和第二个变量的线性积。

基础的归约操作

Eigen提供了而一些归约函数:sum()、prod()、maxCoeff()和minCoeff(),他们对所有元素进行操作。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
  Eigen::Matrix2d mat;
  mat << 1, 2,
         3, 4;
  cout << "Here is mat.sum():       " << mat.sum()       << endl;
  cout << "Here is mat.prod():      " << mat.prod()      << endl;
  cout << "Here is mat.mean():      " << mat.mean()      << endl;
  cout << "Here is mat.minCoeff():  " << mat.minCoeff()  << endl;
  cout << "Here is mat.maxCoeff():  " << mat.maxCoeff()  << endl;
  cout << "Here is mat.trace():     " << mat.trace()     << endl;
}

输出

Here is mat.sum():       10
Here is mat.prod():      24
Here is mat.mean():      2.5
Here is mat.minCoeff():  1
Here is mat.maxCoeff():  4
Here is mat.trace():     5

trace表示矩阵的迹,对角元素的和等价于 a.diagonal().sum()

minCoeff和maxCoeff函数也可以返回结果元素的位置信息。

Matrix3f m = Matrix3f::Random();
  std::ptrdiff_t i, j;
  float minOfM = m.minCoeff(&i,&j);
  cout << "Here is the matrix m:n" << m << endl;
  cout << "Its minimum coefficient (" << minOfM 
       << ") is at position (" << i << "," << j << ")nn";
  RowVector4i v = RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  cout << "Here is the vector v: " << v << endl;
  cout << "Its maximum coefficient (" << maxOfV 
       << ") is at position " << i << endl;

输出

Here is the matrix m:
  0.68  0.597  -0.33
-0.211  0.823  0.536
 0.566 -0.605 -0.444
Its minimum coefficient (-0.605) is at position (2,1)

Here is the vector v:  1  0  3 -3
Its maximum coefficient (3) is at position 2

操作的有效性

Eigen会检测执行操作的有效性,在编译阶段Eigen会检测它们,错误信息是繁冗的,但错误信息会大写字母突出,比如:

Matrix3f m;
Vector4f v;
v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

当然动态尺寸的错误要在运行时发现,如果在debug模式,assertions会触发后,程序将崩溃。

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"

原文链接: https://www.cnblogs.com/houkai/p/6348044.html

欢迎关注

微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍

    Eigen教程(3)

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/248420

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年2月14日 上午2:53
下一篇 2023年2月14日 上午2:54

相关推荐