Использование Eigen's expression templates приводит к тому, что каждый выражение может иметь разный тип. Если вы передаете такое выражение функции, принимающей параметр типа Matrix, ваше выражение неявно будет вычислено в временную Matrix, которая затем будет передана в функцию. Это означает, что вы теряете преимущества expression templates. Конкретно, это имеет два недостатка:

• Вычисление во временную переменную может быть бесполезным и неэффективным;
• Это позволяет функции только считывать данные из выражения, а не записывать в него.

К счастью, все эти многочисленные типы выражений имеют общие шаблонизированные базовые классы. Путем передачи в функцию шаблонизированных параметров этих базовых типов, вы можете позволить им работать вместе с expression templates Eigen.

Некоторые первые примеры

Этот раздел предоставит простые примеры для различных типов объектов, которые предоставляет Eigen. Прежде чем начать с фактических примеров, нам нужно повторить, с какими базовыми объектами мы можем работать (см. также Иерархия классов).

• MatrixBase: Общий базовый класс для всех плотных матричных выражений (в отличие от выражений массива, в отличие от разреженных и специальных классов матриц). Используйте его в функциях, предназначенных для работы только с плотными матрицами.
• ArrayBase: Общий базовый класс для всех плотных выражений массивов (в отличие от матричных выражений и т. д.). Используйте его в функциях, предназначенных для работы только с массивами.
• DenseBase: Общий базовый класс для всех плотных матричных выражений, то есть базовый класс как для MatrixBase , так и для ArrayBase. Он может использоваться в функциях, предназначенных для работы с матрицами и массивами.
• EigenBase: Базовый класс, объединяющий все типы объектов, которые могут быть вычислены в плотные матрицы или массивы, например, специальные классы матриц, такие как диагональные матрицы, перестановочные матрицы и т. д. Его можно использовать в функциях, предназначенных для работы с любым таким общим типом.

Пример EigenBase

Выводит размерности наиболее обобщенного объекта, присутствующего в Eigen. Он может быть любым матричным выражением, любой плотной или разреженной матрицей, и любым массивом.

#include <iostream>
#include <Eigen/Core>
using namespace Eigen;
 
template <typename Derived>
void print_size(const EigenBase<Derived>& b)
{
  std::cout << "size (rows, cols): " << b.size() << " (" << b.rows()
            << ", " << b.cols() << ")" << std::endl;
}
 
int main()
{
    Vector3f v;
    print_size(v);
    // v.asDiagonal() returns a 3x3 diagonal matrix pseudo-expression
    print_size(v.asDiagonal());
}

size (rows, cols): 3 (3, 1)
size (rows, cols): 9 (3, 3)

Пример DenseBase

Выводит подблок плотного выражения. Принимает любое плотное матричное или массивно выражение, но не разреженные объекты и не специальные классы матриц, такие как DiagonalMatrix.

template <typename Derived>
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
{
  std::cout << "block: " << b.block(x,y,r,c) << std::endl;
}

Пример ArrayBase

Выводит максимальный коэффициент массива или массивно-выражения.

template <typename Derived>
void print_max_coeff(const ArrayBase<Derived> &a)
{
  std::cout << "max: " << a.maxCoeff() << std::endl;
}

Пример MatrixBase

Выводит обратную условную число матрицы или матрично-выражения.

template <typename Derived>
void print_inv_cond(const MatrixBase<Derived>& a)
{
  const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType&
    sing_vals = a.jacobiSvd().singularValues();
  std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
}

Пример с несколькими шаблонизированными аргументами

Вычисляет евклидово расстояние между двумя точками.

template <typename DerivedA,typename DerivedB>
typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
{
  return (p1-p2).squaredNorm();
}

Обратите внимание, что мы использовали два параметра шаблона, по одному на аргумент. Это позволяет функции обрабатывать входные данные разных типов, например,

squaredist(v1,2*v2)

где первый аргумент v1 является вектором, а второй аргумент 2*v2 является выражением.

Эти примеры предназначены только для того, чтобы дать читателю первое представление о том, как можно создавать функции, которые принимают обычный и постоянный Matrix или Array аргумент. Они также предназначены для того, чтобы дать читателю представление о наиболее распространенных базовых классах, являющихся оптимальными кандидатами для функций. В следующем разделе мы более подробно рассмотрим пример и различные способы его реализации, обсуждая проблемы и преимущества каждой реализации. Для дальнейшего обсуждения Matrix и Array, а также MatrixBase и ArrayBase могут быть взаимозаменяемы, и все аргументы по-прежнему будут действительными.

Как написать общие, но не шаблонизированные функции?

Во всех предыдущих примерах функции были шаблонными функциями. Этот подход позволяет писать очень общий код, но часто желательно писать нешаблонные функции и сохранять определенный уровень общности, чтобы избежать бессмысленных копий аргументов. Типичным примером является создание функций, принимающих как MatrixXf, так и блок MatrixXf. Именно для этого служит класс Ref. Вот простой пример:

#include <iostream>
#include <Eigen/SVD>
using namespace Eigen;
using namespace std;
 
float inv_cond(const Ref<const MatrixXf>& a)
{
  const VectorXf sing_vals = a.jacobiSvd().singularValues();
  return sing_vals(sing_vals.size()-1) / sing_vals(0);
}
 
int main()
{
  Matrix4f m = Matrix4f::Random();
  cout << "matrix m:" << endl << m << endl << endl;
  cout << "inv_cond(m):          " << inv_cond(m)                      << endl;
  cout << "inv_cond(m(1:3,1:3)): " << inv_cond(m.topLeftCorner(3,3))   << endl;
  cout << "inv_cond(m+I):        " << inv_cond(m+Matrix4f::Identity()) << endl;
}

matrix m:
   0.68   0.823  -0.444   -0.27
 -0.211  -0.605   0.108  0.0268
  0.566   -0.33 -0.0452   0.904
  0.597   0.536   0.258   0.832

inv_cond(m):          0.0562343
inv_cond(m(1:3,1:3)): 0.0836819
inv_cond(m+I):        0.160204

В двух первых вызовах inv_cond копирование не происходит, потому что расположение памяти аргументов соответствует расположению памяти, принимаемому Ref<MatrixXf>. Однако в последнем вызове у нас есть общее выражение, которое автоматически будет вычислено в временную MatrixXf объектом Ref<>.

Объект Ref также может быть изменяемым. Вот пример функции, вычисляющей ковариационную матрицу двух входных матриц, где каждая строка представляет собой наблюдение:

void cov(const Ref<const MatrixXf> x, const Ref<const MatrixXf> y, Ref<MatrixXf> C)
{
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

и вот два примера вызова cov без копирования:

MatrixXf m1, m2, m3
cov(m1, m2, m3);
cov(m1.leftCols<3>(), m2.leftCols<3>(), m3.topLeftCorner<3,3>());

Класс Ref<> имеет два других необязательных параметра шаблона, позволяющих управлять типом расположения памяти, который может быть принят без копирования. См. документацию класса Ref для получения подробной информации.

В каких случаях функции, принимающие обычные аргументы Matrix или Array, работают?

Без использования шаблонных функций и без класса Ref, примитивное реализация предыдущей функции cov может выглядеть так:

MatrixXf cov(const MatrixXf& x, const MatrixXf& y)
{
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  return (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

и, вопреки тому, что можно было бы подумать на первый взгляд, эта реализация хороша, если только вы не требуете универсальной реализации, которая также работает с двойными матрицами, и если вы не беспокоитесь о временных объектах. Почему это так? Где участвуют временные объекты? Как может компилироваться предоставленный ниже код?

MatrixXf x,y,z;
MatrixXf C = cov(x,y+z);

В этом специальном случае пример работает, потому что оба параметра объявлены как const ссылки. Компилятор создает временную переменную и вычисляет выражение x+z в этой временной переменной. После обработки функции временная переменная освобождается, и результат присваивается C.

Примечание: Функции, принимающие const ссылки на Matrix (или Array), могут обрабатывать выражения за счет временных объектов.

В каких случаях функции, принимающие обычные аргументы Matrix или Array, не работают?

Здесь мы рассмотрим немного измененную версию функции, приведенной выше. На этот раз мы не хотим возвращать результат, а передаем дополнительный не-const параметр, который позволяет нам сохранить результат. Первая примитивная реализация может выглядеть следующим образом.

// Note: This code is flawed!
void cov(const MatrixXf& x, const MatrixXf& y, MatrixXf& C)
{
  const float num_observations = static_cast<float>(x.rows());
  const RowVectorXf x_mean = x.colwise().sum() / num_observations;
  const RowVectorXf y_mean = y.colwise().sum() / num_observations;
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

При попытке выполнить следующий код

MatrixXf C = MatrixXf::Zero(3,6);
cov(x,y, C.block(0,0,3,3));

компилятор выдаст ошибку, потому что невозможно преобразовать выражение, возвращаемое MatrixXf::block(), в не-const MatrixXf&. Это происходит потому, что компилятор хочет защитить вас от записи результата во временный объект. В этом конкретном случае эта защита не нужна – мы хотим записывать во временный объект. Итак, как можно решить эту проблему?

Предпочтительное решение на данный момент основано на небольшом обмане. Необходимо передать ссылку const на матрицу, а внутри const-ограничение должно быть снято. Правильная реализация для совместимых с C98 компиляторов будет

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;
 
  const Scalar num_observations = static_cast<Scalar>(x.rows());
 
  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;
 
  const_cast< MatrixBase<OtherDerived>& >(C) =
    (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

Теперь реализация не только работает с временными выражениями, но также позволяет использовать функцию с матрицами произвольных типов скалярных чисел с плавающей запятой.

Примечание: Хак с приведением const будет работать только с шаблонными функциями. Он не будет работать с реализацией MatrixXf, потому что невозможно привести выражение Block к ссылке на Matrix!

Как изменять размер матриц в универсальных реализациях?

Можно подумать, что мы закончили. Это не совсем так, потому что для того, чтобы наша функция ковариации была применимой в общем случае, мы хотим, чтобы следующий код работал:

MatrixXf x = MatrixXf::Random(100,3);
MatrixXf y = MatrixXf::Random(100,3);
MatrixXf C;
cov(x, y, C);

Это больше не так, когда мы используем реализацию, принимающую MatrixBase в качестве параметра. В общем случае Eigen поддерживает автоматическое изменение размера, но сделать это для выражений невозможно. Почему должно быть разрешено изменение размера блока Block матрицы? Это ссылка на подматрицу, и мы определенно не хотим изменять ее размер. Итак, как мы можем включить изменение размера, если мы не можем изменить размер MatrixBase? Решение заключается в изменении размера производного объекта, как в этой реализации.

template <typename Derived, typename OtherDerived>
void cov(const MatrixBase<Derived>& x, const MatrixBase<Derived>& y, MatrixBase<OtherDerived> const & C_)
{
  typedef typename Derived::Scalar Scalar;
  typedef typename internal::plain_row_type<Derived>::type RowVectorType;
 
  const Scalar num_observations = static_cast<Scalar>(x.rows());
 
  const RowVectorType x_mean = x.colwise().sum() / num_observations;
  const RowVectorType y_mean = y.colwise().sum() / num_observations;
 
  MatrixBase<OtherDerived>& C = const_cast< MatrixBase<OtherDerived>& >(C_);
  
  C.derived().resize(x.cols(),x.cols()); // resize the derived object
  C = (x.rowwise() - x_mean).transpose() * (y.rowwise() - y_mean) / num_observations;
}

Эта реализация теперь работает как для параметров, являющихся выражениями, так и для параметров, являющихся матрицами неверного размера. Изменение размера выражений в этом случае не навредит, если только они фактически не требуют изменения размера. Это означает, что передача выражения с неправильными размерами приведет к ошибке во время выполнения (только в режиме отладки), а передача выражения с правильными размерами просто сработает.

Примечание: В вышеприведенном обсуждении термины Matrix и Array и MatrixBase и ArrayBase можно использовать взаимозаменяемо, и все аргументы останутся действительными.

Обзор

• Подводя итог, реализация функций, принимающих не-записываемые (ссылочные const) объекты, не является большой проблемой и не приводит к проблемам при компиляции и запуске вашей программы. Однако, наивная реализация, вероятно, добавит ненужные временные объекты в ваш код. Чтобы избежать создания временных объектов из параметров, передавайте их как (const) ссылки на MatrixBase или ArrayBase (т.е. используйте шаблонную функцию).
• Функции, принимающие изменяемые (не-const) параметры, должны принимать ссылки const и удалять const в теле функции.
• Функции, которые принимают в качестве параметров объекты MatrixBase (или ArrayBase) и потенциально нуждаются в изменении их размера (в случае, если они изменяемые), должны вызывать resize() на производном классе, как возвращает derived().