Eigen学习（三） 矩阵和向量的运算

继续翻译，原文链接：http://eigen.tuxfamily.org/dox/group__TutorialArrayClass.html

这一节主要介绍如何在Eigen中实现矩阵、向量及标量之间的运算。Eigen提供了一些矩阵和向量的数值运算，其中一些是通过通用的C++运算符重载实现，如+,-,*等，另一些通过特殊的方法实现，如dot()，cross()等方法。对于Matrix类，这些操作只支持线性代数运算，比如matrix1*matrix2代表矩阵的乘积，而向量+标量则是不允许的。

加和减

运算符两侧必须有相同的行和列。而且还要有相同的数据类型，因为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];

因此你并不需要担心在Eigen中使用大的算法表达式。它仅仅是给了Eigen更多的优化机会。

转置和共轭

对矩阵的转置、共轭和共轭转置由成员函数transpose(),conjugate(),adjoint()实现

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;
cout << "Here is the matrix a^T\n" << a.transpose() << endl;
cout << "Here is the conjugate of a\n" << 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()等同于reanspose()。

至于基本的数值运算符、transpose()和adjoint()会简单的返回一个中间对象而不是对原对象做处理。比如你做b=a.transpose()，那么a会保存不变。然而当你做a=a.transpose()时Eigen会在转置执行结束前就往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没有禁止时，这种问题会被自动检测到。要避免错误，可以使用in-place转置。类似的还有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，等于a = 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

注：Eigen对矩阵相乘这种情况进行特殊处理，不用担心别名问题，比如m = m*m会被编译成tmp = m*m; m = tmp。

点积和叉积

对于点积和叉积，直接使用dot()和cross()方法

#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

记住叉积仅仅用于尺寸为3的向量！点积可以用于任意尺寸的向量，当使用复数时，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 << ")\n\n";
  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会检查你的操作的有效性，如果可能的话会在编译阶段产生，输出信息。这些信息可能比较长并且难看，不过重要的信息都用了大写的语句标出了，可以很容易找到。比如

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

但是另一些情况下，比如对于动态尺寸的检查不能在编译阶段完成，那么Eigen会在运行期间检查，如果出错，那么程序会崩溃并输出一些信息。

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