Eigen学习（七）规约、迭代和广播

规约

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

矩阵的迹是对角线元素的和，也可以用a.diagnoal().sun()方式求得。

范数计算

向量的平方范数由squaredNorm()获得，等价于向量对自身做点积，也等同于所有元素额平方和。Eigen也提供了norm()范数，返回的是squaredNorm()的根。这些操作也适用于矩阵。如果想使用其他元素级的范数，使用lpNorm<p>()方法，当求无穷范数时，模板参数p可以取特殊值Infinity，得到的是所有元素的最大绝对值。

#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;
  
  m << 1,-2,
       -3,4;
  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;
  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;
}

v.squaredNorm() = 5
v.norm() = 2.23607
v.lpNorm<1>() = 3
v.lpNorm<Infinity>() = 2

m.squaredNorm() = 30
m.norm() = 5.47723
m.lpNorm<1>() = 10
m.lpNorm<Infinity>() = 4

矩阵的1范数和无穷范数也可以用下面的方法计算：

#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
  MatrixXf m(2,2);
  m << 1,-2,
       -3,4;
  cout << "1-norm(m)     = " << m.cwiseAbs().colwise().sum().maxCoeff()
       << " == "             << m.colwise().lpNorm<1>().maxCoeff() << endl;
  cout << "infty-norm(m) = " << m.cwiseAbs().rowwise().sum().maxCoeff()
       << " == "             << m.rowwise().lpNorm<1>().maxCoeff() << endl;
}

1-norm(m)     = 6 == 6
infty-norm(m) = 7 == 7

布尔规约

如下的操作得到的是布尔值

• all()返回真，如果矩阵或数组的所有元素为真
• any()返回真，如果矩阵或数组至少有一个元素为真
• 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

迭代

当需要获得元素在矩阵或数组中的位置时使用迭代。

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
  Eigen::MatrixXf m(2,2);
  
  m << 1, 2,
       3, 4;
  //get location of maximum
  MatrixXf::Index maxRow, maxCol;
  float max = m.maxCoeff(&maxRow, &maxCol);
  //get location of minimum
  MatrixXf::Index minRow, minCol;
  float min = m.minCoeff(&minRow, &minCol);
  cout << "Max: " << max <<  ", at: " <<
     maxRow << "," << maxCol << endl;
  cout << "Min: " << min << ", at: " <<
     minRow << "," << minCol << endl;
}

Max: 4, at: 1,1
Min: 1, at: 0,0

这两个函数同样返回最小或最大值。

部分规约

部分规约指的是对矩阵或数组按行或列进行的操作，对每一列或者行进行规约操作时得到的是一个列或者行向量。如下例子得到矩阵每一列的最大值并存入一个行向量中

#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

同样也可以得到每一行的最大值，返回一个列向量

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

Row's maximum: 
9
7

部分规约和其他操作的结合

使用部分规约操作得到的结果去做其他的操作也是可以的，如下例子用于得到矩阵中元素和最大的一列

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
  MatrixXf mat(2,4);
  mat << 1, 2, 6, 9,
         3, 1, 7, 2;
  
  MatrixXf::Index   maxIndex;
  float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);
  
  std::cout << "Maximum sum at position " << maxIndex << std::endl;
  std::cout << "The corresponding vector is: " << std::endl;
  std::cout << mat.col( maxIndex ) << std::endl;
  std::cout << "And its sum is is: " << maxNorm << std::endl;
}

Maximum sum at position 2
The corresponding vector is: 
6
7
And its sum is is: 13

通过colwise()迭代，应用sum()对每一列规约操作得到一个新的1x4的矩阵，因此如果

那么

最终执行maxCoeff()操作获得元素和最大的列的索引。

广播

广播的概念类似于部分规约，不同之处在于广播通过对向量在一个方向上的复制，将向量解释成矩阵。如下例子将一个列向量加到矩阵的每一列中

#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 2 6 9
4 2 8 3

可以将mat.colwise()+=v理解成两种等价的方式，它将列向量加到矩阵的每一列中；或者将列向量复制4次的得到一个2x4的矩阵，之后进行矩阵的相加运算：

+=、+和-运算符也可以按列或行操作。在数组中也可以用*=、/=、和/运算符执行元素级的按行或列乘除运算。但不能用在矩阵上，如果想用v(0)乘以矩阵的第0列，v(1)乘以矩阵的第1列...使用mat = mat*v.asDiagonal()。

结合广播和其他操作

广播也可以和其他操作结合，比如矩阵或数组的运算、规约和部分规约操作。下面介绍一个更加复杂的例子，演示了在矩阵中找到和给定向量最接近的一列，使用到了欧氏距离。

#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 << m.col(index) << endl;
}

Nearest neighbour is column 0:
1
3

其中 

(m.colwise() - v).colwise().squaredNorm().minCoeff(&index);

这句话做的工作是：

（1）m.colwise()-v是一个广播操作，矩阵m的每一列减去v，得到一个新的矩阵

（2）(m.colwise() - v).colwise().squareNorm()是部分规约操作，按列计算矩阵的平方范数，得到一个行向量

（3）最终minCoeff(&index)根据欧氏距离获得矩阵中最接近v的一列的索引。