特色技术

专业从事预应力结构体系的设计、施工一体化解决方案

资讯分类

线性代数计算库Eigen(四)

  • 分类:程序技术
  • 作者:
  • 来源:
  • 发布时间:2018-11-08 15:32
  • 访问量:

【概要描述】(六)基础操作:  本小节旨在提供关于如何用Eigen执行矩阵、向量、标量运算的一个概述和一些细节。  1、概述  eigen通过重载常见的C+算术运算符(如+、-、*),或通过dot()、cross()等特殊方法提供矩阵、向量的算术运算。对于矩阵类(矩阵和向量),运算符仅被重载以支持线性代数运算。例如,matrix1*matrix2意味着矩阵乘积,而vector+scalar(向量+标量)是不允

线性代数计算库Eigen(四)

【概要描述】(六)基础操作:  本小节旨在提供关于如何用Eigen执行矩阵、向量、标量运算的一个概述和一些细节。  1、概述  eigen通过重载常见的C+算术运算符(如+、-、*),或通过dot()、cross()等特殊方法提供矩阵、向量的算术运算。对于矩阵类(矩阵和向量),运算符仅被重载以支持线性代数运算。例如,matrix1*matrix2意味着矩阵乘积,而vector+scalar(向量+标量)是不允

  • 分类:程序技术
  • 作者:
  • 来源:
  • 发布时间:2018-11-08 15:32
  • 访问量:
详情

  (六)基础操作:

  本小节旨在提供关于如何用Eigen执行矩阵、向量、标量运算的一个概述和一些细节。

  1、概述

  eigen通过重载常见的C+算术运算符(如+、-、*),或通过dot()、cross()等特殊方法提供矩阵、向量的算术运算。对于矩阵类(矩阵和向量),运算符仅被重载以支持线性代数运算。例如,matrix1*matrix2意味着矩阵乘积,而vector+scalar(向量+标量)是不允许的。

  2、加减运算:

  要进行矩阵、向量的加减运算,前提条件是左侧和右侧必须具有相同的行数和列数,它们还必须具有相同的标量类型,因为eigen不进行自动类型升级。这里的操作是:

  (1)双目运算符+:例如a+b

  (2) 双目运算符-:例如a-b

  (3)单目运算符-:例如-a

  (4)组合运算符+=:例如a+=b

  (5)组合运算符-=:例如a-=b

  例题:

  运算结果:

  说明:如果代码无法运行,提示找不到类库,可以按照下图对visual studio进行设置。

  3、标量乘除:

  被标量乘除也是非常简单的,它的操作是:

  (1)双目运算符*:例如matrix*scalar

  (2) 双目运算符*:例如scalar*matrix

  (3)双目运算符/:例如matrix/scalar

  (4)组合运算符*=:例如matrix*=scalar

  (5)组合运算符/=:例如matrix/=scalar

  例题:

  运算结果:

  4、关于表达式模板的一点说明

  这是我们在本节中需要解释的一个高级主题,现在讲解它是非常有用的。运算符(如operator+)本身并不执行任何计算,它们只是返回一个描述要执行计算的“表达式对象”(expression object)。 当计算整个表达式时,实际的计算发生在后面,通常是在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更多的优化机会。

  5、转置和共轭:

  (1)矩阵A的转置(transpose)矩阵AT,共轭(conjugate)矩阵Ᾱ,共轭转置(adjoint)矩阵A*,通过方法函数transpose(), conjugate(),adjoint()得到。

  (2)名词解释:

  共轭复数:实数部分相同而虚数部分互为相反数的两个复数。如1+2i的共轭复数是1-2i。

  矩阵的共轭转置(adjoint):把矩阵转置后,再把每一个数换成它的共轭复数。通常记做A*或者A^H。一般来讲A^H的写法不会有歧义。

  另外,A^*被称为伴随矩阵,用adj(A)表示A的伴随不会有歧义。伴随矩阵(adjugate matrix)是指由矩阵各元素的代数余子式组成的矩阵,与共轭转置矩阵(adjoint matrix)有本质区别。

  A的转置共轭A*(A^H)和A的伴随阵adj(A)没有直接关系。

  共轭矩阵究竟如何定义取决于其中的“共轭”是如何定义的:

  a.如果共轭是指“复共轭”,那么共轭矩阵就是指相应的复共轭矩阵,记为A*

  b.如果共轭是指“厄米共轭”,那么共轭矩阵就是指相应的复共轭转置矩阵(一些线性代数书也会将共轭转置矩阵记为A*或A^H),在量子力学的书里会将之记成A^+

  共轭是一个很广泛的概念,不同书的作者会对之有不同的定义,不同的记法,目前尚无统一的定义。在具体使用的时候只需对所说的“共轭”以及相应的记号作适当的说明即可。

  c.伴随矩阵有两种,这是由于不同的英文翻译成相同的中文所造成的:

  Adjugate matrix,即A的余子矩阵的转置矩阵

  Adjoint matrix,即A的复共轭转置矩阵,在这种翻译情况下,伴随矩阵与厄米共轭矩阵是指同一样东西。

  (3)例题:

  运行结果:

  (4)对实数矩阵,conjugate()函数不改变原矩阵,也就是无操作,adjoint()函数与transpose()函数等效。

  (5) 对于基本算术运算符,transspose()和adjoint()只是返回一个代理对象,而不执行实际的换位操作。如果你做了b=a.trspose(),那么转置将在结果写入b的同时进行计算。然而,这里有一个复杂的问题。如果您执行a=a.transpose(),那么Eigen在完成转置计算之前就开始将结果写入到a中。因此,指令a=a.transspose()并没有像人们预期的那样用转置来替换a。在最新版的Eigen中,已经禁止这样操作。

  为了实现自我转置,可以简单地应用b=b.transposeInPlace()函数。相应地,adjointInPlace()函数也可以被应用在复矩阵中。

  例如:

  运行结果:

  6、矩阵乘法:

  (1)矩阵乘法依然用operator*,因为向量是矩阵的特例,它也是Eigen隐式处理的,所以矩阵向量积实际上只是矩阵积的一个特例,向量外积也是如此。因此,所有这些情况仅由以下两个操作处理:

  A.双目操作符*:a*b

  B.组合操作符*=:a*=b(相当于a=a*b)

  例如:

  运行结果:

  (2) 如果您阅读了关于表达式模板的上一段,并且担心执行m=m*m可能会导致别名问题,那么现在请放心:eigen将矩阵乘法视为一种特例,并在这里引入一个临时函数,因此它将m=m*m编译为:

  tmp=m*m;

  m=tmp;

  (3) 如果您知道您的矩阵产品可以安全地计算到目标矩阵中,而不存在混叠问题,那么您可以使用noalis()函数来避免临时的情况,例如:

  c.noalias()+=a*b;

  7、点积和叉乘:

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

  (2)记住,叉乘只适用于大小为3的向量,点积适用于任何大小的向量。当使用复数时,Eigen的点积在第一个变量中是共轭线性的,在第二个变量中是线性的。

  (3)例如:

  运行结果:

  8、基本算术归约操作:

  (1) eigen还提供了一些约简操作,将给定的矩阵或向量简化为单个值,如其所有系数的和(由sum()计算)、乘积prod()或最大maxCoff()和最小minCoff()。

  (2) 如后面我们将看到的,一个矩阵的迹(由函数trace()返回)是对角线系数的和,也可以像有效地使用a. diagonal ().sum()来计算。

  (3) 还存在minCoff和maxCoff函数的变体,它们通过参数返回各自系数的坐标。

  (4)例如:

  运行结果:

  9、操作校验:

  (1) eigen检查您执行的操作的有效性。如果可能,它会在编译时检查它们,从而产生编译错误。这些错误消息可能很长,而且很难看,但是eigen用UPPERCASE_LETTERS_SO_IT_STANDS_OUT来写重要的消息。例如:

  Matrix3f m;

  Vector4f v;

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

  (2) 当然,在许多情况下,例如在检查动态大小时,不能在编译时执行检查。eigen然后使用运行时断言。这意味着,如果程序在“调试模式”下运行,它将在执行非法操作时带有错误消息而中止,如果断言被关闭,程序可能会崩溃。

  MatrixXf m(3,3);

  VectorXf v(4);

  v = m * v; // Run-time assertion failure here: "invalid matrix product"

  重要说明:本文的程序均在visual studio 2013上调试通过。

  (未完待续)

扫二维码用手机看