Eigen库：（一）密集矩阵和数组操作

1. Matrix类

在 Eigen 中，向量和矩阵都属于 Matrix 类。

1.1 矩阵和向量的定义

// 模板
Matrix<数据类型, 行数, 列数>
// 行数和列数可以是静态的，也可以是动态的
Matrix<float, 3, 3>
Matrix<float, 3, Dynamic>

// Matrix类typedef了很多常用类型
// 静态矩阵和向量
typedef Matrix<float, 4, 4> Matrix4f;
typedef Matrix<float, 3, 1> Vector3f;
typedef Matrix<int, 1, 2> RowVector2i;
// 动态矩阵和向量
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;

1.2 构造器

// 默认构造器，不分配内存，不初始化系数
Matrix3f a;
MatrixXf b;

// 分配内存（指定大小）的构造器，不初始化系数
MatrixXf a(10,15);
VectorXf b(30);

// 小向量构造器，初始化系数，最长为4
Vector2d a(5.0, 6.0);
Vector3d b(5.0, 6.0, 7.0);
Vector4d c(5.0, 6.0, 7.0, 8.0);

1.3 （）访问Matrix对象元素

// 矩阵 -> (行下标, 列下标)
MatrixXd m(2,2);
m(0,0) = 3;

// 向量 -> (下标)
VectorXd v(2);
v(0) = 4;

1.4 逗号初始化

Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;

1.5 属性：行数、列数、大小

MatrixXd m(2,5);
std::cout << "The matrix m is of size " << m.rows() << "x" << m.cols() << std::endl;
std::cout << "It has " << m.size() << " coefficients" << std::endl;

1.6 修改动态Matrix大小

// resize()方法
MatrixXd m(2,5);
m.resize(4,3);

// 赋值方法
MatrixXf a(2,2);
std::cout << "a is of size " << a.rows() << "x" << a.cols() << std::endl;
MatrixXf b(3,3);
a = b;
std::cout << "a is now of size " << a.rows() << "x" << a.cols() << std::endl;

2. 矩阵和向量运算

Eigen 不支持自动类型转换，因此我们对不同数据类型的Matrix进行运算的时候，必须作强制类型转换。

// 强制类型转换
Matrix<double,2,3> matrix_23 = Matrix<double,2,3>::Random();
Vector3f v_3d = Vector3f::Random();
Vector2f result = matrix_23.cast<float>() * v_3d;
cout<<result<<endl;

Eigen支持复合运算符。

2.1 加法和减法

// 必须对两个相同大小的矩阵或向量进行加减操作
binary operator + as in a+b
binary operator - as in a-b
unary operator - as in -a
compound operator += as in a+=b
compound operator -= as in a-=b

2.2 乘法和除法

binary operator * as in matrix*scalar
binary operator * as in scalar*matrix
binary operator / as in matrix/scalar
compound operator *= as in matrix*=scalar
compound operator /= as in matrix/=scalar

2.3 转置和共轭

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;

// 转置 -> transpose()
cout << "Here is the matrix a^T\n" << a.transpose() << endl;

// 共轭 -> conjugate()
cout << "Here is the conjugate of a\n" << a.conjugate() << endl;

// 共轭转置 -> adjoint()
cout << "Here is the matrix a^*\n" << a.adjoint() << endl;

这里有一个混叠问题需要注意：

// 当我们将矩阵的转置直接赋值给原矩阵时，发生混叠现象，报错
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;

# 因此我们对原矩阵进行转置时，用in-place的操作
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;

// There is also the adjointInPlace() function for complex matrices.

// 当我们将矩阵的转置直接赋值给另一个矩阵时，没有任何问题

2.4 矩阵-向量乘法 和 矩阵-矩阵乘法

binary operator * as in a*b
compound operator *= as in a*=b (this multiplies on the right: a*=b is equivalent to a = a*b)

2.5 内积和外积

Vector3d v(1,2,3);
Vector3d w(0,1,2);
cout << "Dot product: " << v.dot(w) << endl;;
cout << "Cross product:\n" << v.cross(w) << endl;

2.6 reduction操作

1、累加、累积、平均值、最小最大值、迹

#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;

    // minCoeff 和 maxCoeff 还可以接收位置参数
    Matrix3f m = Matrix3f::Random();
    std::ptrdiff_t i, j; // 这里可以直接用int来定义
    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;

}

2、范数

norm() // l2范数
squaredNorm() // l2范数的平方
lpNorm<p>() // lp范数

3、 Boolean reductions

#include <Eigen/Dense>
#include <iostream>
 
using namespace std;
using namespace Eigen;
 
int main()
{
  ArrayXXf a(2,2);
  
  a << 1,2,
       3,4;
 
  cout << "(a > 0).all()   = " << (a > 0).all() << endl;
  cout << "(a > 0).any()   = " << (a > 0).any() << endl;
  cout << "(a > 0).count() = " << (a > 0).count() << endl;
  cout << endl;
  cout << "(a > 2).all()   = " << (a > 2).all() << endl;
  cout << "(a > 2).any()   = " << (a > 2).any() << endl;
  cout << "(a > 2).count() = " << (a > 2).count() << endl;
}

结果：
(a > 0).all()   = 1
(a > 0).any()   = 1
(a > 0).count() = 4

(a > 2).all()   = 0
(a > 2).any()   = 1
(a > 2).count() = 2

4、Partial reductions

reduction 和colwise() or rowwise() 配合使用。

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
  std::cout << "Column's maximum: " << std::endl
   << mat.colwise().maxCoeff() << std::endl;
}

结果：
Column's maximum: 
3 2 7 9

2.7 广播机制

广播机制需要我们显式的实现。

#include <iostream>
#include <Eigen/Dense>
 
using namespace std;
int main()
{
  Eigen::MatrixXf mat(2,4);
  Eigen::VectorXf v(2);
  
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
         
  v << 0,
       1;
       
  //add v to each column of m
  mat.colwise() += v;
  
  std::cout << "Broadcasting result: " << std::endl;
  std::cout << mat << std::endl;
}

结果：
Broadcasting result: 
 1  3  8 12
 3  2  9  5

3. Array类和系数操作

由于Matrix类没有逐像素操作，不太方便，Array类正是为了实现逐像素操作的。

3.1 定义Array数组

// Array类和Matrix类一样的定义机制
Array<数据类型, 行数, 列数>

//typedef，与Matrix有一点区别
typedef Array<float,Dynamic,1> ArrayXf 
typedef Array<float,3,1> Array3f 
typedef Array<double,Dynamic,Dynamic> ArrayXXd 
typedef Array<double,3,3> Array33d 

3.2 （）访问Array数组元素

ArrayXXf  m(2,2);
m(0,0) = 1.0;
m(0,1) = 2.0;

3.3 逗号初始化

ArrayXXf a(3,3);
a << 1,2,3,
     4,5,6,
     7,8,9;

3.4 加法和减法

// 两个同样大小Array的加减法，和Matrix一样
a+b
a-b

// Array和Scalar的加减法，广播机制
a + 2
2 - a

3.5 逐像素乘法

// 与Matrix不一样，*是逐像素乘法，不是矩阵乘法
#include <Eigen/Dense>
#include <iostream>
 
using namespace Eigen;
using namespace std;
 
int main()
{
  ArrayXXf a(2,2);
  ArrayXXf b(2,2);
  a << 1,2,
       3,4;
  b << 5,6,
       7,8;
  cout << "a * b = " << endl << a * b << endl;
}

结果：
a * b = 
 5 12
21 32

3.6 其他逐像素操作

#include <Eigen/Dense>
#include <iostream>
 
using namespace Eigen;
using namespace std;
 
int main()
{
  ArrayXf a = ArrayXf::Random(5);
  a *= 2;
  cout << "a =" << endl 
       << a << endl;
  cout << "a.abs() =" << endl 
       << a.abs() << endl;
  cout << "a.abs().sqrt() =" << endl 
       << a.abs().sqrt() << endl;
  cout << "a.min(a.abs().sqrt()) =" << endl 
       << a.min(a.abs().sqrt()) << endl; //取两个相同大小数组的较小值
}

3.7 Matrix和Array的转换

// Matrix的array()方法：将matrix转换为array
// Array的matrix()方法：将array转换为matrix
// Matrix的cwiseProduct()方法：逐像素乘法
#include <Eigen/Dense>
#include <iostream>
 
using namespace Eigen;
using namespace std;
 
int main()
{
  MatrixXf m(2,2);
  MatrixXf n(2,2);
  MatrixXf result(2,2);
 
  m << 1,2,
       3,4;
  n << 5,6,
       7,8;
 
  result = m * n;
  cout << "-- Matrix m*n: --" << endl << result << endl << endl;
  result = m.array() * n.array();
  cout << "-- Array m*n: --" << endl << result << endl << endl;
  result = m.cwiseProduct(n);
  cout << "-- With cwiseProduct: --" << endl << result << endl << endl;
  result = m.array() + 4;
  cout << "-- Array m + 4: --" << endl << result << endl << endl;
}

4. 块操作

块操作可以作为左值和右值。

对于Matrix和Array都适用。

4.1 block操作

4.2 col() 和 row()

4.3 向量的块操作

5. 进阶初始化

5.1 逗号 + 块 初始化

// 向量的块初始化
RowVectorXd vec1(3);
vec1 << 1, 2, 3;
std::cout << "vec1 = " << vec1 << std::endl;
 
RowVectorXd vec2(4);
vec2 << 1, 4, 9, 16;
std::cout << "vec2 = " << vec2 << std::endl;
 
RowVectorXd joined(7);
joined << vec1, vec2;
std::cout << "joined = " << joined << std::endl;

// 矩阵的块初始化
MatrixXf matA(2, 2);
matA << 1, 2, 3, 4;
MatrixXf matB(4, 4);
matB << matA, matA/10, matA/10, matA;
std::cout << matB << std::endl;

5.2 特殊初始化

采用Matrix和Array的静态方法进行初始化。

Zero([size]);
Identity([size]);
LinSpaced(size, low, high); // 向量方法
Random([size])
Constant(size,val)

未完待续。。。