关于几个数学工具，矩阵运算速度的初步测试

gps99 · 发表于 2022-4-4 18:49:57

本帖最后由 gps99 于 2022-4-4 19:00 编辑

喔噻，，我居然有耐心等这个4k的结果出来溜溜的跑了十几分钟。、。

是的，十几分钟一致在想，这也许是今后自己手写数值计算的下场。

著名前辈尚且如此，我等千万不要自己造轮子！！！千万千万！！

zoziha · 发表于 2022-4-4 20:18:57

本帖最后由 zoziha 于 2022-4-4 20:28 编辑

你是使用纯Fortran的Lapack吗？从我的认知上看，Numpy后端矩阵求解使用的好像是OpenBLAS，在Fortran中也能方便地使用OpenBLAS，理论上Fortran写法(OpenBLAS)和Numpy写法性能表现应该一致。
Python代码越少，OpenBLAS调用相关代码越多，Python的Numpy写法效率应该更高。

同样的Intel的Fortran编译器使用自家的MKL性能应该也是飞起。（据我所知，matlab使用的是intel家的MKL库，这谁也打不过啊，intel家软硬件结合，cpu届应该无敌）
自实现算法与使用著名函数库（如OpenBLAS、MKL、IMSL、Petsc、Numpy、fftw等），后者往往是许多人集体工作的结晶，也普遍存在大量优化，是典型的商业或社区驱动“产品”，而自实现算法往往只能在代码层面优化，难以完美掌控底层编译器行为，性能一般是不如著名函数库。但是，这也是跟着具体的业务和需求走的，不需要长时间大量计算的代码，也不一定极致追求调用性能最好的库，盲目引入外部库也可能是一种不好的实践。

互联网时代，Python自带优势，借着这东风，已经很强了，我能看到很多c/cpp/fortran代码被Python社区（从scipy/numpy源码中可见）很好地利用起来了，但编译型语言始终代表性能top1行列。希望Fortran能迎头赶上。

gps99 · 发表于 2022-4-4 23:12:10

zoziha 发表于 2022-4-4 20:18
你是使用纯Fortran的Lapack吗？从我的认知上看，Numpy后端矩阵求解使用的好像是OpenBLAS，在Fortran中也能 ...

是的，纯fortran的Lapack。在 netlib-Lapack 3.9.1 下载源码，按包内的makefile本地编译得到三个.a静态库.

从几个软件所各自代表的LinearAlgabra包角度看，实测的差距感觉比预想的大。纯矩阵效率上自实现代码 <- 专用包Lapack... <- 商业专业包MKL Matlab包... 三个台阶泾渭分明，居然能有数量级的差距的确很意外。当然实际场景中综合各种运算，不会有这类白手套测试的差距明显，正如您所说盲目引入外部包不一定是明智的选择，结构需综合考虑。

另外今天还大致跑了其他几个测试，包括纯for迭代的整形、浮点速度，fibnaci的递归效率，。。有一个明显的感觉就是：编译语言的各项性能比较稳定（gcc，gfortran .VS. matlab， python+numba）。在递归、for迭代、整形、浮点的效率上，c/fortran比动态语言有几十倍的压倒性优势。综合性能上尤其是C，绝对实至名归，一个字稳。

唐汉 · 发表于 2022-4-17 07:49:33

Julia那个应该是把时间花在类型转换上了。类似的“benchmark”网上有很多，但是大多数都是疑点重重的，关键在于怎么用不同的语言写出等效的代码。比如说你的julia代码，M乘以Mi的时候是一个整数矩阵乘以一个浮点数矩阵，每一个元素相乘的时候都要做一次类型转换。但是如果你提前把M转化成浮点数的矩阵，速度就上来了。python或者是numpy应该是把这一步给包起来了。我测试了一下

（1）Python：

[Python] 纯文本查看 复制代码

import numpy as np

matrix_size = 4000
int_range = (1, 99)
M = np.random.randint(*int_range, (matrix_size, matrix_size))
Mi = np.linalg.inv(M)
Re = np.dot(M, Mi)

（2）Julia：（论坛里貌似没有Julia的高亮）

[Python] 纯文本查看 复制代码

matrix_size = (4000, 4000)
int_range = 1:99
M = rand(int_range, matrix_size...)
Mᵣ = Array{Float32}(M)
Mᵢ = inv(Mᵣ)
Rₑ = Mᵣ * Mᵢ

结果是差不多的julia稍微快一点：
$ time python test.py

real 0m3.473s
user 0m22.469s
sys 0m2.250s

$time julia test.jl

real 0m2.202s
user 0m8.245s
sys 0m2.246s

如果换成Float64 julia速度会降下来一点但是也就是4秒多一点，绝对不至于跌出一个数量级来。然后调用OpenBLAS的时候矩阵超过一定大小应该就会自动超线程了，Julia和python在这个层面上是否做了一样的事情我也不太清楚。

大部分情况下如果代码耗时在同一个硬件上出现了数量级的差距，那基本上都是算法上不等效（当然你如果用纯python不调用任何库那肯定会慢出天际……）。

Fortran那个代码我暂时还没测试，有空了我再来玩玩~

gps99 · 发表于 2022-4-20 19:25:16

本帖最后由 gps99 于 2022-4-20 19:28 编辑

唐汉发表于 2022-4-17 07:49
Julia那个应该是把时间花在类型转换上了。类似的“benchmark”网上有很多，但是大多数都是疑点重重的，关键 ...

感谢您对Julia的指导。
前面的Jl代码没有指定类型，确实影响了julia的测试速度。指定Array{Float64}后，Julia速度提高近20倍。核Numpy基本一致，甚至略快。。果然Julia是需要注入技巧，才有速度。
同是 float64格式下：
Julia 3.3 s
Numpy 3.4 s

Julia 代码如下：

[Fortran] 纯文本查看 复制代码

@time begin
global M = rand(1:99, 4000, 4000)
global M = Array{Float64}(M)
Mi = inv(M)
Re = M * Mi
println(Re[11, 11])
end

gps99 · 发表于 2022-4-20 20:16:28

本帖最后由 gps99 于 2022-4-20 20:20 编辑

Numpy 和 Julia在纯向量的优化程度上，如此的接近，应该不是巧合。也许是目前商业库能达到的某种特征水平。

过几天再把openBLAS试试，看看是一个什么结果。。。估计达不到4 s以内。

论坛内大咖云集，哪位有时间的话 c++或fortran调一下openBLAS，并优化一下。。。。。。。最后放一个EXE的附件。

测试目的只有一个： 即便是一个看似很简单的调用，我们能否实现商业库的效率？？尽力追，，能追到一个什么程度？

[线性代数] 关于几个数学工具，矩阵运算速度的初步测试