Eigen入门教程（3）矩阵的算术运算

myleader

矩阵的四则运算

Matrix类重载了+、-、*、/、这4个运算符以及+=、-=、*=、/=这4个运算符，使其符合矩阵运算的规则。两个矩阵只要尺寸匹配都可以完成+-*的操作，而除法要求矩阵有逆，否则会返回一个由null组成的矩阵，但是不空。

1. int size=4；
   
2. MatrixXf A = MatrixXf::Random(size,size);
   
3. MatrixXf B = MatrixXf::Random(size,size);
   
4. Z = A*B;

复制代码

要注意的是，Eigen的矩阵类+-操作不支持标量，这一点和MATLAB不同。而Eigen的矩阵类支持和标量的乘法。也就是说如果你希望一个矩阵都加上5，你需要

1. A=A+Matrix<double, A.rows(), A.cols()>::Constant(5);

复制代码

另外Eigen也支持简写

1. A+=B;  // equal to A=A+B

复制代码

对于加减法，这没有太多好说的。但是对于乘法，实际上Eigen必须要在内存中额外开辟一块内存空间用以存放临时变量，计算完成后再复制回来，那么这种简写操作非但不能改善性能、节约内存，反倒会带来性能下降和内存占用量的增大甚至内存的碎片化，所以楼主不建议对乘法使用这种方式。同样的，也不建议把矩阵乘法的结果返回原矩阵，Eigen手册中提到可以使用noalias()函数来处理这个问题。

1. c.noalias() += a * b;

复制代码

不过楼主本人是不用的，觉得这个地方内存管理不太好弄。

对于除法Eigen不如MATLAB，因为它没有提供左除，而仅支持右除。

1. C=A/B;  //ok
   
2. C=A\B;  //can't compile

复制代码

如果你需要左除的话，可以使用求逆运算，Eigen提供了一个函数来求矩阵的逆inverse()

1. C=A*B.inverse();    // equal to C=A/B
   
2. C=A.inverse()*B;    // equal to MATLAB C=A\B

复制代码

和乘法类似，不推荐简写计算符，也不推荐返回原矩阵。

关于求逆，下一篇准备介绍方程求解，到时候再详细讲。

一些常用运算

点乘与叉乘

点乘就是数量积乘法，叉乘就是向量积乘法，这两个算法是针对矢量的，不适用于矩阵。Eigen提供了2个函数来完成dot()和cross()，使用起来非常简单。

1. double x=vec1.dot(vec2);
   
2. vec1=vec1.cross(vec2);

复制代码

点乘返回的是一个标量，不可能和矢量自身共用一个变量名，但是叉乘返回的是矢量，可以和参与计算的原始矢量共用一个变量名，这就涉及到一些内存管理的问题，Eigen已经解决了这个问题，即使你让叉乘返回原始矢量也是没问题的。实际上Eigen是在内存中开辟了临时区域用于计算，完成后返回等号左边的变量，毕竟叉乘仅能支持3个元素的矢量，内存占用不大。所以说不论使用者是否使用原始变量作为返回值，内存和计算的消耗都是差不多的，使用者也就不必在意了。后面讲到转置还有内存管理的问题。

求模

求矢量的模是很常用的，Eigen提供了norm()函数

1. double x=vec1.norm();

复制代码

norm()函数会把所有的元素都乘方、求和、再开方。如果是复数矢量，那么就等于所有元素的实部、虚部都乘方、求和、再开方。

单位化

有时候求模之后，需要把矢量按比例缩放成单位矢量，这个可以用除以模的方法实现，不过代码会比较难看，Eigen提供了单位化函数normalize()

1. vec2=vec1.normalize();

复制代码

因为这个操作返回的是矢量，所以可以返回给原矢量，这就涉及到内存管理的问题。如果你需要直接在原矢量上操作，就必须使用normalized()函数，注意和上面的函数比多了一个d，而且不需要返回值。

1. vec1.normalized();

复制代码

平方和

考虑到平方和也是常用的，而且求平方和也是求模的一个步骤，Eigen对矢量也提供了求平方和的函数squaredNorm()

1. double x=vec1.squaredNorm();

复制代码

注意不要使用norm()再平方的算法，那样只是平白浪费计算罢了。

多次方和
多次方和楼主觉得不是很常用，不过Eigen官方手册却把它列为常用，那就为大家介绍一下吧

1. double x=vec1.lpNorm<p>();
   
2. double x=vec1.lpNorm<Infinity>();
   

复制代码

其中p是次数，而且必须是常数，不允许变量，如果你需要变量的，请看下面的逐元素操作。另外变量p还支持Infinity内置变量，也就是你所使用的CPU和编译器所支持的最大值，这个东西在数学理论上可以用来计算无穷小的无穷大次方，但是在计算机数值运算上似乎没什么用处，笔者觉得唯一的用处就是看这个最大值是奇数还是偶数，比如说你把vec1的变量设置为-1，看最后的结果是-1还是+1


求和

求和也是很常用的，Eigen提供了sum()函数

1. double x=mat1.sum();

复制代码

Eigen的sum()和MATLAB的不同，MATLAB的sum只进行一层求和，而Eigen的sum()却是直接把整个矩阵全部元素求和。

求平均数

求平均数其实很简单，Eigen提供了mean()

1. double x=mat1.mean();

复制代码

求积
求积也是常用的运算，Eigen里面用prod()

1. double x=mat1.prod();

复制代码

求最小值和最大值
最小值和最大值也算是很常用了，Eigen提供了minCoeff()和maxCoeff()，这两个函数还支持2个参数，这两个参数必须是指向整数的指针，被指向的整数会被设定为最大值（或最小值）在矩阵中的行、列，也就是位置。如果是矢量，那就只返回一个位置参数。

1. s = R.minCoeff()              // min(R(:))
   
2. s = R.maxCoeff()              // max(R(:))
   
3. s = R.minCoeff(&r, &c)    // [aa, bb] = min(R); [cc, dd] = min(aa);
   
4.                           // r = bb(dd); c = dd; s = cc
   
5. s = R.maxCoeff(&r, &c)    // [aa, bb] = max(R); [cc, dd] = max(aa);

复制代码

求最大最小值仅支持实数矩阵，不支持复数矩阵

求矩阵的迹

迹是特征值的和，也是主对角元素的和，求特征值是很复杂的计算，所以求迹最好还是直接求主对角元素的和比较简单，Eigen提供了trace()。至于特征值的计算，那就复杂了，这里先不讲，以后再慢慢添加内容。

1. double x=mat1.trace();

复制代码

上面说的都是矩阵的数值运算，下面讲一讲矩阵自身的操作。数值运算对矩阵、矢量都是只读的，除非返回参与计算的原始矩阵，否则不会涉及到内存管理的问题。而对矩阵自身的操作就要涉及到内存管理的内容，要稍微复杂一些。

逻辑运算

Eigen提供了all()和any()两个逻辑运算函数，返回bool值，用来判断是否全部为非零元素或者有任何非零元素。另外还提供了一个count()函数用来计算有多少个元素返回了真，在MATLAB中，可以对逻辑矩阵进行sum运算来计算真的个数，这里使用count()函数。

1. if((array1 > 0).all()) ... // if all coefficients of array1 are greater than 0 ...
   
2. if((array1 < array2).any()) ... // if there exist a pair i,j such that array1(i,j) < array2(i,j) ...
   
3. cout<<(array1 < array2).count()<<endl;

复制代码

这里使用了Array类，而非Matrix类，下面会讲Array类到底是什么

共轭

复数矩阵的共轭也是常用的，Eigen提供了conjugate()

1. matcx1=matcx1.conjugate();
   
2. matcx1.conjugate();
   

复制代码

共轭操作可以直接在原矩阵上完成，甚至不需要返回值，这不会带来更多的内存消耗或者更大的计算量

转置

转置是方阵常用的计算，Eigen提供了transpose()

1. mat2=mat1.transpose();

复制代码

和共轭不同，转置不是简单地加上正负号这么简单，需要在内存中调整数据的排列方式。看一下下面这段错误代码

1. Matrix2i a; a << 1, 2, 3, 4;
   
2. cout << "Here is the matrix a:\n" << a << endl;
   
3. a = a.transpose(); // !!! do NOT do this !!!
   
4. 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

很明显转置后的a(0,1)——别忘了序列从0开始——不是我们期待的结果。这里涉及到transpose()函数的内存操作，它是按存储顺序逐个去挪动元素的，所以如果你让转置返回原矩阵就会导致错误，所以必须要使用另外一个矩阵来保存转置后的结果。

如果你坚持使用原矩阵的内存空间来存储转置矩阵就必须使用transposeInPlace()函数

1. a.transposeInPlace();

复制代码

这样可以避免转置错误，会省一半的内存，但是会付出内存操作的代价，速度会慢。选择空间还是时间要由程序员自己来把握。

伴随矩阵

伴随矩阵笔者觉得并不常用，不过Eigen还是将其列为常用，并提供了adjoint()函数

1. b=a.adjoint();

复制代码

和转置类似Eigen也提供了一个adjointInPlace()函数，使用方法和优缺点也一样。

矩阵的逐元素运算、逐行运算、逐列运算

矩阵运算中对乘法和除法的定义并不是对应元素相乘或相除，Eigen重载的*/操作符也都是针对矩阵运算的，如果你需要逐元素运算，类似于MATLAB的.*、./可以使用其自带的函数

1. mat1.cwiseProduct(mat2)
   
2. mat1.cwiseQuotient(mat2)

复制代码

cwiseProduct()函数是逐元素相乘，cwiseQuotient()是逐元素相除，Eigen称其为Coefficient-wise operation。

除了这两个必须要重新定义的操作意外，Eigen还提供了更多的函数进行逐元素运算

1. mat1.cwiseMin(mat2)
   
2. mat1.cwiseMax(mat2)
   
3. mat1.cwiseAbs2()
   
4. mat1.cwiseAbs()
   
5. mat1.cwiseSqrt()

复制代码

如果是无参函数，那么你需要将其返回正确的类型，如果是有参函数，那么实参将会被赋值为返回值。

其中的Abs函数是求绝对值，对于复数而言是求复数的模，而Abs2函数是绝对值的平方，对于复数而言就是模的平方。也就是说如果是一个m*n的复矩阵进行了cwiseAbs()操作，它将返回一个m*n的实矩阵，而不再是复矩阵了。

另外Eigen还提供了逐行和逐列操作，通常只有配合更进一步的运算才会有意义

1. R.prod()                  // prod(R(:))
   
2. R.colwise().prod()          // prod(R)
   
3. R.rowwise().prod()          // prod(R, 2) or prod(R')'
   
4. R.all()                   // all(R(:))
   
5. R.colwise().all()         // all(R)
   
6. R.rowwise().all()         // all(R, 2)
   
7. R.any()                   // any(R(:))
   
8. R.colwise().any()         // any(R)
   
9. R.rowwise().any()         // any(R, 2)
   
10. R.colwise()+=vec1;
    
11. R.colwise().maxCoeff(&index);
    
12. cout<<R.col(index)<<endl;
    

复制代码

colwise().sum()的操作会返回一个RowVector，而rowwise().sum()操作会返回一个Vector，其他操作依此类推。


Array类

Array类其实和Matrix类差别不大，其主要差别就是，Array类的算数符号都是针对逐元素运算的，而非矩阵和矢量的运算法则。也就是说跟Matrix类的cwise操作差不多。

Array类的加减法和Matrix类一样都是逐元素的，乘除法就已经不同了

1. int size=4
   
2. ArrayXd arr1=ArrayXd::Random(size);
   
3. ArrayXd arr2=ArrayXd::Random(size);
   
4. ArrayXd arr3;
   
5. arr3=arr1+arr2;  //equal to mat1+mat2
   
6. arr3=arr1*arr2;  //equal to mat1.cwiseProduct(mat2)
   
7. arr3=arr1/arr2;  //equal to mat1.cwiseQuotient(mat2)

复制代码

对于Matrix类而言，任何一项逐元素操作都需要单独定义函数，这就使得Matrix类在面对逐元素运算时能力非常有限（你不可能预见到所有的逐元素操作并定义所有的函数，即使你能预见到，太多的成员函数也会带来效率的降低）。Array类支持更多的针对逐元素的运算。比如比较运算、三角函数等

1. array1 < array2     //array1 > array2     array1 < scalar     array1 > scalar
   
2. array1 <= array2    //array1 >= array2    array1 <= scalar    array1 >= scalar
   
3. array1 == array2    //array1 != array2    array1 == scalar    array1 != scalar
   
4. 
   
5. array1.min(array2)            
   
6. array1.max(array2)            
   
7. array1.abs2()
   
8. array1.abs() //std::abs(array1)
   
9. array1.sqrt()                 //std::sqrt(array1)
   
10. array1.log()                  //std::log(array1)
    
11. array1.exp()                  //std::exp(array1)
    
12. array1.pow(exponent)          //std::pow(array1,exponent)
    
13. array1.square()
    
14. array1.cube()
    
15. array1.inverse()
    
16. array1.sin()                  //std::sin(array1)
    
17. array1.cos()                  //std::cos(array1)
    
18. array1.tan()                  //std::tan(array1)
    
19. array1.asin()                 //std::asin(array1)
    
20. array1.acos()                 //std::acos(array1)
    

复制代码

确实这种操作可以用for循环来做，不过Array类可以调用OpenMP多线程进行加速，会稍微快一点点，所以还是推荐Array类。

Array类使用了模板技术来把只有一个输入参数的函数逐个映射到元素上，所以基本上只要是只有一个输入参数的数学函数都可以映射上去。只有乘方的运算还涉及到一个次数，不过因为它是一个常数，所以Array类内部已经把它处理掉了。

有一点遗憾atan2函数需要2个输入参数，所以没办法用到Array类上，大家编程时尽量避免atan2函数吧。如果实在避免不了，那就只能用for循环来做了。

Array和Matrix的转换
在实际使用中如果针对某个矩阵有时需要进行矩阵算法，有时需要逐元素的算法，就需要在Array类和Matrix类之间进行转换，很方便：

1. arr1=mat1.array();
   
2. mat1=arr1.matrix();

复制代码

因为这个转换会在内存中创建一个新的对象，所以对于大矩阵而言其实是非常耗费资源的，所以尽量使用cwise操作，而不要通过在Matrix和Array类之间的转换来完成。
矩阵的一些高级操作

矩阵的拼接

MATLAB中矩阵的拼接只要用方括号括起来就自动完成了，C++中稍微复杂点。如果是把A和B拼接后的结果存入C的话，只要创建一个大一点的矩阵C，然后用block()函数操作C的子矩阵将其赋值为A和B就好了。如果是把B拼接到A后面的话，可以先把A扩充，扩充部分用block()函数赋值为B就可以。如果是把B拼接到A的前面，这个真没办法，只能先把A和B拼接为C，然后把C赋值给A才行。在这3种方式中，楼主最推荐的是第1种，其次是第2种，实在不行才用第3种。

下面代码演示了在行方向拼接的方法

1. //conbine A and B to C
   
2. MatrixXd C(A.rows(), A.cols()+B.cols());
   
3. C.leftCols(A.cols())=A;
   
4. C.rightCols(B.cols())=B;
   
5. 
   
6. //append B to A
   
7. A.conservativeResize(Eigen::NoChange, A.cols()+B.cols());
   
8. A.rightCols(B.cols())=B;
   
9. 
   
10. //prepend B to A
    
11. MatrixXd C(A.rows(), A.cols()+B.cols());
    
12. C.leftCols(A.cols())=A;
    
13. C.rightCols(B.cols())=B;
    
14. A=C;

复制代码

复制矩阵 repmat

有时候我们需要把一个小矩阵复制好多次然后拼成一个大矩阵以利用矢量加速完成计算，其实就是用空间换时间，至于时间和空间的平衡要由程序员自己把握。Eigen提供了replicate

1. MatrixXi m = MatrixXi::Random(2,3);
   
2. cout << "Here is the matrix m:" << endl << m << endl;
   
3. cout << "m.replicate<3,2>() = ..." << endl;
   
4. cout << m.replicate<3,2>() << endl;
   
5. cout << m.replicate(3,2) << endl;

复制代码

这个函数的功能类似于MATLAB中的repmat，和block函数类似，replicate()提供了2种调用方法，尖括号版本要求复制倍数是常数，这样在编译时就可以优化。圆括号版本不能在编译时优化，但是有更好的适应性。括号里第1个参数是行方向复制的倍数，第2个参数是列方向的复制倍数。上面这段代码的输出结果是

Here is the matrix m:
7  6  9
-2  6 -6
m.replicate<3,2>() = ...
7  6  9  7  6  9
-2  6 -6 -2  6 -6
7  6  9  7  6  9
-2  6 -6 -2  6 -6
7  6  9  7  6  9
-2  6 -6 -2  6 -6

网格矩阵 meshgrid
MATLAB中提供了meshgrid函数生成网格点阵的矩阵，如果是Eigen的话，可以这样

1. X = RowVectorXd::LinSpaced(3,1,3).replicate(5,1);
   
2. Y = VectorXd::LinSpaced(5,10,14).replicate(1,3);

复制代码

重新排序 reshape
有些时候我们需要把矩阵改变一下排序，比如说2*3的矩阵改成3*2的，但是里面的元素还是按顺序存储。MATLAB中提供了reshape函数，Eigen并未提供这项功能，不过我们可以通过Map来实现。

1. B = Map<MatrixXd>(A.data(),m,n);

复制代码

repmat、meshgrid和reshape功能是楼主在实际使用中摸索出来的，官方教程隐藏在非常隐蔽的地方，这里写出来给大家一点方便。

到这里一般的矩阵运算应该都清楚了，下一篇计划内容是方程求解，敬请期待

caoer

写的真好，期待期待

caoer

interested in the perforamce of array and vector. so did following test.

start vector adding ...
time taken = 0.029
start array adding ...
time taken = 0.036
请按任意键继续. . .

1. // check the computing speed between matrix and array
   
2. 
   
3. #include <iostream>
   
4. #include <ctime>
   
5. #include "../include/Eigen/Core"
   
6. 
   
7. using namespace std;
   
8. using namespace Eigen;
   
9. 
   
10. int main()
    
11. {
    
12.     int size=2048*16, N=1000;
    
13. 
    
14.     clock_t start, end;
    
15. 
    
16.     VectorXd A = VectorXd::Random(size);
    
17.     VectorXd B = VectorXd::Random(size);
    
18.     VectorXd Y;
    
19. 
    
20.     cout <<"start vector adding ..."<<endl;
    
21.     start = clock();
    
22.     for(int i=0; i<N; ++i)
    
23.         Y = A+B;  //  or Z = A+B+C ... etc
    
24.     end = clock();
    
25.     cout << "time taken = " << 1.0*(end-start)/CLOCKS_PER_SEC << endl;
    
26. 
    
27. 
    
28.     ArrayXd C = ArrayXd::Random(size);
    
29.     ArrayXd D = ArrayXd::Random(size);
    
30.     ArrayXd Z;
    
31.     cout <<"start array adding ..."<<endl;
    
32.     start = clock();
    
33.     for(int i=0; i<N; ++i)
    
34.         Z = C+D;  //  or Z = A+B+C ... etc
    
35.     end = clock();
    
36.     cout << "time taken = " << 1.0*(end-start)/CLOCKS_PER_SEC << endl;
    
37. 
    
38. 
    
39.     system("pause");
    
40.     return 0;
    
41. }
    

复制代码

myleader

caoer 发表于 2012-12-9 09:51
interested in the perforamce of array and vector. so did following test.

start vector adding ...

Matrix和Array的速度我还真没对比过，不过你的代码有点问题，VectorXd生成的是Matrix(size,1)，而ArrayXd生成的是Array(size,size)，也就是说Array类的尺寸更大、数据更多，所以才会导致计算时间变长的。

我在奔腾G620上的实验，Matrix类耗时0.046s，Array类耗时0.047秒，相差不是很大。

caoer

thanks, i thought array is a one D vector.

hoop247

要注意的是，Eigen的矩阵类+-操作不支持标量，这一点和MATLAB不同。而Eigen的矩阵类支持和标量的乘法。也就是说如果你希望一个矩阵都加上5，你需要

A=A+Matrix<double, A.rows(), A.cols()>::Constant(5);
复制代码


这里也可以换成

A=A.array()+5;

Matrix类是不支持，但是Array类支持，可以用 .array() 转换matrix到array再加
-----------------------------------------------------------------------------------------------

后面看到了array的用法，