Tutorial page 7 - Reductions, visitors and broadcasting

This tutorial explains Eigen's reductions, visitors and broadcasting and how they are used with matrices and arrays .

Table of contents

Reductions

In Eigen, a reduction is a function taking a matrix or array, and returning a single scalar value. One of the most used reductions is .sum() , returning the sum of all the coefficients inside a given matrix or array.

Example:Output: 
#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;
}

 




The trace of a matrix, as returned by the function trace(), is the sum of the diagonal coefficients and can equivalently be computed a.diagonal().sum().



Norm computations


The (Euclidean a.k.a. $\ell^2$) squared norm of a vector can be obtained squaredNorm() . It is equal to the dot product of the vector by itself, and equivalently to the sum of squared absolute values of its coefficients.


Eigen also provides the norm()  method, which returns the square root of squaredNorm() .


These operations can also operate on matrices; in that case, a n-by-p matrix is seen as a vector of size (n*p), so for example the norm()  method returns the "Frobenius" or "Hilbert-Schmidt" norm. We refrain from speaking of the $\ell^2$ norm of a matrix because that can mean different things.


If you want other $\ell^p$ norms, use the lpNnorm<p>()  method. The template parameter p can take the special value Infinity if you want the $\ell^\infty$ norm, which is the maximum of the absolute values of the coefficients.


The following example demonstrates these methods.





Example:Output: 




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

  
 





Boolean reductions


The following reductions operate on boolean values:



These are typically used in conjunction with the coefficient-wise comparison and equality operators provided by Array. For instance, array > 0 is an Array of the same size as array , with true at those positions where the corresponding coefficient of array is positive. Thus, (array > 0).all() tests whether all coefficients of array are positive. This can be seen in the following example:





Example:Output: 




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

  
 





User defined reductions


TODO


In the meantime you can have a look at the DenseBase::redux() function.



Visitors


Visitors are useful when one wants to obtain the location of a coefficient inside a Matrix or Array. The simplest examples are maxCoeff(&x,&y)  and minCoeff(&x,&y), which can be used to find the location of the greatest or smallest coefficient in a Matrix or Array.


The arguments passed to a visitor are pointers to the variables where the row and column position are to be stored. These variables should be of type Index  (FIXME: link ok?), as shown below:





Example:Output: 




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

  
 




Note that both functions also return the value of the minimum or maximum coefficient if needed, as if it was a typical reduction operation.



Partial reductions


Partial reductions are reductions that can operate column- or row-wise on a Matrix or Array, applying the reduction operation on each column or row and returning a column or row-vector with the corresponding values. Partial reductions are applied with colwise()  or rowwise() .


A simple example is obtaining the maximum of the elements in each column in a given matrix, storing the result in a row-vector:





Example:Output: 




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

  
 




The same operation can be performed row-wise:





Example:Output: 




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

  
 




Note that column-wise operations return a 'row-vector' while row-wise operations return a 'column-vector'



Combining partial reductions with other operations


It is also possible to use the result of a partial reduction to do further processing. Here is another example that aims to find the column whose sum of elements is the maximum within a matrix. With column-wise partial reductions this can be coded as:





Example:Output: 




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

  
 




The previous example applies the sum()  reduction on each column though the colwise()  visitor, obtaining a new matrix whose size is 1x4.


Therefore, if 



\[ \mbox{m} = \begin{bmatrix} 1 & 2 & 6 & 9 \\ 3 & 1 & 7 & 2 \end{bmatrix} \]



then



\[ \mbox{m.colwise().sum()} = \begin{bmatrix} 4 & 3 & 13 & 11 \end{bmatrix} \]



The maxCoeff()  reduction is finally applied to obtain the column index where the maximum sum is found, which is the column index 2 (third column) in this case.



Broadcasting


The concept behind broadcasting is similar to partial reductions, with the difference that broadcasting constructs an expression where a vector (column or row) is interpreted as a matrix by replicating it in one direction.


A simple example is to add a certain column-vector to each column in a matrix. This can be accomplished with:





Example:Output: 




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

  
 




It is important to point out that the vector to be added column-wise or row-wise must be of type Vector, and cannot be a Matrix. If this is not met then you will get compile-time error. This also means that broadcasting operations can only be applied with an object of type Vector, when operating with Matrix. The same applies for the Array class, where the equivalent for VectorXf is ArrayXf.


Therefore, to perform the same operation row-wise we can do:





Example:Output: 




#include <iostream>
#include <Eigen/Dense>

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

  
 





Combining broadcasting with other operations


Broadcasting can also be combined with other operations, such as Matrix or Array operations, reductions and partial reductions.


Now that broadcasting, reductions and partial reductions have been introduced, we can dive into a more advanced example that finds the nearest neighbour of a vector v within the columns of matrix m. The Euclidean distance will be used in this example, computing the squared Euclidean distance with the partial reduction named squaredNorm() :





Example:Output: 




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

  
 




The line that does the job is 


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

We will go step by step to understand what is happening:












re_vision

Author(s): Dorian Galvez-Lopez

autogenerated on Sun Jan 5 2014 11:33:48