Манипулирование матрицами с помощью нульарных выражений

Основное назначение класса 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, который нам нужно создать с помощью: 1 - подходящего circulant_functor для хранения входного вектора и реализации соответствующего доступа к коэффициентам operator(i,j) 2 - экземпляра шаблона класса Matrix, предоставляющего информацию о времени компиляции, такую как тип скаляра, размеры и предпочтительный порядок хранения.

Называя ArgType типом входного вектора, мы можем создать эквивалентный квадратный тип Matrix следующим образом:

template<class ArgType>
struct circulant_helper {
  typedef Matrix<typename ArgType::Scalar,
                 ArgType::SizeAtCompileTime,
                 ArgType::SizeAtCompileTime,
                 ColMajor,
                 ArgType::MaxSizeAtCompileTime,
                 ArgType::MaxSizeAtCompileTime> MatrixType;
};

Эта небольшая вспомогательная структура поможет нам реализовать нашу функцию makeCirculant следующим образом:

template <class ArgType>
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 (см. эту страницу для получения более подробной информации). Затем объект CwiseNullaryOp создается статическим методом DenseBase::NullaryExpr с соответствующими размерами во время выполнения.

Затем нам нужно реализовать наш 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() (Index row, Index col) const {
    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 Matrix<typename ArgType::Scalar,
                 RowIndexType::SizeAtCompileTime,
                 ColIndexType::SizeAtCompileTime,
                 ArgType::Flags&RowMajorBit?RowMajor: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() (Index row, Index col) const {
    return m_arg(m_rowIndices[row], m_colIndices[col]);
  }
};

Затем давайте создадим функцию indexing(A,rows,cols), создающую нульарное выражение:

template <class ArgType, class RowIndexType, class ColIndexType>
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);
  Array3i ri(1,2,1);
  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;
#if EIGEN_COMP_CXXVER >= 11
  B =  mat_indexing(A, ArrayXi::LinSpaced(13,0,12).unaryExpr([](int x){return x%4;}), 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;
#endif

и результатом является:

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