zoukankan      html  css  js  c++  java
  • Eigen中的noalias(): 解决矩阵运算的混淆问题

    作者:@houkai
    本文为作者原创,转载请注明出处:http://www.cnblogs.com/houkai/p/6349990.html


    目录

    混淆
    例子
    解决混淆问题
    混淆和component级的操作。
    混淆和矩阵的乘法
    总结

    整理下Eigen库的教程,参考:http://eigen.tuxfamily.org/dox/index.html

    混淆

    在Eigen中,当变量同时出现在左值和右值,赋值操作可能会带来混淆问题。这一篇将解释什么是混淆,什么时候是有害的,怎么使用做。

    例子

    MatrixXi mat(3,3); 
    mat << 1, 2, 3,   4, 5, 6,   7, 8, 9;
    cout << "Here is the matrix mat:
    " << mat << endl;
    // This assignment shows the aliasing problem
    mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
    cout << "After the assignment, mat = 
    " << 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) 在bottomRightCorner(2,2)和topLeftCorner(2,2)都存在。赋值结果中mat(2,2)本应该赋予操作前mat(1,1)的值=5。但是,最终程序结果mat(2,2)=1。原因是Eigen使用了lazy evaluation(懒惰评估),上面等价于

    mat(1,1) = mat(0,0);
    mat(1,2) = mat(0,1);
    mat(2,1) = mat(1,0);
    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:
    " << a << endl;
    a = a.transpose(); // !!! do NOT do this !!!
    cout << "and the result of the aliasing effect:
    " << a << endl;
    
    Here is the matrix a:
    1 2
    3 4
    and the result of the aliasing effect:
    1 2
    2 4
    

    我们可以通过EIGEN_NO_DEBUG宏,在编译时关闭运行时的断言。

    解决混淆问题

    Eigen需要把右值赋值为一个临时matrix/array,然后再将临时值赋值给左值,便可以解决混淆。eval()函数实现了这个功能。

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

    同样: a = a.transpose().eval(); ,当然我们最好使用 transposeInPlace()。如果存在xxxInPlace函数,推荐使用这类函数,它们更加清晰地标明了你在做什么。提供的这类函数:

    OriginIn-place
    MatrixBase::adjoint() MatrixBase::adjointInPlace()
    DenseBase::reverse() DenseBase::reverseInPlace()
    LDLT::solve() LDLT::solveInPlace()
    LLT::solve() LLT::solveInPlace()
    TriangularView::solve() TriangularView::solveInPlace()
    DenseBase::transpose() DenseBase::transposeInPlace()

    而针对vec = vec.head(n)这种情况,推荐使用conservativeResize()

    混淆和component级的操作。

    组件级是指整体的操作,比如matrix加法、scalar乘、array乘等,这类操作是安全的,不会出现混淆。

    MatrixXf mat(2,2); 
    mat << 1, 2,  4, 7;
    cout << "Here is the matrix mat:
    " << mat << endl << endl;
    mat = 2 * mat;
    cout << "After 'mat = 2 * mat', mat = 
    " << mat << endl << endl;
    mat = mat - MatrixXf::Identity(2,2);
    cout << "After the subtraction, it becomes
    " << mat << endl << endl;
    ArrayXXf arr = mat;
    arr = arr.square();
    cout << "After squaring, it becomes
    " << arr << 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
    

    混淆和矩阵的乘法

    在Eigen中,矩阵的乘法一般都会出现混淆。除非是方阵(实质是元素级的乘)。

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

    其他的操作,Eigen默认都是存在混淆的。所以Eigen对矩阵乘法自动引入了临时变量,对的matA=matA*matA这是必须的,但是对matB=matA*matA这样便是不必要的了。我们可以使用noalias()函数来声明这里没有混淆,matA*matA的结果可以直接赋值为matB。

    matB.noalias() = matA * matA;
    

    从Eigen3.3开始,如果目标矩阵resize且结果不直接赋值给目标矩阵,默认不存在混淆。

    MatrixXf A(2,2), B(3,2);
    B << 2, 0,  0, 3, 1, 1;
    A << 2, 0, 0, -2;
    A = (B * A).cwiseAbs();//cwiseAbs()不直接赋给目标
    //A = (B * A).eval().cwiseAbs()
    cout << A;
    

    当然,对于任何混淆问题,都可以通过matA=(matB*matA).eval() 来解决。

    总结

    当相同的矩阵或array在等式左右都出现时,很容易出现混淆。

    1. compnent级别的操作不用考虑混淆。
    2. 矩阵相乘,Eigen默认会解决混淆问题,如果你确定不会出现混淆,可以使用noalias()来提效。
    3. 混淆出现时,可以用eval()和xxxInPlace()函数解决。
  • 相关阅读:
    Linq聚合操作之Aggregate,Count,Sum,Distinct源码分析
    Linq分区操作之Skip,SkipWhile,Take,TakeWhile源码分析
    Linq生成操作之DefautIfEmpty,Empty,Range,Repeat源码分析
    Linq基础操作之Select,Where,OrderBy,ThenBy源码分析
    PAT 1152 Google Recruitment
    PAT 1092 To Buy or Not to Buy
    PAT 1081 Rational Sum
    PAT 1084 Broken Keyboard
    PAT 1077 Kuchiguse
    PAT 1073 Scientific Notation
  • 原文地址:https://www.cnblogs.com/defe-learn/p/7456778.html
Copyright © 2011-2022 走看看