Eigen入门之密集矩阵 6 - Reductions, visitors and broadcasting

简介

本文介绍一下Dense Matrix的3中操作： reduction， visitor, broadcasting.

归约计算reduction.

Eigen的归约计算是这样的一类计算，它是对矩阵或者数组进行的计算，并返回一个标量。可能其中最常用的就是sum()计算，它返回的是系数的和。

比如,前面曾经用到的示例：

#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

范式计算（Norm comptation）

以下列出常用的向量范数和矩阵范数的定义，以及对应的MATLAB的函数：

1、向量范数

1-范数：，即向量元素绝对值之和，matlab调用函数norm(x, 1) 。
2-范数：，Euclid范数（欧几里得范数，常用计算向量长度），即向量元素绝对值的平方和再开方，matlab调用函数norm(x, 2)。
∞-范数：，即所有向量元素绝对值中的最大值，matlab调用函数norm(x, inf)。
-∞-范数：，即所有向量元素绝对值中的最小值，matlab调用函数norm(x, -inf)。
p-范数：，即向量元素绝对值的p次方和的1/p次幂，matlab调用函数norm(x, p)。

2、矩阵范数

1-范数：，列和范数，即所有矩阵列向量绝对值之和的最大值，matlab调用函数norm(A, 1)。
2-范数：，谱范数，即A’A矩阵的最大特征值的开平方。matlab调用函数norm(x, 2)。
∞-范数：，行和范数，即所有矩阵行向量绝对值之和的最大值，matlab调用函数norm(A, inf)。
F-范数：，Frobenius范数，即矩阵元素绝对值的平方和再开平方，matlab调用函数norm(A, ’fro‘)。

在Eigen中，提供了计算范式的一些方法。

squaredNorm 向量的平方范式，2-范式，ℓ2\ell^2ℓ2的范式，可以使用方法squaredNorm()，即所有系数的平方之和，数学上是一个向量与自身的点乘。
norm 即称欧拉范式，squaredNorm结果再开平方，得到平方根。

当这些计算作用于矩阵时，一个 nXp的矩阵，可以当做一个大小为（n*p）的向量，比如"Frobenius" or "Hilbert-Schmidt" norm。

我们避免提矩阵的ℓ2\ell^2ℓ2的范式（欧拉平方范式）运算，这样会导致歧义。

Eigen也额外提供了一些对系数进行欧拉范式操作方法，如ℓp\ell^pℓp范式，可以使用lpNorm<p>()方法。如果要计算ℓ∞\ell^\inftyℓ∞ 范式，模板参数p取特殊值Infinity
即可，这得到系数绝对值的最大值。

示例：

//matrix_norm1.cpp#include <Eigen/Dense>
#include <iostream>using namespace std;
using namespace Eigen;int main()
{VectorXf v(2);MatrixXf m(2,2), n(2,2);v << -1,2;cout << "vector v: "<<endl<<v<<endl;cout << "v.squaredNorm() = " << v.squaredNorm() << endl;cout << "v.norm() = " << v.norm() << endl;cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl;cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl;cout << endl;m << 1,-2,-3,4;cout << "matrix m: "<<endl<<m<<endl;cout << "m.squaredNorm() = " << m.squaredNorm() << endl;cout << "m.norm() = " << m.norm() << endl;cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << endl;cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Infinity>() << endl;
}

执行结果：

$ g++   -I /usr/local/include/eigen3 matrix_norm1.cpp -o matrix_norm1
$ ./matrix_norm1
vector v:
-12
v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2matrix m: 1 -2
-3  4
m.squaredNorm() = 30
m.norm() = 5.47723
m.lpNorm<1>() = 10
m.lpNorm<Infinity>() = 4

boolean 归约

Eigen提供了一些针对Boolean运算的归约操作。

all() : 如果系数全部满足Boolean计算结果，则结果为true；否则为false。
any() : 如果系数有任何满足Boolean计算的结果，则结果为true；否则为false。
count() : 返回满足条件的系数的数。

给一个示例：

#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

访问子 visitor

访问子visitor用于访问矩阵或数组内的系数及系数的位置。最简单的是前面曾经使用到的maxCoeff(&x,&y) , minCoeff(&x,&y)。

示例请看下面的Partial reduction中的例子程序。

Partial reductions 局部归约

Eigen提供了针对矩阵的列或者行的归约操作，它们对矩阵的每行/每列进行操作，并返回结果向量。
对应的函数为colwise() , rowwise()。
Note colwise() 返回的是一个行向量，而rowwise()返回的是一个列向量。

示例：

//matrix_norm2.cpp#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;std::cout << "Row's maximum: " << std::endl<< mat.rowwise().maxCoeff() << std::endl;}

执行：

$ g++   -I /usr/local/include/eigen3 matrix_norm2.cpp -o matrix_norm2
$ ./matrix_norm2
Column's maximum:
3 2 7 9
Row's maximum:
9
7

broadcasting

广播broadcasting比Partial reduction多了一个操作。其把一个向量赋值多次，结合成一个矩阵，可以和以及矩阵进行运算。

下面是一个简单的示例程序：

//matrix_broadcast.cpp
#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 mmat.colwise() += v;std::cout << "Broadcasting result1: " << std::endl;std::cout << mat << std::endl;// -------------Eigen::MatrixXf mm(2,2);mm<<1,2,3,4;v << 1,1;mm.rowwise() += v.transpose();std::cout << "Broadcasting result2: " << std::endl;std::cout << mm << std::endl;}

执行的情况类似：
[12693172]+[00001111]=[12694283].\begin{bmatrix} 1 & 2 & 6 & 9 \\ 3 & 1 & 7 & 2 \end{bmatrix} + \begin{bmatrix} 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 2 & 6 & 9 \\ 4 & 2 & 8 & 3 \end{bmatrix}. [13216792]+[01010101]=[14226893].

[1234]+[1111]=[2345].\begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix} + \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} = \begin{bmatrix} 2 & 3 \\ 4 & 5 \end{bmatrix}. [1324]+[1111]=[2435].

执行结果：

$ g++   -I /usr/local/include/eigen3 matrix_broadcast.cpp -o matrix_broadcast
$
$ ./matrix_broadcast
Broadcasting result1:
1 2 6 9
4 2 8 3
Broadcasting result2:
2 3
4 5

综合：结合broadcasting和其他的操作及运算

//matrix_broadcast2.cpp
#include <iostream>
#include <Eigen/Dense>using namespace std;
using namespace Eigen;int main()
{Eigen::MatrixXf m(2,4);Eigen::VectorXf v(2);m << 1, 23, 6, 9,3, 11, 7, 2;v << 2,3;MatrixXf::Index index;// find nearest neighbour(m.colwise() - v).colwise().squaredNorm().minCoeff(&index);cout << "Nearest neighbour is column " << index << ":" << endl;cout << "Nearest column vector: " << endl << m.col(index) << endl;
}

计算m.colwise() - v：
[1236931172]−[22223333]=[−12147084−1]\begin{bmatrix} 1 & 23 & 6 & 9 \\ 3 & 11 & 7 & 2 \end{bmatrix} - \begin{bmatrix} 2 & 2 & 2 & 2 \\ 3 & 3 & 3 & 3 \end{bmatrix} = \begin{bmatrix} -1 & 21 & 4 & 7 \\ 0 & 8 & 4 & -1 \end{bmatrix} [1323116792]−[23232323]=[−10218447−1]
计算(m.colwise() - v).colwise().squaredNorm():

(m.colwise() - v).colwise().squaredNorm() = [ \begin{bmatrix} 1 & 505 & 32 & 50 \end{bmatrix} ]

minCoeff(&index) 用来获得在m矩阵内，最小的系数位置索引，该列向量与向量v具有最短欧拉距离。

执行结果：

$ g++   -I /usr/local/include/eigen3 matrix_broadcast2.cpp -o matrix_broadcast2
$ ./matrix_broadcast2
Nearest neighbour is column 0:
Nearest column vector:
1
3

[431311431311431311]\begin{bmatrix} 4 & 3 & 13 & 11 \\ 4 & 3 & 13 & 11\\4 & 3 & 13 & 11 \end{bmatrix} ⎣⎡444333131313111111⎦⎤

[431311]\begin{bmatrix} 4 & 3 & 13 & 11 \end{bmatrix} [431311]