Eigen向量和矩阵的用法1（C++）

Posted 2023-05-03

参考技术A 在Eigen中，所有的矩阵Matrix和向量Vector都是由Matrix类构造的。向量只不过是矩阵的特殊形式，只有一列（列向量）或者一行。

Matrix模板类有6个参数，其中前三个参数是必须的。前三个参数如下：

Matrix<typename Scalar,int RowsAtCompileTime,int ColsAtCompileTime >

Scalar 是 标量类型，取值可以是 float ,int double 等。

RowsAtCompileTime 和 ColsAtCompileTime 是在程序编译时就已经知道的矩阵的行数和列数。

Eigen 提供了一些常用的 定义好的类型。比如：

typedef Matrix<float,4,4> Matrix4f .

在Eigen中，列向量是默认向量，在不特别说明的情况下，向量Vector就是指的列向量。在Eigen中定义了列向量：

typedef Matrix<float,3,1> Vector3f ;

Eigen也定义了行向量：

typedef Matrix<int ,1,2 > RowVector2i ;

如果矩阵的尺寸在编译的时候是不确定的，而在运行的时候才能确定，Eigen提供了定义动态大小的方法。比如非常好用的：

typedef Matrix<double ,Dynamic,Dynamic > MatrixXd;

MatrixXd定义了任意行数和列数的矩阵，可以在运行时确定。

类似地，对于向量有：

typedef Matrix<int ,Dynamic ,1> VectorXi ;

也可以对于一个维度确定，而指定另外一个维度是动态大小的。

Matrix<float,3,Dynamic> 矩阵的行数是 3，列数不确定。

矩阵的构造，Eigen提供了默认构造函数。

Matrix3f a;

MatirxXf b;

a 是一个3 x 3的矩阵，每个元素都是未初始化的float。

b 目前是一个 0 x 0 的矩阵。

带参数的构造函数，对于矩阵，行数在列数前面，对于向量，只有向量的大小。

MatrixXf a(10,15);

VectorXf b(30);

a 是一个 10 x 15的动态大小的矩阵，分配了内存但是没有初始化。

b 是一个动态大小的向量，大小是30,分配了内存但是没有初始化。

对于维度在 4 以下的矩阵和向量，都定义了固定大小的类型。

例如可以使用，

Vector2d ; Vector3d ; Vector4d; 等来定义向量。

Matrix2f ;Matrix3f ; Matrix4f ； 等定义矩阵。

可以使用逗号初始化方式给矩阵和向量赋值。例如：

这样就将上述值赋给了矩阵，在Eigen中矩阵默认的存储方式是行优先，就是先存储行。

Eigen支持对动态大小的矩阵和向量重新指定大小。

rows() , cols() , size() 分别返回行数，列数和 元素的个数。

resize() 则可以重新指定矩阵大小。

实例如下;

矩阵乘法性能 numpy 和 eigen c++

【中文标题】矩阵乘法性能 numpy 和 eigen c++

【英文标题】：Matrix multiplication performance numpy and eigen c++

【发布时间】：2017-12-19 01:53:00

【问题描述】：

我正在尝试使用 C++ 和 numpy 比较 eigen 的矩阵乘法性能。

这里是矩阵乘法的 c++ 代码

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

using namespace Eigen;
using namespace std;


    int main()

    time_t begin,end;
    double difference=0;
    time (&begin);
    for(int i=0;i<500;++i)
    
    MatrixXd m1 = MatrixXd::Random(500,500);
    MatrixXd m2 = MatrixXd::Random(500,500);
    MatrixXd m3 = MatrixXd::Zero(500,500);
    m3=m1*m2;
    
    time (&end); 
    difference = difftime (end,begin);
    std::cout<<"time = "<<std::setprecision(10)<<(difference/500.)<<" seconds"<<std::endl;

    return 0;

使用g++ -Wall -Wextra -I "path-to-eigen-directory" prog5.cpp -o prog5 -O3 -std=gnu++0x编译

输出：

time = 0.116 seconds

这是python代码。

import timeit
import numpy as np

start_time = timeit.default_timer()
for i in range(500):

    m1=np.random.rand(500,500)
    m2=np.random.rand(500,500)
    m3=np.zeros((500,500))
    m3=np.dot(m1,m2)

stop_time = timeit.default_timer()
print('Time =  seconds'.format((stop_time-start_time)/500))

输出：

Time = 0.01877937281645333 seconds

看起来 C++ 代码比 python 慢 6 倍。有人可以提供见解我是否在这里遗漏了什么？

我正在使用 Eigen 3.3.4、g++ 编译器 (MinGW.org GCC-6.3.0-1) 6.3.0、python 3.6.1、numpy 1.11.3。使用 spyder ide 运行的 Python。使用 Windows。

更新：

根据答案和cmets，我更新了代码。

使用g++ -Wall -Wextra -I "path-to-eigen-directory" prog5.cpp -o prog5 -O3 -std=gnu++0x -march=native 编译的C++ 代码。我无法让-fopenmp 工作 - 如果我使用此标志，似乎没有输出。

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

using namespace Eigen;
using namespace std;

int main()

    time_t begin,end;
    double difference=0;
    time (&begin);
    for(int i=0;i<10000;++i)
    
    MatrixXd m1 = MatrixXd::Random(500,500);
    MatrixXd m2 = MatrixXd::Random(500,500);
    MatrixXd m3 = MatrixXd::Zero(500,500);
    m3=m1*m2;
    
    time (&end); // note time after execution
    difference = difftime (end,begin);
    std::cout<<"Total time = "<<difference<<" seconds"<<std::endl;
    std::cout<<"Average time = "<<std::setprecision(10)<<(difference/10000.)<<" seconds"<<std::endl;

    return 0;

输出：

Total time = 328 seconds
Average time = 0.0328 seconds

Python 代码：

import timeit
import numpy as np

start_time = timeit.default_timer()
for i in range(10000):

    m1=np.random.rand(500,500)
    m2=np.random.rand(500,500)
    m3=np.zeros((500,500))
    m3=np.dot(m1,m2)

stop_time = timeit.default_timer()
print('Total time =  seconds'.format(stop_time-start_time))
print('Average time =  seconds'.format((stop_time-start_time)/10000))

使用 spyder IDE 使用runfile('filename.py') 命令运行。

输出：

Total time = 169.35587796526667 seconds
Average time = 0.016935587796526666 seconds

现在使用 eigen 的性能更好，但不等于或快于 numpy。可能是-fopenmp 可以解决问题，但不确定。但是，我没有在 numpy 中使用任何并行化，除非它隐式地这样做。

【问题讨论】：

让测试运行几分钟。您提供的时间在统计上并不合理。

更新了问题。

矩阵乘法很大程度上取决于您使用的 BLAS/Linpack 库。您必须添加 numpy 和 eigen 使用哪一个。 （例如英特尔 MKL）

【参考方案1】：

您的基准有几个问题：

您正在对系统rand() 功能进行基准测试，这是非常昂贵的！ 您缺少编译器 -march=native 以获得 AVX/FMA 提升 您缺少-fopenmp 来启用多线程。

在我的 quad i7 2.6GHz CPU 上，我得到：

initial code:                       0.024s
after replacing `Random` by `Ones`: 0.018s
adding `-march=native`:             0.006s
adding `-fopenmp`:                  0.003s

矩阵有点太小，无法获得良好的多线程优势。

【讨论】：

(1) 我在 python 代码中使用 rand() 等价物，所以它也应该很昂贵（如果真的如此的话）。 (2) 将问题更新为-march=native。它使代码更快(3) 不幸的是我无法开始工作-fopenmp

但 python 的 rand() 可能更快，请检查通过将 Random 替换为 Ones 以消除此来源，然后请指定您的 CPU。我的已经 4 岁了，所以你应该能够获得类似的性能。我们目前观察到一个大于 5 的因子（即使没有 openmp），所以您这边一定有问题。