Eigen中的Aliasing问题

Eigen中的Aliasing问题

在看贺一家博士的代码中发现,在矩阵的计算中他写入了一个.noalias(),觉得比较奇怪,所以就去搜索一下eigen中alias的含义

文章主要参考eigen官方的文档

Aliasing产生的原因与问题

Aliasing问题产生最直接的情况就是在赋值的左边和右边同时存在同一个矩阵,比如以下的例子

1
2
3
4
5
6
7
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;

其中mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);目的是让左上角的部分赋值给同一个矩阵的右下角.这段代码的输出是:

1
2
3
4
5
6
7
8
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

可以发现,在坐标为(2 , 2)的位置,我们希望是5,但是实际上是1.

在执行上述的操作时,实际上Eigen只是在对同一个矩阵在做运算,如下:

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

但我们实际上需要的是将原始数据存在一个temp中,然后再从temp中拿出需要的对原始矩阵做一个填补

这个问题同样出现在矩阵的缩小以及矩阵转置中

如下代码所示

1
2
3
4
5
6
7
8
9
10
// shrink problem 
mat = mat.block(0,0,2,2);
cout << "After shrink, mat = \n" << mat << endl;

// transpose problem
Eigen::MatrixXi a(3,3);
a << 1,2,3,4,5,6,7,8,9;
cout <<"Here is the matrix a:\n" << a << endl;
a = a.transpose();
cout << "After transpose matrtix a: \n" << a << endl;

实际的输出结果为

1
2
3
4
5
6
7
8
9
10
11
After shrink, mat = 
12805328 2
0 1
Here is the matrix a:
1 2 3
4 5 6
7 8 9
After transpose matrtix a:
1 2 3
2 5 6
3 6 9

这些问题和之前所述是一致的

解决Aliasing问题

Eigen有自己的方式去使用一个temp来保证原矩阵的完整性,也就是使用.eval()来实现

1
2
3
4
5
6
7
8
9
10
// 使用eval,解决赋值问题
Eigen::MatrixXi mat2(3,3);
mat2 << 1,2,3,4,5,6,7,8,9;
cout <<"Here is the matrix mat2:\n" << mat2 << endl;

mat2.bottomRightCorner(2,2) = mat2.topLeftCorner(2,2).eval();
cout << "After the assignment, mat = \n" << mat2 << endl;

mat2 = mat2.block(0,0,2,2).eval();
cout << "After shrink, mat = \n" << mat2 << endl;

在变换之后的矩阵后添加eval,可以保证原矩阵的完整性.输出结果为

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

和我们希望的答案是一致的

同样,使用eval可以解决transpose的问题.但是eigen提供了一种更好的方式,也就是xxxInPlace()这种函数

1
2
3
4
5
6
7
8
// 使用transposeinplace
Eigen::MatrixXi b(3,3);
b << 1,2,3,4,5,6,7,8,9;
cout <<"Here is the matrix b:\n" << b << endl;
// transpose problem
// b = b.transpose().eval()
b.transposeInPlace();
cout << "After transpose matrtix b: \n" << b << endl;

输出结果

1
2
3
4
5
6
7
8
Here is the matrix b:
1 2 3
4 5 6
7 8 9
After transpose matrtix b:
1 4 7
2 5 8
3 6 9