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

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

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

Eigen也额外提供了一些对系数进行欧拉范式操作方法，如${\ell }^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: 
-1
 2
v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2

matrix 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 m
    mat.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;

}

执行的情况类似：
$$\left[\begin{array}{cccc}1& 2& 6& 9\\ 3& 1& 7& 2\end{array}\right]+\left[\begin{array}{cccc}0& 0& 0& 0\\ 1& 1& 1& 1\end{array}\right]=\left[\begin{array}{cccc}1& 2& 6& 9\\ 4& 2& 8& 3\end{array}\right].$$

$$\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]+\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right]=\left[\begin{array}{cc}2& 3\\ 4& 5\end{array}\right].$$

执行结果：

$ 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：
  $$\left[\begin{array}{cccc}1& 23& 6& 9\\ 3& 11& 7& 2\end{array}\right]-\left[\begin{array}{cccc}2& 2& 2& 2\\ 3& 3& 3& 3\end{array}\right]=\left[\begin{array}{cccc}-1& 21& 4& 7\\ 0& 8& 4& -1\end{array}\right]$$

• 计算(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

$$\left[\begin{array}{cccc}4& 3& 13& 11\\ 4& 3& 13& 11\\ 4& 3& 13& 11\end{array}\right]$$

$$\left[\begin{array}{cccc}4& 3& 13& 11\end{array}\right]$$