[Eigen中文文档] 矩阵与向量运算

简介: 本文章旨在提供有关如何使用 Eigen 在矩阵、向量和标量之间执行算术操作的概述和一些详细信息。

专栏总目录

英文原文(Matrix and vector arithmetic)

本文章旨在提供有关如何使用 Eigen 在矩阵、向量和标量之间执行算术操作的概述和一些详细信息。

介绍

Eigen 通过重载常见的 C++ 算术运算符(如 +-*)或通过特殊方法(如 dot()cross() 等)提供矩阵/向量算术运算。对于 Matrix 类(矩阵和向量),重载运算符仅支持线性代数运算。例如,matrix1 * matrix2 代表矩阵乘法,vector + 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 -1
 2  0
Doing a += b;
Now a =
3 5
4 8
-v + w - v =
-1
-4
-6

标量的标量乘法与除法

标量的乘法和除法也非常简单。目前支持的运算符示例如下:

运算符 示例
二元运算符 * matrix * scalar
(矩阵 * 标量)
二元运算符 * scalar * matrix
(标量 * 矩阵)
二元运算符 / matrix / scalar
(矩阵 / 标量)
复合运算符 *= matrix *= scalar
(矩阵 *= 标量)
复合运算符 /= matrix /= scalar
(矩阵 /= 标量)

代码示例:

#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中,诸如+之类的算术运算符,他们自己不执行任何操作,只是返回一个表达式对象,该对象描述了将要执行的计算操作。实际的计算发生在后面整个表达式被求值的时侯,比如使用=运算符时。虽然这听起来很繁琐,但任何现代优化编译器都能优化掉这种抽象,从而得到完美优化代码。例如:

VectorXf a(50), b(50), c(50), d(50);
...
a = 3*b + 4*c + 5*d;

Eigen会把上述表达式编译成一个循环,这个数组只遍历一次。数组循环如下所示:

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

因此,你不要害怕使用相对较大的运算表达式,这只会给Eigen更多机会进行优化。

转置与共轭

矩阵或向量 $a$ 的转置($a^T$)、共轭($\overline{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()不会得到预期的结果。

示例如下:

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

上述的问题就是所谓的混叠问题,在debug模式下,当assertion打开,这个问题可以自动检测到。(g++编译默认是debug模式,关闭需要使用-DNDEBUG选项)。

对于就地转置,可以使用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=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 -0
 2  0
Let's multiply mat by itself
Now mat is mat:
 7 10
15 22

注意:如果你阅读过上面的关于表达式模板的段落并且担心 m = m * m 会引发混淆问题,这里请放心,Eigen把矩阵乘法作为一个特殊的例子,并在此引入了一个临时变量,所以它会编译为:

tmp = m*m;
m = tmp;

如果你知道你的矩阵乘法可以安全的计算并且没有混淆问题,那么你可以使用noalias()函数来避免编译临时变量,例如:

c.noalias() += a * b;

更多细节请参考 aliasing

注意:对于担心性能的 BLAS 用户,表达式如:c.noalias() -= 2 * a.adjoint() * b;可以完全的优化并触发一个类似矩阵乘法的函数调用。

点积和叉积

对于点积和叉积,需要使用 dot()cross() 方法。当然,点积也可以像 u.adjoint()*v 一样得到一个1x1的矩阵。

示例如下:

#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 scalar
  std::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
-2
 1

注意,叉积仅适用于大小为 3 的向量。点积适用于任何大小的向量。使用复数时,Eigen的点积在第一个变量中是共轭线性的,在第二个变量中是线性的。

基本算术的简化运算

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

矩阵的迹(对角线系数的总和)可以通过函数trace()计算,也可以使用更高效的方法a.diagonal().sum()

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

  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.536
 0.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

这里的输出,自己测试和我官网给出的有出入:

  1. 浮点型数字,官网截取了前面几位,而且做了四舍五入(前文已提到过)
  2. 整型数字,和官网完全不一样(暂未发现原因,如你知道欢迎留言)

关于这里的Random()函数,源码中有解释:

Numbers are uniformly spread through their whole definition range for integer types, and in the [-1:1] range for floating point scalar types.

对于整型,在整个定义范围内均匀分布;对于浮点型,分布在[-1,1]的范围内。所以,从代码来看,官网的结果 Here is the vector v: 1 0 3 -3 是不符合的。

!!!注意,以上两个问题,整篇都有,下文不再赘述!!!

操作的有效性

Eigen会检查操作的有效性,如果有错误,它会在编译的时候产生错误提示。这些错误提示可能又长又难看,但Eigen会把重要的信息写成大写,以使其更加显眼,例如:

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

当然,在很多情况下,如检查动态矩阵的大小时,无法在编译时进行检查,Eigen会使用运行时的断言。这意味如果程序在debug模式下运行,遇到非法操作时会终止运行并打印出错误信息。如果关闭断言,程序可能会崩溃。

MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"
相关文章
|
缓存 测试技术 编译器
[Eigen中文文档] 稠密矩阵分解函数对比
本文介绍了 Eigen 为各种方阵和过约束问题提供的稠密矩阵分解的速度比较。
133 0
|
XML 并行计算 算法
[Eigen中文文档] 求解稀疏线性系统
在Eigen中,有多种方法可用于求解稀疏系数矩阵的线性系统。由于此类矩阵的特殊表示,必须特别小心以获得良好的性能。本文列出了Eigen中可用的稀疏求解器。还介绍了所有这些线性求解器共同的主要步骤。根据矩阵的属性、所需的准确度,最终用户可以调整这些步骤以提高其代码的性能。请注意,并不需要深入了解这些步骤背后的内容:最后一节介绍了一个基础例程,可轻松使用以获取所有可用求解器的性能洞察。
340 0
|
存储 算法 NoSQL
[Eigen中文文档] 稀疏矩阵操作
在许多应用中(例如,有限元方法),通常要处理非常大的矩阵,其中只有少数系数不为零。在这种情况下,可以通过使用仅存储非零系数的特殊表示来减少内存消耗并提高性能。这样的矩阵称为稀疏矩阵。
463 0
|
6月前
|
算法 Java
行列式表示只需3步,快速掌握矩阵运算!
行列式表示只需3步,快速掌握矩阵运算!
72 0
|
机器学习/深度学习 前端开发 rax
第3章 数组与矩阵——3.4 矩阵运算(1)
第3章 数组与矩阵——3.4 矩阵运算(1)
|
机器学习/深度学习 资源调度 算法
第3章 数组与矩阵——3.4 矩阵运算(2)
第3章 数组与矩阵——3.4 矩阵运算(2)
[Eigen中文文档] 无矩阵求解器
本文介绍Eigen的无矩阵求解器。
145 0
|
测试技术
[Eigen中文文档] 线性代数与分解
本节将说明如何求解线性系统,计算各种分解,如 LU、QR、SVD、特征分解……
228 0
|
机器学习/深度学习 人工智能
线性代数基础--矩阵
线性代数基础--矩阵
1228 0
【线性代数-基础理解】对换行列式两行(列),行列式变号
【线性代数-基础理解】对换行列式两行(列),行列式变号
747 1
【线性代数-基础理解】对换行列式两行(列),行列式变号