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

Введение

Eigen предлагает операции матричной/векторной арифметики, либо через перегрузки общих операторов арифметики C++, таких как +, -, *, либо через специальные методы, такие как dot(), cross() и т.д. Для класса Matrix (матрицы и векторы), операторы перегружены только для поддержки операций линейной алгебры. Например, matrix1 * matrix2 означает произведение матриц, а vector + scalar просто не разрешено. Если вы хотите выполнить все виды операций с массивами, не связанные с линейной алгеброй, обратитесь к следующей странице.

Сложение и вычитание

Левая и правая части, конечно, должны иметь одинаковое количество строк и столбцов. Они также должны иметь одинаковый Scalar тип, так как Eigen не выполняет автоматическое повышение типа. Операторы, имеющиеся в распоряжении:

• бинарный оператор + как в a+b
• бинарный оператор - как в a-b
• унарный оператор - как в -a
• составной оператор += как в a+=b
• составной оператор -= как в a-=b

#include <iostream>
#include <Eigen/Dense>
 
using namespace Eigen;
 
int main()
{
  Matrix2d a;
  a << 1, 2,
       3, 4;
  MatrixXd b(2,2);
  b << 2, 3,
       1, 4;
  std::cout << "a + b =\n" << a + b << std::endl;
  std::cout << "a - b =\n" << a - b << std::endl;
  std::cout << "Doing a += b;" << std::endl;
  a += b;
  std::cout << "Now a =\n" << a << std::endl;
  Vector3d v(1,2,3);
  Vector3d w(1,0,0);
  std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}

a + b =
3 5
4 8
a - b =
-1 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

Умножение и деление на скаляр

Умножение и деление на скаляр также очень просто. Операторы, имеющиеся в распоряжении:

• бинарный оператор * как в matrix*scalar
• бинарный оператор * как в scalar*matrix
• бинарный оператор / как в matrix/scalar
• составной оператор *= как в matrix*=scalar
• составной оператор /= как в matrix/=scalar

#include <iostream>
#include <Eigen/Dense>
 
using namespace Eigen;
 
int main()
{
  Matrix2d a;
  a << 1, 2,
       3, 4;
  Vector3d v(1,2,3);
  std::cout << "a * 2.5 =\n" << a * 2.5 << std::endl;
  std::cout << "0.1 * v =\n" << 0.1 * v << std::endl;
  std::cout << "Doing v *= 2;" << std::endl;
  v *= 2;
  std::cout << "Now v =\n" << v << std::endl;
}

a * 2.5 =
2.5   5
7.5  10
0.1 * v =
0.1
0.2
0.3
Doing v *= 2;
Now v =
2
4
6

Примечание об шаблонах выражений

Это продвинутая тема, которую мы объясняем на этой странице, но полезно упомянуть её сейчас. В Eigen, арифметические операторы, такие как operator+ не выполняют вычислений сами по себе, они просто возвращают «объект выражения», описывающий вычисления, которые необходимо выполнить. Фактические вычисления выполняются позже, когда всё выражение вычисляется, как правило, в operator=. Хотя это может показаться сложным, любой современный оптимизирующий компилятор способен оптимизировать это абстрагирование, и результатом является идеально оптимизированный код. Например, когда вы делаете:

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;

Eigen компилирует его всего лишь в один цикл, так что массивы обрабатываются только один раз. Упрощая (например, игнорируя оптимизации SIMD), этот цикл выглядит так:

for(int i = 0; i < 50; ++i)
  a[i] = 3*b[i] + 4*c[i] + 5*d[i];

Таким образом, вы не должны бояться использования относительно больших арифметических выражений с Eigen: это просто предоставляет Eigen больше возможностей для оптимизации.

Транспонирование и сопряжение

Транспонированная \( a^T \), сопряжённая \( \bar{a} \), и сопряжённо-транспонированная (т.е. адъюнкт) \( a^* \) матрицы или вектора \( a \) получаются с помощью функций-членов transpose(), conjugate() и adjoint() соответственно.

MatrixXcf a = MatrixXcf::Random(2,2);
cout << "Here is the matrix a\n" << a << endl;
 
cout << "Here is the matrix a^T\n" << a.transpose() << endl;
 
 
cout << "Here is the conjugate of a\n" << a.conjugate() << endl;
 
 
cout << "Here is the matrix a^*\n" << a.adjoint() << endl;
 
 

Here is the matrix a
 (-0.211,0.68) (-0.605,0.823)
 (0.597,0.566)  (0.536,-0.33)
Here is the matrix a^T
 (-0.211,0.68)  (0.597,0.566)
(-0.605,0.823)  (0.536,-0.33)
Here is the conjugate of a
 (-0.211,-0.68) (-0.605,-0.823)
 (0.597,-0.566)    (0.536,0.33)
Here is the matrix a^*
 (-0.211,-0.68)  (0.597,-0.566)
(-0.605,-0.823)    (0.536,0.33)

Для вещественных матриц conjugate() — это бесполезная операция, и поэтому adjoint() эквивалентно transpose().

Что касается основных арифметических операторов, transpose() и adjoint() просто возвращают прокси-объект без выполнения фактического транспонирования. Если вы сделаете b = a.transpose(), то транспонирование будет вычислено одновременно с записью результата в b. Однако здесь есть проблема. Если вы сделаете a = a.transpose(), то Eigen начнёт запись результата в a до завершения вычисления транспонирования. Поэтому инструкция a = a.transpose() не заменяет a своей транспонированной, как можно было бы ожидать:

Matrix2i a; a << 1, 2, 3, 4;
cout << "Here is the matrix a:\n" << a << endl;
 
a = a.transpose(); // !!! do NOT do this !!!
cout << "and the result of the aliasing effect:\n" << a << endl;

Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4

Это так называемая проблема алиасинга. В «отладочном режиме», т.е. когда утверждения не были отключены, подобные общие подводные камни автоматически обнаруживаются.

Для непосредственного транспонирования, например, в a = a.transpose(), просто используйте функцию transposeInPlace():

MatrixXf a(2,3); a << 1, 2, 3, 4, 5, 6;
cout << "Here is the initial matrix a:\n" << a << endl;
 
 
a.transposeInPlace();
cout << "and after being transposed:\n" << a << endl;

Here is the initial matrix a:
1 2 3
4 5 6
and after being transposed:
1 4
2 5
3 6

Также существует функция adjointInPlace() для комплексных матриц.

Умножение матриц и матриц-векторов

Умножение матриц выполняется с помощью operator*. Поскольку векторы являются частным случаем матриц, они также обрабатываются неявно, поэтому произведение матрицы-вектора — это всего лишь частный случай произведения матриц, а также внешнее произведение вектор-вектор. Таким образом, все эти случаи обрабатываются всего двумя операторами:

• бинарный оператор * как в a*b
• составной оператор *= как в a*=b (это умножение справа: a*=b эквивалентно a = a*b)

#include <iostream>
#include <Eigen/Dense>
 
using namespace Eigen;
int main()
{
  Matrix2d mat;
  mat << 1, 2,
         3, 4;
  Vector2d u(-1,1), v(2,0);
  std::cout << "Here is mat*mat:\n" << mat*mat << std::endl;
  std::cout << "Here is mat*u:\n" << mat*u << std::endl;
  std::cout << "Here is u^T*mat:\n" << u.transpose()*mat << std::endl;
  std::cout << "Here is u^T*v:\n" << u.transpose()*v << std::endl;
  std::cout << "Here is u*v^T:\n" << u*v.transpose() << std::endl;
  std::cout << "Let's multiply mat by itself" << std::endl;
  mat = mat*mat;
  std::cout << "Now mat is mat:\n" << mat << std::endl;
}

Here is mat*mat:
 7 10
15 22
Here is mat*u:
1
1
Here is u^T*mat:
2 2
Here is u^T*v:
-2
Here is u*v^T:
-2 -0
 2  0
Let's multiply mat by itself
Now mat is mat:
 7 10
15 22

Примечание: если вы прочитали вышеупомянутый абзац о шаблонах выражений и беспокоитесь о том, что выполнение m=m*m может вызвать проблемы с алиасингом, успокойтесь пока: Eigen обрабатывает умножение матриц как частный случай и позаботится о создании временного объекта, поэтому он будет компилировать m=m*m как:

tmp = m*m;
m = tmp;

Если вы знаете, что ваше произведение матриц может быть безопасно вычислено в целевой матрице без проблем с алиасингом, тогда вы можете использовать функцию noalias(), чтобы избежать временного объекта, например:

c.noalias() += a * b;

Для получения дополнительной информации по этой теме см. страницу алиасинг.

Примечание: для пользователей BLAS, которые беспокоятся о производительности, выражения, такие как c.noalias() -= 2 * a.adjoint() * b; полностью оптимизированы и вызывают один вызов функции типа gemm.

Скалярное произведение и векторное произведение

Для скалярного произведения и векторного произведения вам нужны методы dot() и cross(). Конечно, скалярное произведение также можно получить как 1x1 матрицу как u.adjoint()*v.

#include <iostream>
#include <Eigen/Dense>
 
using namespace Eigen;
using namespace std;
int main()
{
  Vector3d v(1,2,3);
  Vector3d w(0,1,2);
 
  cout << "Dot product: " << v.dot(w) << endl;
  double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
  cout << "Dot product via a matrix product: " << dp << endl;
  cout << "Cross product:\n" << v.cross(w) << endl;
}

Dot product: 8
Dot product via a matrix product: 8
Cross product:
 1
-2
 1

Помните, что векторное произведение используется только для векторов размера 3. Скалярное произведение используется для векторов любого размера. При использовании комплексных чисел скалярное произведение Eigen является сопряжённо-линейным относительно первой переменной и линейным относительно второй.

Основные арифметические операции редукции

Eigen также предоставляет некоторые операции редукции для редукции заданной матрицы или вектора до одного значения, например, сумма (вычисляется с помощью sum()), произведение (prod()), или максимальное (maxCoeff()) и минимальное (minCoeff()) из всех его коэффициентов.

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

Here is mat.sum():       10
Here is mat.prod():      24
Here is mat.mean():      2.5
Here is mat.minCoeff():  1
Here is mat.maxCoeff():  4
Here is mat.trace():     5

След матрицы, возвращаемый функцией trace(), — это сумма диагональных коэффициентов, а также может быть эффективно вычислена с помощью a.diagonal().sum(), как мы увидим позже.

Также существуют варианты функций minCoeff и maxCoeff возвращающих координаты соответствующего коэффициента через аргументы:

  Matrix3f m = Matrix3f::Random();
  std::ptrdiff_t i, j;
  float minOfM = m.minCoeff(&i,&j);
  cout << "Here is the matrix m:\n" << m << endl;
  cout << "Its minimum coefficient (" << minOfM 
       << ") is at position (" << i << "," << j << ")\n\n";
 
  RowVector4i v = RowVector4i::Random();
  int maxOfV = v.maxCoeff(&i);
  cout << "Here is the vector v: " << v << endl;
  cout << "Its maximum coefficient (" << maxOfV 
       << ") is at position " << i << endl;

Here is the matrix m:
  0.68  0.597  -0.33
-0.211  0.823  0.536
 0.566 -0.605 -0.444
Its minimum coefficient (-0.605) is at position (2,1)

Here is the vector v:  1  0  3 -3
Its maximum coefficient (3) is at position 2

Действительность операций

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

Matrix3f m;
Vector4f v;
v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES

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

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"

Для получения более подробной информации по этой теме, см. эту страницу.