使用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;
}
如果组合所有片段，则会产生以下输出，表明程序按预期工作：
1 8 4 2
2 1 8 4
4 2 1 8
8 4 2 1
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
