Eigen学习笔记3:矩阵和向量运算

介绍

Eigen通过对常见的C ++算术运算符(+,-等)的重载或提供了诸如dot()、cross()等特殊方法以实现矩阵和向量算术运算。对于Matrix类(矩阵和向量),对运算符进行重载以支持线性代数运算。 For example, matrix1 * matrix2 means matrix-matrix product, and vector + scalar is just not allowed. If you want to perform all kinds of array operations, not linear algebra, see the next page.

加减

左侧和右侧必须具有相同数量的行和列,且具有相同的Scalar类型(因为Eigen不会自动进行类型升级)。

  • binary operator + as in a+b
  • binary operator - as in a-b
  • unary operator - as in -a
  • compound operator += as in a+=b
  • compound operator -= as in 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

标量乘法和除法

标量的乘法和除法也非常简单。这里的操作是:

  • binary operator * as in matrix*scalar
  • binary operator * as in scalar*matrix
  • binary operator / as in matrix/scalar
  • compound operator *= as in matrix*=scalar
  • compound operator /= as in 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

A note about expression templates

在Eigen中,算术运算符(例如:operator+)自身并不执行任何计算,它们仅返回描述要执行的计算的“表达式对象”。当计算整个表达式时,实际的计算将在稍后进行,通常在operator=中。

例如,当您这样做时:

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;
Eigen将其编译为一个for循环,因此数组仅被遍历一次。
for(int i = 0; i < 50; ++i)
  a[i] = 3*b[i] + 4*c[i] + 5*d[i];

因此,不必担心Eigen使用相对较大的算术表达式;它只会为Eigen提供更多优化机会。

转置和共轭矩阵

矩阵或向量的转置,共轭矩阵或者伴随矩阵(即共轭转置)分别可以通过成员函数transpose()conjugate()adjoint()获得。

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;


int main()
{
    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.211234,0.680375) (-0.604897,0.823295)
  (0.59688,0.566198) (0.536459,-0.329554)
Here is the matrix a^T
(-0.211234,0.680375)   (0.59688,0.566198)
(-0.604897,0.823295) (0.536459,-0.329554)
Here is the conjugate of a
(-0.211234,-0.680375) (-0.604897,-0.823295)
  (0.59688,-0.566198)   (0.536459,0.329554)
Here is the matrix a^*
(-0.211234,-0.680375)   (0.59688,-0.566198)
(-0.604897,-0.823295)   (0.536459,0.329554)

对于实数矩阵, conjugate() is a no-operation, and so adjoint() is equivalent to transpose()

transpose()adjoint()简单地返回一个对象但没有做实际的转换。b = a.transpose()then the transpose is evaluated at the same time as the result is written into b。但是,这里有一个复杂的问题。如果令a = a.transpose(),则starts writing the result into a before the evaluation of the transpose is finished。因此, the instruction a = a.transpose() does not replace a with its transpose, as one would expect:

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;


int main()
{
    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;

}

输出:

21:12:03: Starting /home/junjun/workspace/testEigen/build-testEigen-Imported_Kit-Default/testEigen ...
Here is the matrix a:
1 2
3 4
testEigen: /usr/include/eigen3/Eigen/src/Core/Transpose.h:378: static void Eigen::internal::checkTransposeAliasing_impl<Derived, OtherDerived, MightHaveTransposeAliasing>::run(const Derived&, const OtherDerived&) [with Derived = Eigen::Matrix<int, 2, 2>; OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2> >; bool MightHaveTransposeAliasing = true]: Assertion `(!check_transpose_aliasing_run_time_selector <typename Derived::Scalar,blas_traits<Derived>::IsTransposed,OtherDerived> ::run(extract_data(dst), other)) && "aliasing detected during transposition, use transposeInPlace() " "or evaluate the rhs into a temporary using .eval()"' failed.
21:12:03: 程序异常结束。
21:12:03: The process was ended forcefully.
21:12:03: /home/junjun/workspace/testEigen/build-testEigen-Imported_Kit-Default/testEigen crashed.

This is the so-called aliasing issue. 在“调试模式”下,未禁用断言时,会自动检测到此类常见的陷阱。

For in-place transposition, as for instance in a = a.transpose(), simply use the transposeInPlace() function:

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;
using namespace std;


int main()
{
    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()函数。

Matrix-matrix and matrix-vector乘法

使用矩阵矩阵乘法operator*。由于向量是矩阵的一种特殊情况,因此它们也被隐式地处理,所以矩阵向量乘积实际上只是矩阵-矩阵乘积的一种特殊情况,向量-向量外积也是如此。因此,所有这些情况仅由两个运算符处理:

  • binary operator * as in a*b
  • compound operator *= as in a*=b (this multiplies on the right: a*=b is equivalent to 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;
}

注意:如果担心这样做m=m*m可能会导致混叠问题,请立即放心:Eigen将矩阵乘法视为一种特殊情况,并在此处引入了一个临时值,因此它将编译m=m*m为:

tmp = m * m;
m = tmp;

If you know your matrix product can be safely evaluated into the destination matrix without aliasing issue, then you can use the noalias() function to avoid the temporary, e.g.:

c.noalias()+ = a * b;
or more details on this topic, see the page on  aliasing.

注意: for BLAS users worried about performance, expressions such as c.noalias() -= 2 * a.adjoint() * b; are fully optimized and trigger a single gemm-like function call.

点积和叉积

对于点积和叉积,您需要dot()cross()方法。当然,点积也可以作为u.adjoint()* v的1x1矩阵获得。

#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的向量。点积适用于任何大小的向量。When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the second variable.

Basic arithmetic reduction operations

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

The trace of a matrix, as returned by the function trace(), is the sum of the diagonal coefficients and can also be computed as efficiently using a.diagonal().sum(), as we will see later on.

还存在minCoeffmaxCoeff函数的变体,它们通过参数返回相应系数的坐标:

#include <iostream>
#include <Eigen/Dense>

using namespace std;
using namespace Eigen;

int main()
{
    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.680375   0.59688 -0.329554
-0.211234  0.823295  0.536459
 0.566198 -0.604897 -0.444451
Its minimum coefficient (-0.604897) is at position (2,1)

Here is the vector v:  115899597  -48539462  276748203 -290373134
Its maximum coefficient (276748203) is at position 2

操作的有效性

Eigen检查您执行的操作的有效性。如果可能,它将在编译时检查它们,从而产生编译错误。这些错误消息可能很长很丑,但是Eigen将重要消息写在UPPERCASE_LETTERS_SO_IT_STANDS_OUT中。例如:

Matrix3f m;
Vector4f v;
v = m * v; //编译时错误:YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
当然,在许多情况下,例如在检查动态大小时,无法在编译时执行检查。然后,Eigen使用运行时断言。这意味着,如果程序在“调试模式”下运行,则在执行非法操作时该程序将中止并显示一条错误消息,并且如果断言被关闭,它可能会崩溃。
MatrixXf m(3,3);
VectorXf v(4);
v = m * v; //此处的运行时断言失败:“无效的矩阵乘积”
For more details on this topic, see  this page.
已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页