2. 矩阵和向量的基本数值计算
本章介绍对Matrix类型的矩阵、向量的数值计算方法。
2.1 加减运算符和标量乘除运算符
根据矩阵加减法的运算规则，我们知道参与运算的两个矩阵应当具有相同的行数和列数。 此外，由于C++的语法要求，两个Matrix对象的元素应当具有相同的数据类型。Eigen提供了如下的运算符：
- '+', 二元运算符，a + b
- '-', 二元运算符，a - b
- '-', 单元运算符，-a
- '+=' 复合运算符，a += b
- '-=' 复合运算符，a -= b
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
MatrixXd b(2,2);
a << 1, 2,
3, 4;
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中像'+'这样的算术运算符并不进行计算， 它只是返回一个“表达式对象(expression object)”来记录需要进行的计算。 当解析完整个表达式，执行赋值操作(通常是'='运算符)时再具体计算整个表达式的结果。 这种机制称为'expression templates'看起来很麻烦，但现代的编译器可以做一些优化。 例如，当我们需要计算如下的表达式时：
VectorXf a(50), b(50), c(50), d(50);
...
a = 3 * b + 4 * c + 5 * d;
实际上可以只用一个循环就完成所有的操作了：
for (int i = 0; i < 50; i++)
a[i] = 3*b[i] + 4*c[i] + 5*d[i];
因此，我们可以把一个表达式写的很长，这样Eigen可以由更多的优化可能。这种机制我们还会在后续的章节中予以详细的介绍。
2.2 转置和共轭运算
矩阵\(a\)的转置\(a^T\)，共轭\(\bar{a}\)，以及共轭转置\(a^*\) 可以分别通过成员函数'transpose()', 'conjugate()'和'adjoint()'获得。
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
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;
cout << "========" << endl;
MatrixXf b = MatrixXf::Random(2,2);
cout << "Here is the matrix b\n" << b << endl;
cout << "Here is the matrix b^T\n" << b.transpose() << endl;
cout << "Here is the conjugate of b\n" << b.conjugate() << endl;
cout << "Here is the matrix b^*\n" << b.adjoint() << endl;
}
其输出结果如下：
Here is the matrix a
(-0.211234,0.680375) (-0.604897,0.823295)
(0.59688,0.566198) (0.536459,-0.329554)
Here is the matrix a^T
(-0.211234,0.680375) (0.59688,0.566198)
(-0.604897,0.823295) (0.536459,-0.329554)
Here is the conjugate of a
(-0.211234,-0.680375) (-0.604897,-0.823295)
(0.59688,-0.566198) (0.536459,0.329554)
Here is the matrix a^*
(-0.211234,-0.680375) (0.59688,-0.566198)
(-0.604897,-0.823295) (0.536459,0.329554)
========
Here is the matrix b
-0.444451 -0.0452059
0.10794 0.257742
Here is the matrix b^T
-0.444451 0.10794
-0.0452059 0.257742
Here is the conjugate of b
-0.444451 -0.0452059
0.10794 0.257742
Here is the matrix b^*
-0.444451 0.10794
-0.0452059 0.257742
对于实数矩阵，共轭就是它本身所以conjugate()不做任何操作，而adjoint()与transpose()等价。 对于基本的算术运算符，transpose()和adjoint()只是返回一个代理对象并不做任何实际的转换。 只有在执行像 "b = a.transpose()"这样的操作时，才会立即计算a的转置并把结果赋予矩阵b。 但是当执行"a = a.transpose()"这样的操作时，Eigen就会在执行转置运算的同时写矩阵a的元素， 这样会导致a中的结果不是我们想要的a的转置，应当避免这样的操作。 这就是所谓的引用错误"aliasing issue"，在debug模式下，如果没有关闭断言"assertions"将会报错。 Eigen提供了'transposeInPlace()'函数来实现这类原地转换。类似的，也有adjointInPlace()
#define EIGEN_NO_DEBUG
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
Matrix2i a;
a << 1, 2, 3, 4;
cout << "Here is the initial matrix a:\n" << a << endl;
cout << "Here is the matrix a^T:\n" << a.transpose() << endl;
// 如下的操作将不能够得到a的转置
a = a.transpose();
cout << "After the Operation 'a = a.transpose()', the matrix a:\n" << a << endl;
a << 1, 2, 3, 4;
a.transposeInPlace();
cout << "After the Operation 'a.transposeInPlace()', the matrix a:\n" << a << endl;
}
其输出结果如下：
Here is the initial matrix a:
1 2
3 4
Here is the matrix a^T:
1 3
2 4
After the Operation 'a = a.transpose()', the matrix a:
1 2
2 4
After the Operation 'a.transposeInPlace()', the matrix a:
1 3
2 4
2.3 矩阵乘法运算
矩阵的乘法运算也是通过运算符*实现的。因为向量也是一种特殊的矩阵，所以这里讲的矩阵的乘法对向量同样适用。 Eigen提供了两个运算符：
- '*', 二元运算符，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
细心的读者可能会对"mat = mat*mat"产生疑问，按照Eigen的'expression templates'机制，也存在类似于例程3中的'aliasing issues'。 实际上，Eigen对矩阵乘法做了特殊的处理引入了一个临时变量，所以 'm = m*m'等价于：
tmp = m*m;
m = tmp;
如果我们可以确定不存在'aliasing issue'的问题时，可以通过方法'noallias()'不引入临时变量，如下：
c.noalias() += a * b;
关于aliasing issue的问题，我们会在后续章节中予以详细介绍。
2.4 点乘和叉乘运算
我们可以通过dot()和cross()方法实现点乘和叉乘运算。当然点乘运算也可以看作时一个1×1的矩阵u.ajoint()*v。 需要注意的时叉乘运算只针对维度为3的向量，点乘则适用于任何维度的向量。在处理复数时...
#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;
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
2.5 基本的arithmetic reduction运算
Eigen也提供了一些reduction运算，把一个矩阵后者向量reduce到一个数， 例如求所有元素的和'sum()'，最大值'maxCoeff()'，最小值'minCoeff()'等。 此外，Eigen还提供了用来求取minCoeff和maxCoeff元素的索引的方法：
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
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;
cout << "\n====\n" << endl;
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 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
====
Here is the matrix m:
0.680375 0.59688 -0.329554
-0.211234 0.823295 0.536459
0.566198 -0.604897 -0.444451
Its minimum coefficient (-0.604897) is at position (2,1)
Here is the vector v: 115899597 -48539462 276748203 -290373134
Its maximum coefficient (276748203) is at position 2