特色技术
专业从事预应力结构体系的设计、施工一体化解决方案
线性代数计算库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上调试通过。
(未完待续)
扫二维码用手机看