Eigen3学习笔记——密集矩阵和数组操作(2)

矩阵和算数运算

Eigen 通过重载C++的算数运算符或特殊方法来进行矩阵的算数运算

加减

参与加减的矩阵要满足行列情况，并且元素数据类型也要相同，Eigen 不会自己进行数据类型转换。

//这里简单给出加减法运算，注意数据矩阵的数据类型
Matrix2d a;
a << 1, 2,
	3, 4;
MatrixXd b(2,2);
b << 2, 3,
	1, 4;
//双目加法：a+b
//双目减法：a-b
//单目减法：-a
//复合加法：a+=b			//a = a + b;
//复合减法：a-=b			//a = a - b;

标量与矩阵的乘除法

Matrix2d a;
a << 1, 2,
	3, 4;
//双目乘法：a*2.5
//双目乘法：2.5*a
//双目除法：a/2.5
//复合乘法：a*=2.5
//复合除法：a/=2.5

关于表达式模板的说明

在Eigen 中，算术运算符本身是不执行任何计算，他只返回一个描述要执行的“表达式对象”，实际的计算是编译器进行优化后进行的。

说白了，不要去想象矩阵运算是怎么进行的，是按照数学上直接相加还是怎么样，不关心也没必要去知道。

转置、共轭和伴随

注意：

1. transpose()并adjoint()简单地返回一个代理对象没有做实际的换位。
2. 尽量避免a = a.transpose()这样的操作，这种操作a 将不会发生变化，这种现象称为aliasing issue。

Matrix2d a;
a << 1, 2,
	3, 4;
Matrix2d b;
a << 5, 2,
	7, 5;
a.transpose();			//转置，a取转置，但a本身不变化
a.conjugate();			//共轭
a.adjoint();			//伴随
a.transposeInPlace();	//对a 进行立即取转置，a 直接产生变化
b = a.transpose();		//结果在写入的时候同时产生转置，之前的a.transpose()只是简单的返回代理对象，并没有进行转换
a = a.transpose();		//不允许这样做，a的值将不变化。在写入a的同时还要对a进行转置，这不难为人么

矩阵乘法

注意：

1. Eiegn将矩阵乘法当成特殊情况，不会出现aliasing issue，因为Eigen会添加一个临时变量temp = a*a

Matrix2d a;
a << 1, 2,
	3, 4;
Matrix2d b;
b << 2, 3,
	4, 5;
//双目乘法：a*b
//复合乘法：a*=b			//等价于a = a*b，注意a，b位置
//允许：a = a*a		 //可以这样写，不会产aliasing issue

点积和叉积

注意：

1. 点积使用任意大小向量；
2. 叉积只适用于大小为3的向量；

Vector3d v(1,2,3);
Vector3d w(0,1,2);
\\点积：v.dot(w);
\\叉积：v.cross(w):

Array类

Array 类与Matrix 类很相似，但是我认为：

1. Array 更加强调元素运算，元素对应元素。
2. Matrix 更加强调线性代数运算，也就是更加符合我们认知的数学意义。

Array的定义

Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>	//与Matrix一样

Array的操作

自带定义
Array<float,Dynamic,1> 			typedef ArrayXf 
Array<float,3,1> 				typedef Array3f 
Array<double,Dynamic,Dynamic> 	typedef ArrayXXd 
Array<double,3,3> 				typedef Array33d 
    
索引
ArrayXXf  m(2,2);
m(0,0) = 1.0;	//索引，与Matrix一样

加减法
ArrayXXf a(3,3);
ArrayXXf b(3,3);
a << 1,2,3,
	4,5,6,
	7,8,9;
b << 1,2,3,
	1,2,3,
	1,2,3;
//a + b 的结果是：2，4，6；5，7，8；8，10，12；对应元素进行操作！！！！
//a + 2 :Array类可以与标量进行运算

乘法
ArrayXXf a(2,2);
ArrayXXf b(2,2);
a << 1,2,
	3,4;
b << 5,6,
	7,8;
//a * b 的结果是：5，12；21，32；对应元素进行乘法！！！！

其他运算
a.abs()			//对每个元素取绝对值
a.sqrt()		//对每个元素取算术平方根

Array与Matrix之间的转换

MatrixXf m(2,2);
ArrayXXf n(2,2);
m.array();		//Matrix转Array
m.matrix();		//Array转Matrix

块操作

就是取Matrix或者Array的一个子块：

matrix.block(i,j,p,q);
matrix.block<p,q>(i,j);
array.block(i,j,p,q);
array.block<p,q>(i,j);
//i,j为起始行列;
//p,q块的大小；！！！！

行列块

matrix.row(i);		//取低i行
matrix.col(j);		//取低j列

其他块操作函数，针对矩阵

//取左上角p行q列的矩阵
matrix.topLeftCorner(p,q);
matrix.topLeftCorner<p,q>();
//取左下角p行q列的矩阵
matrix.bottomLeftCorner(p,q);
matrix.bottomLeftCorner<p,q>();
//取右上角p行q列的矩阵
matrix.topRightCorner(p,q);
matrix.topRightCorner<p,q>();
//取右下角p行q列的矩阵
matrix.bottomRightCorner(p,q);
matrix.bottomRightCorner<p,q>();
//取从顶部开始的q行
matrix.topRows(q);
matrix.topRows<q>();
//取从底部开始的q行
matrix.bottomRows(q);
matrix.bottomRows<q>();
//取从左部开始的p列
matrix.leftCols(p);
matrix.leftCols<p>();
//取从右部开始的p列
matrix.rightCols(q);
matrix.rightCols<q>();
//取从i行开始的q行
matrix.middleCols(i,q);
matrix.middleCols<q>(i);
//取从i列开始的q列
matrix.middleRows(i,q);
matrix.middleRows<q>(i);

其他的块操作函数，针对向量

//取向量前n个元素
vector.head(n);
vector.head<n>();
//取向量后n个元素
vector.tail(n);
vector.tail<n>();
//取从i开始的n个元素
vector.segment(i,n);
vector.segment<n>(i);

切片与索引(3.4引入)

对于Matrix 和Array 的索引不单单局限于之前提到的角标或者块的方式，Eigen 中支持利用**“切片进行索引”，可以在操作符”()“**，中添加切片完成对Matrix 或者Array 子集索引。很像matlab中的索引方式

这种利用**”操作符()“**进行的索引可以如下表示：

DenseBase::operator()(const RowIndices&, const ColIndices&) 	//官方文档的表示，只要是“()”操作，则都可以用切片索引

其中两个参数可以是如下的任何一种形式：

1. 单个整形数，可以是符号
2. Eigen::all，Eigen::last等表示的集合
3. Eigen::seq,Eigen::seqN,Eigen::lastN构造的集合
4. 任意一维向量或数组，Matrix/Array/C数组。

基本切片

seq (firstIdx,lastIdx);				//从firstIdx开始，每次递增1，到lastIdx结束
seq (firstIdx,lastIdx,incr);		//从firstIdx开始，每次递增incr，到lastIdx结束
seqN (firstIdx,size);				//从firstIdx开始，每次递增1，递增size个
seqN (firstIdx,size,incr);			//从firstIdx开始，每次递增incr，递增size个

Eigen::all;							//从0开始每次递增1，直到结束
Eigen::last;							//表示结尾，不是集合

注：后面有些内容看不懂，是关于自定义切片类型的部分，感觉用不到，就不写了

高级初始化方式

利用逗号进行初始化

逗号初始化是：从左到右，从上到下，一行一行进行填充。

//用逗号直接构造一个矩阵
Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;

//逗号还可以用小矩阵来构建大矩阵
RowVectorXd vec1(3);	//小向量
vec1 << 1, 2, 3;
RowVectorXd vec2(4);	//小向量
vec2 << 1, 4, 9, 16;
RowVectorXd joined(7);	//大向量
joined << vec1, vec2;	//逗号构造！！！！

MatrixXf matA(2, 2);	//小矩阵
matA << 1, 2, 3, 4;
MatrixXf matB(4, 4);	//大矩阵
matB << matA, matA/10, matA/10, matA;	//逗号构造！！！

//逗号构建矩阵的子集
Matrix3f m;
m.row(0) << 1, 2, 3;
m.block(1,0,2,2) << 4, 5, 7, 8;
m.col(2).tail(2) << 6, 9;  

此外：

1. 有很多预定义好的函数可以直接得到特殊矩阵，比如单位矩阵等等，详细见官方快速参考手册
2. 可以利用特殊矩阵作为表达式来辅助构造需要的矩阵，这样不会产生额外开销。

//m数据结构：Matrix或者Array,一下是常量矩阵和向量，动态矩阵和向量需要指定行和列
m::Zero();							//零矩阵
m::Ones();							//全1矩阵
m::Constant(value);					//常量矩阵
m::Random();						//随机矩阵
m::LinSpaced(size, low, high);		//生成size大小，从low到high的向量，只能是向量！！！！
m::Identity();
//x是对象，Matrix或者Array对象
x.setZero();
x.setOnes();
x.setConstant(value);
x.setRandom();
x.setLinSpaced(size, low, high);
//利用表达式进行构造
MatrixXd m = MatrixXd::Random (3,3);
m = (m + MatrixXd::Constant (3,3,1.2)) * 50		//用函数生成的常量矩阵来构造m，减小开销。

Reductions, visitors and broadcasting

Reductions

就是：根据矩阵元素计算一个标量，比如矩阵的迹等等

Eigen::Matrix2d mat;
mat << 1, 2,
	3, 4;
mat.sum();			//矩阵元素求和
mat.prod();			//矩阵元素求积
mat.mean();			//矩阵元素就平均值
mat.minCoeff();		//最小元素
mat.maxCoeff();		//最大元素
mat.trace();		//矩阵的迹

Norm computations（范数计算）

MatrixXf m(2,2);
m << 1,-2,
	-3,4;
m.squaredNorm();			//二范数
m.norm();					//二范数的平方根
m.lpNorm<1>();				//p范数，<1>表示一范数
m.lpNorm<Infinity>();		//无穷范数

Boolean reductions

注：通常与Array的比较运算符联合使用，并且C++中非零都算是真哦

ArrayXXf a(2,2);
a << 1,2,
    3,4;
(a > 0).all();			//元素是否都为真
(a > 0).any();			//元素是否有至少一个真
(a > 0).count();		//返回真的个数
a.all();				//也可以，非零都算是真

Visitors

有时需要获取元素的位置信息，常用比如获取最大最小元素的位置和值，用到.maxCoeff和.minCoeff函数。

Eigen::MatrixXf m(2,2);
m << 1, 2,
	3, 4;

MatrixXf::Index maxRow, maxCol;					//Index类型，查找后是std::ptrdiff_t
float max = m.maxCoeff(&maxRow, &maxCol);		//max,maxRow,maxCol均是输入

MatrixXf::Index minRow, minCol;
float min = m.minCoeff(&minRow, &minCol);

Partial reductions

可以对矩阵的部分元素执行计算。

mat.colwise().maxCoeff();		//mat.colwize()本身返回的是VectorwiseOp类对象，maxCoeff()是这个类的函数
mat.rowwise().maxCoeff();		//直观理解就是将矩阵按每行或每列切片，个人理解啊！！
//官方给的一个复杂的应用
MatrixXf::Index   maxIndex;
float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);
//先把矩阵按每列切片，然后按行求和，此时得到的是一个向量了
//然后找向量中最大元素

Broadcasting

类似于Partial reductions，但是进行表达式的算术运算。

m << 1, 23, 6, 9,
	3, 11, 7, 2;
v << 2,
	3;
MatrixXf::Index index;
(m.colwise() - v).colwise().squaredNorm().minCoeff(&index);	
//!!!!! m.colwise() 返回的是类
//!!!!! m.colwise()-v 返回的是矩阵。。很有意思吧

注：这里是我根据源文件查看的信息，自己的理解，我尝试取直接输出 m.colwise()和 m.colwise()-v，前者会编译错误，后者输出矩阵。