极限学习机的matlab程序_Julia, Python(Numba), C++, MATLAB 性能测试之“冰雹数”

最新推荐文章于 2024-03-08 14:30:08 发布

weixin_39797324

最新推荐文章于 2024-03-08 14:30:08 发布

阅读量325

点赞数

文章标签：极限学习机的matlab程序

本专栏由"MATLAB"更名为"MATLAB or Julia", 注意这里的"or"并不是"live or die"里面做选择的意思, 而是数学里面的并集的意思.

然后我搜索了一下关键词"MATLAB", 发现我的专栏由原先的第二, 下降到了第七, 第六名是"1关注, 0文章"的专栏, 知乎的排名规则很迷啊!!!

搜索了一下关键词"Julia", 发现我的专栏竟然排到了第一, 尽管我的专栏没有写过一篇Julia主题的文章,知乎的排名规则很迷啊!!!

为了对得起这个Julia排名, 我决定写一篇关于Julia的文章. 先声明啊, 我是从11月28日开始学Julia的, 写的Julia程序肯定是比较幼稚的, 希望熟悉Julia的读者包涵与指正.

想想我学Julia的初衷是因为看到了一个Julia的官方性能测试报告, 我被它的高速与代码简洁所吸引, 然后就开始了学习Julia之旅. 既有C/C++的性能, 又有MATLAB或Python的代码简洁度, 这完全是我梦寐以求的编程语言啊!

性能测试的链接:

Julia Micro-Benchmarksjulialang.org

Julia排名第二, 而MATLAB的排名是中间靠后的, MATLAB的这个排名符合我的使用经验.

然后一对比Julia与MATLAB的测试代码, 我被震惊到了, 两者的代码相似度极高!

也就是说, 如果我会使用MATLAB, 那么我能快速学习使用Julia.

成本低(学习成本), 收益很高(因为运行速度快), 有什么理由不学Julia呢?

以上是我学习Julia的理由.

-----------------------正题分割线-----------------------------------------------------

作为一个做事严谨的人(疑心病患者晚期)对Julia官网的benchmark的还是有怀疑的, 怀疑官方特意挑选了对Julia有利的测试例子, 又或者故意将其他语言的速度拖慢.

于是, 我想要对Julia的性能进行亲手动手测试.

测试例子怎么选呢?

我首先想到就是Project Euler第14题.

Problem 14 - Project Eulerprojecteuler.net

Longest Collatz sequence
The following iterative sequence is defined for the set of positive integers: n → n/2 ( n is even) n → 3 n + 1 ( n is odd)
Using the rule above and starting with 13, we generate the following sequence:
13 → 40 → 20 → 10 → 5 → 16 → 8 → 4 → 2 → 1
It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.
Which starting number, under one million, produces the longest chain? NOTE: Once the chain starts the terms are allowed to go above one million.

这题我思考了很久, 代码也优化了很多次, 还在知乎上进行了讨论:

Matlab做project euler 14怎样提速？www.zhihu.com

@Falccm 的代码差不多是MATLAB做这题的性能极限了.(我试图提速他的这个代码, 每次都无法成功)

于是, 我将他的代码作为MATLAB的benchmark.

Falccm的MATLAB代码:

function [f,lmax] = pe14(n)
[l(n), l(1), lmax] = deal(0, 1, 0);
for k = 2:n
    a = k; lk = 0;
    while a >= k
        lk = lk + 1;
        if a == floor(a/2)*2, a = a*.5;
        else                  a = a*3 + 1; end
    end
    lk = lk + l(a);
    l(k) = lk;
    if lk > lmax, lmax = lk; f = k; end
end

测试一下速度(注意, 这里用的是1e8, 而不是原题要求的1e6, 原因是后者运算时间大概为0.06 秒, 太短了, 每次运行的相对浮动较大):

在MATLAB的R2018b上运行.

>> tic;
[f,lmax] = pe14(1e8)
toc;


f =
    63728127
lmax =
   950
时间已过 4.134375 秒。

我将Falccm的MATLAB直接改写为Julia语言:

function euler14(N::Int)
    len_star = 1
    n_star = 1
    dp = zeros(Int16, N)
    dp[1] = 1
    for n=2:N
        t = n
        cnt = 0
        while t >= n
            if t % 2 == 0
                t = div(t, 2)
            else
                t = 3*t+1
            end
            cnt += 1
        end
        dp[n] = dp[t] + cnt
        if dp[n] > len_star
            len_star = dp[n]
            n_star = n
        end
    end
    n_star, len_star
end

测试(在Julia 1.0.1上运行):

julia> @time euler14(10^8)
  1.920117 seconds (10 allocations: 190.735 MiB, 1.67% gc time)
(63728127, 950)


julia> using BenchmarkTools

julia> @benchmark euler14(10^8)
BenchmarkTools.Trial:
  memory estimate:  190.74 MiB
  allocs estimate:  5
  --------------
  minimum time:     1.880 s (0.62% GC)
  median time:      1.915 s (3.53% GC)
  mean time:        1.915 s (3.00% GC)
  maximum time:     1.951 s (4.76% GC)
  --------------
  samples:          3
  evals/sample:     1

Julia的速度是MATLAB的4.134375/1.920117 = 2.15倍.

代码复杂度差不多, 却获得了速度的提升.

然后我使用位运算, 以及其他一些技巧, 进一步提速.

function euler0014(N::Int)
    len_star = 1
    n_star = 1
    dp = zeros(Int16, N)
    dp[1] = 1
    dp[2] = 2
    for n=3:2:N
        t = n
        t += (t + 1) >> 1
        cnt = 2
        while t >= n
            if t & 1 == 0
                t >>= 1
                cnt += 1
            else
                t += (t + 1) >> 1
                cnt += 2
            end
        end
        dp[n] = dp[t] + cnt
        dp[n+1] = dp[div(n+1, 2)] + 1
        if dp[n] > len_star
            len_star = dp[n]
            n_star = n
        end
    end
    n_star, len_star
end

测试(在Julia 1.0.1上运行):

julia> @time euler0014(10^8)
  1.438625 seconds (10 allocations: 190.735 MiB, 2.77% gc time)
(63728127, 950)

julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
  memory estimate:  190.74 MiB
  allocs estimate:  5
  --------------
  minimum time:     1.375 s (0.83% GC)
  median time:      1.438 s (2.92% GC)
  mean time:        1.432 s (3.19% GC)
  maximum time:     1.477 s (4.75% GC)
  --------------
  samples:          4
  evals/sample:     1

Julia的速度是MATLAB的4.134375/1.438625 = 2.87倍.

千万别问我为什么不对MATLAB版的程序进行位运算提速, 读者有兴趣动手试一下就明白了, MATLAB的位运算速度感人.

最后, 上大杀器: C++.

使用的是上面知乎讨论帖 @rsa 的C++代码, 我添加了计时功能, 并且修改了2处bug(变量"l"的初始值由0改为了1, "f[1]"的初始值由0改为了1):

#include<cstdio>
#include<time.h>

const int n = 100000000;
int f[n + 1];

int main() {
	clock_t tic, toc;
	tic = clock();
	int a = 1, l = 1;
	f[1] = 1;
	for (int i = 2;i <= n;i++) {
		long long x = i;
		int t = 0;
		while (x >= i) {
			if (x & 1)x += x + 1 >> 1, t += 2;
			else x >>= 1, t++;
		}
		f[i] = t += f[x];
		if (t>l)a = i, l = t;
	}
	toc = clock();
	printf("%d, %dn", a, l);
	printf("Use Time:%f secondsn", 1.0 * (toc - tic) / CLOCKS_PER_SEC);
	getchar();
	return 0;
}

算法和我优化的Julia代码差不多, 都考虑使用了位运算, 奇数时, 直接进行两步(因为奇数变换后, 肯定是偶数, 可以再除以2).

运行结果(VS2017社区版, release模式下):

发现这个例子里面Julia的运算速度比C++的还要快!

为了严谨, 我将我优化的Julia算法直接修改为C++代码:

#include<cstdio>
#include<time.h>

const int n = 100000000;
int f[n + 10];

int main() {
	freopen("output.txt", "w", stdout);
	clock_t tic, toc;
	int a = 1, l = 1;
	f[1] = 1;
	f[2] = 2;
	tic = clock();
	for (int i = 3;i <= n;i+=2) {
		long long x = i;
		x += (x + 1) >> 1;
		int t = 2;
		while (x >= i) {
			if (x & 1) {
				x += (x + 1) >> 1;
				t += 2;
			}
			else {
				x >>= 1;
				t++;
			}
		}
		f[i] = t + f[x];
		f[(i + 1)] = f[(i + 1) / 2] + 1;
		if (f[i] > l) {
			a = i;
			l = f[i];
		}
	}
	toc = clock();
	printf("%d, %dn", a, l);
	printf("Use Time:%f secondsn", 1.0 * (toc - tic) / CLOCKS_PER_SEC);
	return 0;
}

运行结果为:

63728127, 950
Use Time:1.869000 seconds

运行时间总结(单位是秒):

MATLAB: 4.134375
Julia(没有位运算): 1.920117
Julia(位运算): 1.438625
C++(rsa的代码): 2.18
C++(我写的代码): 1.869

为什么Julia会比C++快呢?

我也搞不太懂, Julia我是初学, C++也是个半吊子, 不好意思. 是因为自带了并行计算? 好像说不通, 因为这个算法的循环当前值是依赖于历史值的, 属于迭代算法, 并不能并行计算.

知道答案的读者可以在评论区留言.

经过一些网友的测试, 发现他们用其他C++的编译器运行速度是要快于Julia的, 所以, 我这运行速度慢, 应该是VS2017的锅, 或者是我的配置不太对.

总结:

Julia的优势是既短又快.

本文通过一个简单的例子说明了这点.

本来是想让Julia和MATLAB比较速度的, 结果发现比C++还快, 出乎的我的意料.

---2019年1月1日更新------------------------------------------------

使用宏命令"@code_warntype", 会发现变量"len_star"的类型并不"stable"

因此, 将代码改为如下:

function euler0014(N::Int)
    len_star = Int16(1)
    n_star = 1
    dp = Vector{Int16}(undef, N)
    dp[1] = 1
    dp[2] = 2
    for n=3:2:N
        t = n
        t += (t + 1) >> 1
        cnt = Int16(2)
        while t >= n
            if t & 1 == 0
                t >>= 1
                cnt += Int16(1)
            else
                t += (t + 1) >> 1
                cnt += Int16(2)
            end
        end
        dp[n] = dp[t] + cnt
        dp[n+1] = dp[div(n+1, 2)] + Int16(1)
        if dp[n] > len_star
            len_star = dp[n]
            n_star = n
        end
    end
    n_star, len_star
end

主要变化为: 将一些字面整数转化为Int16类型的, 因为在64为操作系统里面, 字面整数默认为Int64类型.

改完以后, 在使用"@code_warntype"检查一下, 发现类型全部都是"stable"了.

性能测试:

julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
  memory estimate:  190.73 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.201 s (0.98% GC)
  median time:      1.232 s (1.11% GC)
  mean time:        1.236 s (3.06% GC)
  maximum time:     1.273 s (6.39% GC)
  --------------
  samples:          5
  evals/sample:     1

添加宏命令"@inbounds":

function euler0014(N::Int)
    len_star = Int16(1)
    n_star = 1
    dp = Vector{Int16}(undef, N)
    dp[1] = 1
    dp[2] = 2
    for n=3:2:N
        t = n
        t += (t + 1) >> 1
        cnt = Int16(2)
        while t >= n
            if t & 1 == 0
                t >>= 1
                cnt += Int16(1)
            else
                t += (t + 1) >> 1
                cnt += Int16(2)
            end
        end
        @inbounds dp[n] = dp[t] + cnt
        @inbounds dp[n+1] = dp[div(n+1, 2)] + Int16(1)
        @inbounds if dp[n] > len_star
            len_star = dp[n]
            n_star = n
        end
    end
    n_star, len_star
end

性能测试:

julia> @benchmark euler0014(10^8)
BenchmarkTools.Trial:
  memory estimate:  190.73 MiB
  allocs estimate:  2
  --------------
  minimum time:     1.112 s (1.06% GC)
  median time:      1.123 s (1.16% GC)
  mean time:        1.142 s (3.33% GC)
  maximum time:     1.185 s (7.09% GC)
  --------------
  samples:          5
  evals/sample:     1

性能总结:

运行时间总结(单位是秒):

MATLAB: 4.134375
Julia(没有位运算): 1.920117
Julia(位运算): 1.438625
Julia(类型stable) 1.236
Julia(@inbounds) 1.142
C++(rsa的代码): 2.18
C++(我写的代码): 1.869

求与C++的对比, 我用的是VS2017, 速度慢, 我估计这次能真的超越C++了.

---2019年7月4日更新------------------------------------------------

最近使用Python比较多, 顺便测试了一下Python.

代码参照Julia最快的那个版本, 稍微改写一下:

import numpy as np
def euler0014(N):
    len_star = 1
    n_star = 1
    dp = np.zeros(N + 1, dtype=int)
    dp[1] = 1
    dp[2] = 2
    for n in range(3, N + 1, 2):
        t = n
        t += (t + 1) >> 1
        cnt = 2
        while t >= n:
            if t & 1 == 0:
                t >>= 1
                cnt += 1
            else:
                t += (t + 1) >> 1
                cnt += 2
        dp[n] = dp[t] + cnt
        dp[n + 1] = dp[(n + 1) >> 1] + 1
        if dp[n] > len_star:
            len_star = dp[n]
            n_star = n
    return n_star, len_star

发现N=10^8, 运行时间过长了, 手工中断, 然后, 试了一下N=10^6:

借助numba中的jit技术, 可以提速:

from numba import jit
import numpy as np

@jit
def euler0014(N):
    len_star = 1
    n_star = 1
    dp = np.zeros(N + 1, dtype=int)
    dp[1] = 1
    dp[2] = 2
    for n in range(3, N + 1, 2):
        t = n
        t += (t + 1) >> 1
        cnt = 2
        while t >= n:
            if t & 1 == 0:
                t >>= 1
                cnt += 1
            else:
                t += (t + 1) >> 1
                cnt += 2
        dp[n] = dp[t] + cnt
        dp[n + 1] = dp[(n + 1) >> 1] + 1
        if dp[n] > len_star:
            len_star = dp[n]
            n_star = n
    return n_star, len_star

对比上面两份代码, 差别就在于后者多了两行代码:

"from numba import jit"和"@jit".

测试一下性能:

简单的增加了两行代码, 速度竟然提高这么多倍, 从941 ms 到 7.42 ms.

然后测试一下N = 10^8:

1.48秒, 速度接近于Julia版的, 我已经很满意了.

运行时间总结(单位是秒), 取每种语言速度最快的那个版本:

MATLAB: 4.134375
Julia(@inbounds) 1.142
C++(我写的代码): 1.869
Python(@jit) 1.48

以后谁敢说Python速度慢?!

借助numba的@jit 很方便的提速都了Julia, C++有一个数量级的速度.

还没有考虑numba其他的高级用法, 还有Cython.

从代码简洁程度来说, Python版最简洁.

人生苦短, 我用Python!

2019-12-2更新:

大概浏览了一下numba的文档, 感觉"Python + numpy + numba"模式最适合我, 理由是:

1 Python有完整的数据分析工具链, 生态丰富(相比MATLAB, Julia), 使用方便(相比C++).

2 长时间使用MATLAB, 已经形成了向量化思维, 而numpy非常适合做向量化运算.

pandas也可以进行向量化运算, 但是numba文档明确说了, numba适合对numpy进行提速, 而对pandas的提速不理想.

1.1. A ~5 minute guide to Numbanumba.pydata.org

当然, pandas也有适合的应用, 对性能不敏感的时候或者不需要大量循环运算的时候, 用pandas就非常适合了.

3 原生态的Python速度太慢, numpy的向量化运算速度虽然快, 但是适用范围有限, 比如本文这个问题, 必须要迭代运算(当期状态依赖于之前的状态), 这个时候, numba就弥补了numpy的缺陷(不擅长迭代运算)

总结一下:

1 Python: 语言生态丰富, 使用方便

2 numpy: 向量化运算

3 numba: 迭代运算

OK, talk is cheap, show me the codes.

下面代码是我通读numba文档以后, 对上次python代码的改进(不一定是性能极限了):

from numba import njit
import numpy as np

@njit
def euler0014(N):
    len_star = 1
    n_star = 1
    dp = np.zeros(N + 1, dtype=np.int16)
    dp[1] = 1
    dp[2] = 2
    idx_last = 1
    t0_last = 2
    for n in range(3, N + 1, 2):
        t0_last += 3
        t = t0_last
        cnt = 2
        while True:
            if t & 1:
                t += (t >> 1) + 1
                cnt += 2
            else:
                t >>= 1
                cnt += 1
                if t < n:
                    break
        dp[n] = dp[t] + cnt
        idx_last += 1
        dp[n + 1] = dp[idx_last] + 1
        if dp[n] > len_star:
            len_star = dp[n]
            n_star = n
    return n_star, len_star

相比上次的代码, 主要改动点在于:

1 dp变量的类型变为np.int16, 对于N <= 10^8来说, 足够了.

Julia里面也用到了这个技巧.

2 "@jit"改为"@njit"

关于"@njit", 参考:

1.1. A ~5 minute guide to Numbanumba.pydata.org

3 其他小的加速.

运行时间总结(单位是秒), 取每种语言速度最快的那个版本:

MATLAB: 4.134375
Julia(@inbounds) 1.142
C++(我写的代码): 1.869
Python(numba) 1.13

haitao：[高性能实战案例] MATLAB,Julia, Python(Numba)zhuanlan.zhihu.com

weixin_39797324

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
极限学习机的matlab程序_Julia, Python(Numba), C++, MATLAB 性能测试之“冰雹数”

本专栏由"MATLAB"更名为"MATLAB or Julia", 注意这里的"or"并不是"live or die"里面做选择的意思, 而是数学里面的并集的意思.然后我搜索了一下关键词"MATLAB", 发现我的专栏由原先的第二, 下降到了第七, 第六名是"1关注, 0文章"的专栏, 知乎的排名规则很迷啊!!!搜索了一下关键词"Julia", 发现我的专栏竟然排到了第一, 尽管我的专栏没有写过一...
复制链接

扫一扫