[Eigen中文文档] 扩展/自定义Eigen（二）

文档总目录

本文目录

- 使用nullary-expressions操作矩阵
- - 示例 1：循环矩阵
  - 示例 2：索引行和列

使用nullary-expressions操作矩阵

英文原文(Matrix manipulation via nullary-expressions)

CwiseNullaryOp 类的主要目的是定义过程矩阵，例如由Ones()、Zero()、Constant()、Identity()和Random()方法返回的常量或随机矩阵。然而，通过一些想象力，可以用最小的努力实现非常复杂的矩阵操作，因此很少需要实现新的表达式。

示例 1：循环矩阵

为了探索这些可能性，让我们从实现新表达式的循环矩阵示例开始。回想一下，循环矩阵是一个矩阵，其中每列与左侧的列相同，只是向下循环移位。例如，这是一个4x4的循环矩阵：
$\begin{bmatrix} 1&8&4&2 \\ 2&1&8&4 \\ 4&2&1&8 \\ 8&4&2&1 \end{bmatrix}$
循环矩阵由其第一列唯一确定。我们希望编写一个函数 makeCirculant，在给定第一列的情况下，返回一个表示循环矩阵的表达式。

在这个练习中，makeCirculant的返回类型将是一个CwiseNullaryOp，我们需要用以下内容进行实例化：

一个适当的circulant_functor，其中存储输入向量并实现适当的系数访问运算符(i，j)
一个矩阵类的模板实例化，传递编译时信息，如标量类型，大小和首选存储布局。

将 ArgType 称为输入向量的类型，我们可以构造等效的平方 Matrix 类型，如下所示：

template <class ArgType>
Eigen::CwiseNullaryOp<circulant_functor<ArgType>, typename circulant_helper<ArgType>::MatrixType>
makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
  typedef typename circulant_helper<ArgType>::MatrixType MatrixType;
  return MatrixType::NullaryExpr(arg.size(), arg.size(), circulant_functor<ArgType>(arg.derived()));
}

这个小辅助结构将帮助我们实现 makeCirculant 函数，如下所示：

template <class ArgType>
Eigen::CwiseNullaryOp<circulant_functor<ArgType>, typename circulant_helper<ArgType>::MatrixType>
makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
  typedef typename circulant_helper<ArgType>::MatrixType MatrixType;
  return MatrixType::NullaryExpr(arg.size(), arg.size(), circulant_functor<ArgType>(arg.derived()));
}

函数采用 MatrixBase 作为参数（有关更多详细信息，请参阅此页）。然后，通过 DenseBase::NullaryExpr 静态方法以足够的运行时大小构造 CwiseNullaryOp 对象。

接下来，需要实现 circulant_functor，这是一个简单的练习：

template<class ArgType>
class circulant_functor {
  const ArgType &m_vec;
public:
  circulant_functor(const ArgType& arg) : m_vec(arg) {}
 
  const typename ArgType::Scalar& operator() (Eigen::Index row, Eigen::Index col) const {
    Eigen::Index index = row - col;
    if (index < 0) index += m_vec.size();
    return m_vec(index);
  }
};

现在我们已经准备好尝试我们的新功能了：

int main()
{
    Eigen::VectorXd vec(4);
    vec << 1, 2, 4, 8;
    Eigen::MatrixXd mat;
    mat = makeCirculant(vec);
    std::cout << mat << std::endl;
}

如果组合所有片段，则会产生以下输出，表明程序按预期工作：

makeCirculant 的这种实现比从头开始定义新表达式要简单得多。

示例 2：索引行和列

这里的目标是模仿 MatLab，使用两个索引向量分别引用要选取的行和列来索引矩阵，如下所示：

A =
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 
A([1 2 1], [3 2 1 0 0 2]) =
 0  1 -6 -2 -2  1
 9  0 -3  6  6  0
 0  1 -6 -2 -2  1

为此，首先编写一个零元函数对象，存储对输入矩阵和两个索引数组的引用，并实现所需的operator() (i，j)：

template<class ArgType, class RowIndexType, class ColIndexType>
class indexing_functor {
  const ArgType &m_arg;
  const RowIndexType &m_rowIndices;
  const ColIndexType &m_colIndices;
public:
  typedef Eigen::Matrix<typename ArgType::Scalar,
                 RowIndexType::SizeAtCompileTime,
                 ColIndexType::SizeAtCompileTime,
                 ArgType::Flags&Eigen::RowMajorBit?Eigen::RowMajor:Eigen::ColMajor,
                 RowIndexType::MaxSizeAtCompileTime,
                 ColIndexType::MaxSizeAtCompileTime> MatrixType;
 
  indexing_functor(const ArgType& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
    : m_arg(arg), m_rowIndices(row_indices), m_colIndices(col_indices)
  {}
 
  const typename ArgType::Scalar& operator() (Eigen::Index row, Eigen::Index col) const {
    return m_arg(m_rowIndices[row], m_colIndices[col]);
  }
};

然后，创建一个 indexing(A,rows,cols) 函数来创建这个零元表达式(nullary expression)：

template <class ArgType, class RowIndexType, class ColIndexType>
Eigen::CwiseNullaryOp<indexing_functor<ArgType,RowIndexType,ColIndexType>, typename indexing_functor<ArgType,RowIndexType,ColIndexType>::MatrixType>
mat_indexing(const Eigen::MatrixBase<ArgType>& arg, const RowIndexType& row_indices, const ColIndexType& col_indices)
{
  typedef indexing_functor<ArgType,RowIndexType,ColIndexType> Func;
  typedef typename Func::MatrixType MatrixType;
  return MatrixType::NullaryExpr(row_indices.size(), col_indices.size(), Func(arg.derived(), row_indices, col_indices));
}

最后，这是一个如何使用该函数的示例：

Eigen::MatrixXi A = Eigen::MatrixXi::Random(4,4);
Eigen::Array3i ri(1,2,1);
Eigen::ArrayXi ci(6); ci << 3,2,1,0,0,2;
Eigen::MatrixXi B = mat_indexing(A, ri, ci);
std::cout << "A =" << std::endl;
std::cout << A << std::endl << std::endl;
std::cout << "A([" << ri.transpose() << "], [" << ci.transpose() << "]) =" << std::endl;
std::cout << B << std::endl;

这种简单的实现已经非常强大，因为行或列索引数组也可以是执行偏移、取模、跨步、反转等的表达式。

B =  mat_indexing(A, ri+1, ci);
std::cout << "A(ri+1,ci) =" << std::endl;
std::cout << B << std::endl << std::endl;
B =  mat_indexing(A, Eigen::ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), Eigen::ArrayXi::LinSpaced(4,0,3));
std::cout << "A(ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3)) =" << std::endl;
std::cout << B << std::endl << std::endl;

输出如下：

A(ri+1,ci) =
 9  0 -3  6  6  0
 9  3  6  6  6  3
 9  0 -3  6  6  0
 
A(ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), ArrayXi::LinSpaced(4,0,3)) =
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3
-2 -6  1  0
 6 -3  0  9
 6  6  3  9
 7  9 -5 -3