В Eigen алиасинг относится к оператору присваивания, в котором та же матрица (или массив или вектор) появляется слева и справа от операторов присваивания. Такие операторы как mat = 2 * mat; или mat = mat.transpose(); демонстрируют алиасинг. Алиасинг в первом примере безобиден, но алиасинг во втором примере приводит к неожиданным результатам. Эта страница объясняет, что такое алиасинг, когда он вреден и что с этим делать.

Примеры

Вот простой пример, демонстрирующий алиасинг:

MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
 
// This assignment shows the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
cout << "After the assignment, mat = \n" << mat << endl;

Here is the matrix mat:
1 2 3
4 5 6
7 8 9
After the assignment, mat = 
1 2 3
4 1 2
7 4 1

Вывод не такой, как ожидалось. Проблема в присваивании

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);

Это присваивание демонстрирует алиасинг: коэффициент mat(1,1) появляется как в блоке mat.bottomRightCorner(2,2) в левой части оператора присваивания, так и в блоке mat.topLeftCorner(2,2) в правой части. После присваивания элемент (2,2) в правом нижнем углу должен иметь значение mat(1,1) до присваивания, которое равно 5. Однако вывод показывает, что mat(2,2) фактически равен 1. Проблема в том, что Eigen использует ленивое вычисление (см. Шаблоны выражений в Eigen) для mat.topLeftCorner(2,2). Результат аналогичен

mat(1,1) = mat(0,0);
mat(1,2) = mat(0,1);
mat(2,1) = mat(1,0);
mat(2,2) = mat(1,1);

Таким образом, mat(2,2) получает новое значение mat(1,1) вместо старого значения. Следующий раздел объясняет, как решить эту проблему, вызвав eval().

Алиасинг возникает более естественно при попытке уменьшить матрицу. Например, выражения vec = vec.head(n) и mat = mat.block(i,j,r,c) демонстрируют алиасинг.

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

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

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

void Eigen::DenseBase<Derived>::checkTransposeAliasing(const OtherDerived&) const 
[with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]: 
Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other)) 
&& "aliasing detected during transposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed.

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

Решение проблем с алиасингом

Если вы понимаете причину проблемы с алиасингом, то очевидно, что необходимо сделать, чтобы ее решить: Eigen должен полностью вычислить правую часть в временной матрице/массиве, а затем присвоить ее левой части. Функция eval() делает именно это.

Например, вот исправленная версия первого примера выше:

MatrixXi mat(3,3); 
mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
 
// The eval() solves the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
cout << "After the assignment, mat = \n" << mat << endl;

Here is the matrix mat:
1 2 3
4 5 6
7 8 9
After the assignment, mat = 
1 2 3
4 1 2
7 4 5

Теперь mat(2,2) равно 5 после присваивания, как и должно быть.

То же самое решение также работает для второго примера с транспонированием: просто замените строку a = a.transpose(); на a = a.transpose().eval();. Однако в этом общем случае есть лучшее решение. Eigen предоставляет функцию специального назначения 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

Если функция xxxInPlace() доступна, то лучше всего ее использовать, потому что это более ясно показывает, что вы делаете. Это также может позволить Eigen оптимизировать вычисления более агрессивно. Вот некоторые функции xxxInPlace():

В особом случае, когда матрица или вектор уменьшаются с помощью выражения, такого как vec = vec.head(n), можно использовать conservativeResize() .

Алиасинг и покомпонентные операции

Как объяснялось выше, опасно, если та же матрица или массив встречаются как в левой, так и в правой части оператора присваивания, и часто необходимо явно вычислить правую часть. Однако применение покомпонентных операций (таких как сложение матриц, умножение на скаляр и умножение массивов) безопасно.

В следующем примере используются только покомпонентные операции. Поэтому eval() не требуется, даже если та же матрица появляется с обеих сторон операторов присваивания.

MatrixXf mat(2,2); 
mat << 1, 2,  4, 7;
cout << "Here is the matrix mat:\n" << mat << endl << endl;
 
mat = 2 * mat;
cout << "After 'mat = 2 * mat', mat = \n" << mat << endl << endl;
 
 
mat = mat - MatrixXf::Identity(2,2);
cout << "After the subtraction, it becomes\n" << mat << endl << endl;
 
 
ArrayXXf arr = mat;
arr = arr.square();
cout << "After squaring, it becomes\n" << arr << endl << endl;
 
// Combining all operations in one statement:
mat << 1, 2,  4, 7;
mat = (2 * mat - MatrixXf::Identity(2,2)).array().square();
cout << "Doing everything at once yields\n" << mat << endl << endl;

Here is the matrix mat:
1 2
4 7

After 'mat = 2 * mat', mat = 
 2  4
 8 14

After the subtraction, it becomes
 1  4
 8 13

After squaring, it becomes
  1  16
 64 169

Doing everything at once yields
  1  16
 64 169

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

Алиасинг и умножение матриц

Matrix умножение — единственная операция в Eigen, которая по умолчанию предполагает алиасинг при условии, что размерность целевой матрицы не изменяется. Таким образом, если matA — квадратная матрица, то оператор matA = matA * matA; безопасен. Все остальные операции в Eigen предполагают, что проблем с алиасингом нет, либо потому, что результат присваивается другой матрице, либо потому, что это покомпонентная операция.

MatrixXf matA(2,2); 
matA << 2, 0,  0, 2;
matA = matA * matA;
cout << matA;

4 0
0 4

Однако это имеет свои недостатки. При выполнении выражения matA = matA * matA, Eigen вычисляет произведение во временной матрице, которая присваивается matA после вычисления. Это нормально. Но Eigen делает то же самое, когда произведение присваивается другой матрице (например, matB = matA * matA). В этом случае более эффективно вычислить произведение непосредственно в matB вместо того, чтобы сначала вычислить его во временной матрице и скопировать эту матрицу в matB.

Пользователь может указать с помощью функции noalias(), что алиасинг отсутствует, следующим образом: matB.noalias() = matA * matA. Это позволяет Eigen вычислить матричное произведение matA * matA непосредственно в matB.

MatrixXf matA(2,2), matB(2,2); 
matA << 2, 0,  0, 2;
 
// Simple but not quite as efficient
matB = matA * matA;
cout << matB << endl << endl;
 
// More complicated but also more efficient
matB.noalias() = matA * matA;
cout << matB;

4 0
0 4

4 0
0 4

Конечно, вы не должны использовать noalias() при наличии алиасинга. В противном случае вы можете получить неправильные результаты:

MatrixXf matA(2,2); 
matA << 2, 0,  0, 2;
matA.noalias() = matA * matA;
cout << matA;

4 0
0 4

Кроме того, начиная с Eigen 3.3, алиасинг не предполагается, если размерность целевой матрицы изменяется, и произведение не присваивается непосредственно цели. Поэтому следующий пример также неверен:

MatrixXf A(2,2), B(3,2);
B << 2, 0,  0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).cwiseAbs();
cout << A;

4 0
0 6
2 2

Как и в любой проблеме с алиасингом, ее можно решить, явно вычислив выражение перед присваиванием:

MatrixXf A(2,2), B(3,2);
B << 2, 0,  0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).eval().cwiseAbs();
cout << A;

4 0
0 6
2 2

Сводка

Алиасинг возникает, когда одни и те же коэффициенты матрицы или массива появляются в левой и правой частях оператора присваивания.

• Алиасинг безопасен при покомпонентных вычислениях; это включает умножение на скаляр и сложение матриц или массивов.
• При умножении двух матриц Eigen по умолчанию предполагает, что происходит алиасинг. Если вы знаете, что алиасинга нет, то вы можете использовать noalias().
• Во всех остальных случаях Eigen предполагает, что проблемы с алиасингом нет, и поэтому дает неправильный результат, если алиасинг действительно имеет место. Чтобы предотвратить это, вам нужно использовать eval() или одну из функций xxxInPlace().