Eigen学习笔记(3)-矩阵和向量的运算

原文:Eigen官网-Matrix and vector arithmetic
本节内容主要介绍Eigen中关于矩阵、向量、标量之间的数学运算。

1. 引言

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

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

2. 加法和减法

左右两侧的变量都必须具有相同的大小(行和列)和数据类型(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

3. 标量乘法和除法

矩阵或者向量与标量进行乘法或除法操作是比较简单的。操作符如下:

  • 二元运算 * 如matrix*scalar
  • 二元运算 * 如scalar*matrix
  • 二元运算 / 如matrix/scalar
  • 复合运算 *= 如matrix*=scalar
  • 复合运算 /= 如matrix/=scalar

4. 转置和共轭

a

T

a^T

aT表示转置(transpose);

a

ˉ

\\bar{a}

aˉ表示共轭(conjugate);

a

a^*

a表示共轭转置(伴随矩阵)(adjoint)。
上述三个操作分别通过transpose(), conjugate(),以及adjoint()函数实现。
示例如下:

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

using namespace Eigen;

int main()
{
Matrix2d a;
a << 1, 2,
     3, 4;
std::cout << "Here is the matrix a\\n" << a << std::endl;
std::cout << "Here is the transpose of a\\n" << a.transpose() << std::endl;
std::cout << "Here is the conjugate of a\\n" << a.conjugate() << std::endl;
std::cout << "Here is the adjoint of a\\n" << a.adjoint() << std::endl;
return 0;
}

结果如下:

Here is the matrix a
1 2
3 4
Here is the transpose of a
1 3
2 4
Here is the conjugate of a
1 2
3 4
Here is the adjoint of a
1 3
2 4

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

transpose和adjoint会简单的返回一个代理对象并不做真正的转置。如果执行 b=a.transpose() ,a不变,转置结果被赋值给b。如果执行 a=a.transpose() Eigen在转置结束之前结果会开始写入a,所以a的最终结果不一定等于a的转置。这就是所谓的别名问题。在“调试模式”中,也就是说,当断言没有被禁用时,这种常见的陷阱会被自动检测到。

对于a = a.transpose()这种操作,可以使用transposeInPlace()解决,类似的还有adjointInPlace()
示例如下:

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

using namespace Eigen;

int main()
{
Matrix2d a;
a << 1, 2,
     3, 4;
std::cout << "Here is the matrix a\\n" << a << std::endl;

//a = a.transpose(); //!!!do Not do this!!!在这里无法编译通过
//std::cout << "the result of the aliasing effect: \\n" << a << std::endl;

Matrix2d b;
b << 1, 2,
     3, 4;
b.transposeInPlace();
std::cout << " b after being transposed:\\n" << b << std::endl;

return 0;
}

结果如下:

Here is the matrix a
1 2
3 4
 b after being transposed:
1 3
2 4

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

因为向量是一种特殊的矩阵,因此本质上都是矩阵-矩阵的乘法。运算符如下:

  • 二元运算* 如a*b
  • 复合运算 *= 如a*=b (这种右乘操作,a*=b等价于a=a*b)

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

tmp = m*m
m = tmp

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

c.noalias() += a*b

6. 点乘和叉乘

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的点积是第一个变量共轭和第二个变量的线性积。

7. 基础的归约操作

Eigen提供了一些归约函数把一个给定的矩阵或向量化为一个值,比如

  • 求和(sum());
  • 累积乘积(prod()),矩阵中所有元素的乘积;
  • 平均值(mean()),矩阵中所有元素的平均值;
  • 最大值(maxCoeff()),矩阵中所有元素的最大值,也可以返回最大值的位置;
  • 最小值(minCoeff()),矩阵中所有元素的最小值,也可以返回最小值的位置;
  • 迹(trace()),一个矩阵中对角线元素的和,也可以用a.diagonal().sum()进行计算。
    示例如下:
#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;

  ptrdiff_t i,j;
  double minOfMat = mat.minCoeff(&i, &j);
  cout << "Its minimum coefficient (" << minOfMat << ") is at position (" << i << "," << j << ")\\n\\n";

  Eigen::RowVector4i v = Eigen::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;

  return 0;
}

结果如下:

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
Its minimum coefficient (1) is at position (0,0)

Here is the vector v:  730547559 -226810938  607950953  640895091
Its maximum coefficient (730547559) is at position 0

8. 操作的有效性

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"

参考:

“Eigen教程(3)”

© 版权声明
THE END
喜欢就支持一下吧
点赞590 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片

    暂无评论内容