Eigen-矩阵和向量运算

2024-02-29 06:12
文章标签 运算 矩阵 向量 eigen

本文主要是介绍Eigen-矩阵和向量运算,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

矩阵和向量运算

  • 一、概述
  • 二、矩阵加减法
  • 三、标量乘法和除法
  • 四、关于表达式模板的注意事项
  • 五、换位和共轭
  • 六、矩阵-矩阵和矩阵-向量乘法
  • 七、点乘和叉乘
  • 八、基本算术约简运算
  • 九、操作的有效性

一、概述

Eigen 通过重载常见的c++算术运算符(如+,-,*)或通过特殊方法(如 dot(), cross() 等)提供 矩阵/向量 算术运算。对于Matrix类(矩阵和向量),操作符只被重载以支持线性代数操作。我们知道在Eigen中,向量和矩阵都是用的同一个模板 Matrix 类,所以基本相互运算的操作都是一样支持的哈。

例如,matrix1 * matrix2表示矩阵-矩阵乘积,而向量+标量就是不允许的。

这里讨论的全是线性代数的操作,如果是需要对各种数组操作,而不是线性代数,请参阅专栏后一节。

二、矩阵加减法

执行这些操作的前提是,左边和右边的矩阵或向量必须有相同的行数和列数。它们还必须具有相同的Scalar类型,因为Eigen不进行自动类型提升,就是说一个整数向量是不允许和一个浮点数类型向量运算的,这只是库的要求,数学的数全是浮点数哈。这里的操作符是:

  • 二进制运算符 +,如 a+b
  • 二进制运算符 -,如 a-b
  • 一元操作符 - ,如 -a
  • 复合运算符 +=,如 a+=b
  • 复合运算符 -=,如 a-=b
#include <iostream>
#include <Eigen/Dense>int main()
{Eigen::Matrix2d a;a << 1, 2,3, 4;Eigen::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;Eigen::Vector3d v(1,2,3);Eigen::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 -12  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

三、标量乘法和除法

标量的乘法和除法也很简单。这也是满足线性代数的要求。这里的操作符是:

  • 二进制运算符 * 如 矩阵*标量
  • 二进制运算符 * 如 标量*矩阵
  • 二进制运算符 / 如 矩阵/标量
  • 复合运算符 *= 如 矩阵*=标量
  • 复合运算符 /= 如 矩阵/=标量
#include <iostream>
#include <Eigen/Dense>int main()
{Eigen::Matrix2d a;a << 1, 2,3, 4;Eigen::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

四、关于表达式模板的注意事项

这是我们在本页解释的一个高级主题,但现在只是提一下是有用的。在Eigen中,像operator+这样的算术运算符本身不执行任何计算,它们只是返回一个描述要执行的计算的“表达式对象”。实际的计算是在稍后对整个表达式求值时进行的,通常是在operator=中。虽然这可能听起来很沉重,但任何现代优化编译器都能够优化掉这种抽象,结果是完美优化的代码。例如,当你这样做时:

VectorXf a(50)b(50)c(50)d(50);
…
A = 3*b + 4*c + 5*d;

Eigen只将它编译为一个for循环,这样数组只遍历一次。简化(例如忽略SIMD优化),这个循环看起来像这样:

for(int i = 0; i < 50; ++i)a[i] = 3*b[i] + 4*c[i] + 5*d[i];

因此,我们不需要担心在Eigen中使用相对较大的算术表达式:因为Eigen自己会提供更多的优化机会。

五、换位和共轭

矩阵或向量 a 的转置 aT、共轭 a¯ 和 共轭a *(即共轭转置) 分别由成员函数 transpose()、conjugate() 和 adjoint() 得到。

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.211,0.68) (-0.605,0.823)(0.597,0.566)  (0.536,-0.33)
Here is the matrix a^T(-0.211,0.68)  (0.597,0.566)
(-0.605,0.823)  (0.536,-0.33)
Here is the conjugate of a(-0.211,-0.68) (-0.605,-0.823)(0.597,-0.566)    (0.536,0.33)
Here is the matrix a^*(-0.211,-0.68)  (0.597,-0.566)
(-0.605,-0.823)    (0.536,0.33)

对于实矩阵,conjugate()是无运算的,因此 adjoint()等价于 transpose()。

至于基本的算术运算符,transpose() 和 adjoint() 只是返回一个代理对象,而不进行实际的转置。

如果执行 b = a.transpose(),则在将结果写入b的同时求转置。然而,这里有一个复杂的问题。如果执行 a = a.transpose(),那么在转置的计算完成之前,Eigen就开始将结果写入a。因此,指令a = a.transpose() 不会像人们所期望的那样用a的转置替换a:

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;

输出

Here is the matrix a:
1 2
3 4
and the result of the aliasing effect:
1 2
2 4

这就是所谓的混叠问题。在“调试模式”下,也就是说,当断言没有被禁用时,会自动检测到这些常见的缺陷。

对于原地换位,例如在 a = a.transpose() 中,只需使用 transposeInPlace() 函数:

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

六、矩阵-矩阵和矩阵-向量乘法

矩阵与矩阵的乘法也是用运算符*完成的。因为向量是矩阵的一种特殊情况,它们在这里也被隐式处理了,所以矩阵-向量乘积实际上只是矩阵-矩阵乘积的一种特殊情况,向量-向量外积也是。因此,所有这些情况都由两个操作符处理:

  • 二进制运算符 *,如a*b
  • 复合运算符 *= 如 a*=b (在右侧相:a*=b 等价于 a = a*b)
#include <iostream>
#include <Eigen/Dense>int main()
{Eigen::Matrix2d mat;mat << 1, 2,3, 4;Eigen::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;
}

输出

Here is mat*mat:7 10
15 22
Here is mat*u:
1
1
Here is u^T*mat:
2 2
Here is u^T*v:
-2
Here is u*v^T:
-2 -02  0
Let's multiply mat by itself
Now mat is mat:7 10
15 22

注意:如果你阅读了上面关于表达式模板的段落,并且担心 m=m*m 可能会导致混叠问题,现在放心吧:Eigen将矩阵乘法视为一种特殊情况,并在这里引入了一个临时函数,因此它将 m=m*m 编译为:

tmp = m*m;
m = tmp;

如果你知道你的矩阵乘积可以安全地计算成目标矩阵而没有混叠问题,那么你可以使用 noalias() 函数来避免临时的对象产生,例如:

c.noalias() += a * b;

有关此主题的更多详细信息,请参阅有关混叠的页面。

注:对于担心性能的BLAS用户,可以使用 c.noalias() -= 2 * a.adjoint() * b;完全优化并触发单个类似宝石的函数调用。

七、点乘和叉乘

对于点积和叉积,您需要 dot() 和 cross() 方法。当然,这个点积也可以用1x1矩阵表示为 u.adjoint() * v 。

#include <iostream>
#include <Eigen/Dense>int main()
{Eigen::Vector3d v(1,2,3);Eigen::Vector3d w(0,1,2);std::cout << "Dot product: " << v.dot(w) << std::endl;double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalarstd::cout << "Dot product via a matrix product: " << dp << std::endl;std::cout << "Cross product:\n" << v.cross(w) << std::endl;
}

输出

Dot product: 8
Dot product via a matrix product: 8
Cross product:1
-21

记住外积只适用于大小为3的向量。点积适用于任意大小的向量。当使用复数时,艾根的点积在第一个变量上是共轭线性的,在第二个变量上是线性的。

八、基本算术约简运算

Eigen还提供了一些约简操作,将给定的矩阵或向量约简为单个值,例如 sum(由 sum() 计算)、product (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

由函数 trace() 返回的矩阵的轨迹是对角系数的和,也可以使用 a.diagonal().sum() 高效地计算,我们将在后面看到。

也存在minCoeff和maxCoeff函数的变体,它们通过参数返回各自系数的坐标:

  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.68  0.597  -0.33
-0.211  0.823  0.5360.566 -0.605 -0.444
Its minimum coefficient (-0.605) is at position (2,1)Here is the vector v:  1  0  3 -3
Its maximum coefficient (3) is at position 2

九、操作的有效性

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

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

当然,在许多情况下,例如在检查动态大小时,不能在编译时执行检查。然后Eigen使用运行时断言。这意味着,如果程序在“调试模式”下运行,则在执行非法操作时程序将终止并显示错误消息,如果关闭断言,则程序可能会崩溃。

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"

这篇关于Eigen-矩阵和向量运算的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/757821

相关文章

【Prometheus】PromQL向量匹配实现不同标签的向量数据进行运算

✨✨ 欢迎大家来到景天科技苑✨✨ 🎈🎈 养成好习惯,先赞后看哦~🎈🎈 🏆 作者简介:景天科技苑 🏆《头衔》:大厂架构师,华为云开发者社区专家博主,阿里云开发者社区专家博主,CSDN全栈领域优质创作者,掘金优秀博主,51CTO博客专家等。 🏆《博客》:Python全栈,前后端开发,小程序开发,人工智能,js逆向,App逆向,网络系统安全,数据分析,Django,fastapi

uva 575 Skew Binary(位运算)

求第一个以(2^(k+1)-1)为进制的数。 数据不大,可以直接搞。 代码: #include <stdio.h>#include <string.h>const int maxn = 100 + 5;int main(){char num[maxn];while (scanf("%s", num) == 1){if (num[0] == '0')break;int len =

hdu 4565 推倒公式+矩阵快速幂

题意 求下式的值: Sn=⌈ (a+b√)n⌉%m S_n = \lceil\ (a + \sqrt{b}) ^ n \rceil\% m 其中: 0<a,m<215 0< a, m < 2^{15} 0<b,n<231 0 < b, n < 2^{31} (a−1)2<b<a2 (a-1)^2< b < a^2 解析 令: An=(a+b√)n A_n = (a +

hdu 6198 dfs枚举找规律+矩阵乘法

number number number Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others) Problem Description We define a sequence  F : ⋅   F0=0,F1=1 ; ⋅   Fn=Fn

Vector3 三维向量

Vector3 三维向量 Struct Representation of 3D vectors and points. 表示3D的向量和点。 This structure is used throughout Unity to pass 3D positions and directions around. It also contains functions for doin

8. 自然语言处理中的深度学习:从词向量到BERT

引言 深度学习在自然语言处理(NLP)领域的应用极大地推动了语言理解和生成技术的发展。通过从词向量到预训练模型(如BERT)的演进,NLP技术在机器翻译、情感分析、问答系统等任务中取得了显著成果。本篇博文将探讨深度学习在NLP中的核心技术,包括词向量、序列模型(如RNN、LSTM),以及BERT等预训练模型的崛起及其实际应用。 1. 词向量的生成与应用 词向量(Word Embedding)

【Java中的位运算和逻辑运算详解及其区别】

Java中的位运算和逻辑运算详解及其区别 在 Java 编程中,位运算和逻辑运算是常见的两种操作类型。位运算用于操作整数的二进制位,而逻辑运算则是处理布尔值 (boolean) 的运算。本文将详细讲解这两种运算及其主要区别,并给出相应示例。 应用场景了解 位运算和逻辑运算的设计初衷源自计算机底层硬件和逻辑运算的需求,它们分别针对不同的处理对象和场景。以下是它们设计的初始目的简介:

位运算:带带孩子吧,孩子很强的!

快速进制 在聊到位运算之前,不妨先简单过一遍二进制的东西。熟悉二进制和十进制的快速转换确实是掌握位运算的基础,因为位运算直接在二进制位上进行操作。如果不熟悉二进制表示,很难直观理解位运算的效果。 这里主要涉及二进制和十进制之间的互相转换。 十进制转二进制 十进制转二进制可以使用常见的 除2取余法 进行。每次将十进制除以2并记录所得余数,直到商为0,然后再将记录的余数 从下往上排列即

用Python实现时间序列模型实战——Day 14: 向量自回归模型 (VAR) 与向量误差修正模型 (VECM)

一、学习内容 1. 向量自回归模型 (VAR) 的基本概念与应用 向量自回归模型 (VAR) 是多元时间序列分析中的一种模型,用于捕捉多个变量之间的相互依赖关系。与单变量自回归模型不同,VAR 模型将多个时间序列作为向量输入,同时对这些变量进行回归分析。 VAR 模型的一般形式为: 其中: ​ 是时间  的变量向量。 是常数向量。​ 是每个时间滞后的回归系数矩阵。​ 是误差项向量,假

线性代数|机器学习-P35距离矩阵和普鲁克问题

文章目录 1. 距离矩阵2. 正交普鲁克问题3. 实例说明 1. 距离矩阵 假设有三个点 x 1 , x 2 , x 3 x_1,x_2,x_3 x1​,x2​,x3​,三个点距离如下: ∣ ∣ x 1 − x 2 ∣ ∣ 2 = 1 , ∣ ∣ x 2 − x 3 ∣ ∣ 2 = 1 , ∣ ∣ x 1 − x 3 ∣ ∣ 2 = 6 \begin{equation} ||x