Julia的执行效率与主流科学计算语言一个简单的对比

Jake · 2018 年8 月 14 日 02:59

我对比了一下几种常见的科学计算C, Python, Mathematica, Matlab的语言的速度. (顺便也可以给还没入坑的人吹一波Julia)
代码是判断1-500,000的质数(不要吐槽算法,我故意用了个慢的照顾一下C), 考虑到不同语言调用系统I/O可能方式不同,所以都没加输出.
同时还想问下这个Julia的速度/执行效率正常吗?

测试都是单核的, CPU是Intel Xeon E5-1603v3@2.8GHz 以下是结果:

1.C:(GCC: 5.4.0)
Code:

#include <stdio.h>
#define MAX 500000

int prime(int n);

int prime(int n) {
  int i;

  for (i= 2; i < n; i++)
	if (n % i == 0)
	  return 0;

  return 1;
}

int main() {
  int i;

  for (i= 2; i <= MAX; i++) prime(i);

  return 0;
}

用时:

2.Julia:(0.6.4)
Code:

function testp(x)
	for i=2:x-1
		if x%i==0
			return 0
		end
	end
	return 1
end

for i=1:500000
	testp(i)
end

用时:

3.Matlab(R2018a)
Code:

function re=testp(x)
	for i=2:x-1
		if mod(x,i)==0
			re=0;
			return
		end
	end
	re=1;
	return;
end

tic;
for i=1:500000
	testp(i);
end
toc;

用时:

如果用Matlab内置函数isprime(x):
matlab1

4.Mathematica(11.2)
Code:

gcp = Compile[{{x, _Integer}},Module[{i = 2}, While[i <= x - 1 && Mod[x, i] != 0, i++]], CompilationTarget -> "C"]
  
Table[gcp[i], {i, 500000}]; // AbsoluteTiming

用时:

5.Python (3.5.2)
Code:

def testp(x):
	for i in range(2,x):
		if(x%i==0):
			return 0;
	return 1;

for i in range(1,500000):
	testp(i)

用时:

结论: C(35.6s)> Julia(112.9s) > Matlab(232s)> Mathematica(355.2s)>Python(16min57s).

C在执行效率方面还是无与伦比啊,不得不说"你爸爸还是你爸爸". Julia虽然居于第二,但是并没我之前想象中的可以和C的速度比拟. 不知道是不是我的姿势水平太低,能不能继续提升Julia的效率. Mathematica 我在调用函数前预先compile成了C的代码,所以讲道理不是Wolfram language 的速度, 如果不compile,根本算不出来… (不过我代码就一行啊 )Python有点慢的令人发指 .

FYI

内置函数比较:
1.Mathematica的内置函数:

Table[PrimeQ[i], {i, 500000}]; // AbsoluteTiming

MMAbuiltin

2.Matlab内置函数isprime(x):
matlab1

3.Julia 内置函数Primes.isprime(x):
code

Julia prevails!

Roger · 2018 年8 月 14 日 04:29

补充一下各个语言和软件的版本？还有C语言和Julia以及其它的编译指令？

当然如果是以好用为前提的话，那么你按照自己的直觉写出来的第一手代码就是公平的了，这一点上Julia现在确实还不一定能比C，但是不会差这么多（但是这个的前提是你没有任何抽象，等有了抽象以后就还真不一定了）。

Gnimuc · 2018 年8 月 14 日 04:30

虽说Julia是可选类型，但如果要和C比效率就需要留意类型，Julia里的Int类型默认是Int64, 而C的Int是32bit Int，所以C也应该换用长整型(long long int)。或者Julia里用Cint:

julia> function testp(x)
           for i=2:x-1
               x % Cint(i) == 0 && return 0
           end
           return 1
       end
testp (generic function with 1 method)

julia> function foo()
         for i=1:500000
           testp(Cint(i))
         end
       end
foo (generic function with 1 method)

julia> @time foo()
 28.644356 seconds (2.51 k allocations: 132.750 KiB)

Roger · 2018 年8 月 14 日 04:34

不一定Cint吧，其实直接Int32 也可以的

Gnimuc · 2018 年8 月 14 日 04:35

从这个例子也能看出，Julia优化代码的曲线是很平滑的：不需要额外地调这库那库，只需要留意代码本身的问题（这个例子的根本原因是32bit整型比64bit整型的运算快），然后采用相应的解决方法(显式的写类型)。

Gnimuc · 2018 年8 月 14 日 04:36

我个人习惯把和C相关的代码都用Ctype这些alias来做。显式的说明这里是和C中的类型匹配，这样不容易出错。

Roger · 2018 年8 月 14 日 04:38

这个应该会更快一点点你试试（不过其实不特别明显好像，我@benchmark 测试只有maximum time明显下降）

function testp2(x::Int32)::Int32
       for i=Base.OneTo(x-2)
           x % (i + 1) == zero(Int32) && return 0
       end
       1
end

然后这么跑测试

function foo2()
       for i in Base.OneTo(Int32(50000))
           testp2(i)
       end
end

@benchmark foo2()

这个应该要比原来那个更快一些，大家猜猜为什么？

Roger · 2018 年8 月 14 日 04:39

最后一个我有点惊讶啊，我以为Mathematica的优化应该是不错的。你是什么版本啊？

Roger · 2018 年8 月 14 日 04:49

再polish一下，其实上面的代码为了性能加了很多约束，但是其实这样写也许更好

function testp2(x::T)::T where {T <: Integer}
       for i=Base.OneTo(x-2)
           x % (i + 1) == zero(T) && return T(0)
       end
       T(1)
end

Gnimuc · 2018 年8 月 14 日 04:50

解决方法有多种，OP的代码的主要问题是1:50000这个UnitRange默认是Int64, 而经过duck-typing之后，这个类型会传播到%，进而影响计算效率，所以既可以加上type annotation：testp2(x::Int32)::Int32来提早约束类型，也可以保留duck-typing，显式地给UnitRange加类型：

julia> function testp(x)
           for i = UnitRange{Cint}(2,x-1)
               x % i == 0 && return 0
           end
           return 1
       end
testp (generic function with 1 method)

julia> function foo()
         for i=UnitRange{Cint}(1, 500000)
           testp(i)
         end
       end
foo (generic function with 1 method)

julia> @time foo()
 28.564901 seconds (3.22 k allocations: 172.355 KiB)

具体采用那种要根据实际需求来定，我觉得如果解决了Int64和Int32这个本质问题，效率都不会差的太多。

Roger · 2018 年8 月 14 日 04:50

我其实是想说 OneTo 比 UnitRange 更快。因为这里给编译器的约束更加充分（直接写在类型里了）

Gnimuc · 2018 年8 月 14 日 04:50

这样更地道一些，不过OP初学可能得熟悉一下。这种测试还是需要各个语言的code都能做到obvious和concise，如果搞的太复杂，其它语言的用户又要站出来晒一些奇淫巧技了额。。。

Roger · 2018 年8 月 14 日 05:00

嗯，不过我觉得问题在这里：OneTo 也是纯Julia实现的（不同于向量化操作等技巧其实是在跑C）。在类型里加入约束我觉得对了解一些编译原理的人来说是比较直觉的。但是实际上如果要说obvious的代码实现的性能那么就是他原来的这个注意一下Int32就好了。

不过就是给个例子，你要想在Julia里面提高性能要怎么提高，大致是这样：不断地（用Julia自己）给编译器增加约束。还是比较容易的。

shilu1984 · 2018 年8 月 14 日 05:09

约束类型感觉跟cython一样，看来还是学习学习julia了。
另外，在C里，int和long差距并不大，我机子上一个35，一个37秒。

Roger · 2018 年8 月 14 日 05:11

是的，这是几乎所有的语言提高性能的方式。因为编译器也是一些数学推导一样的东西。但是和cython的区别在这些都是纯Julia的类型。cython更加像C一些（个人感觉上）。

实际上正是因为Julia的开发者注意到很多语言无法保留尽可能多的信息给编译器，才有的Julia。

Roger · 2018 年8 月 14 日 05:23

Int64是long long吧？

https://en.cppreference.com/w/c/language/arithmetic_types

按照标准，long 这些类型只规定了最低多少位，应该用 stdint.h 里的 int64_t 和 int32_t

#include <stdio.h>
#include <stdint.h>
#define MAX 50000

int prime(int32_t n);

int prime(int32_t n) {
  int32_t i;

  for (i= 2; i < n; i++)
	if (n % i == 0)
	  return 0;

  return 1;
}

int main() {
  int32_t i;

  for (i= 2; i <= MAX; i++) prime(i);

  return 0;
}

这个时候32和64就是有差别的。

shilu1984 · 2018 年8 月 14 日 07:07

我也来凑个热闹。

机子上只有matlab 2018a, python3.6 + numba0.39 + cython0.28.5，vs 2010，其它没有装，主要对比的是numba和cython还有Julia
cpu为i7-4712h

代码稍微改动一下，要不然我的c和numba偷懒不干活，直接跳过去了，不知道是不是main函数循环对结果没影响被优化了。

c：

int main() {
    int i;
    int n = 0;

    for (i= 2; i <= MAX; i++) {
        n += prime(i);
    }
    return n;
}

cython:

cdef int testp(int x):
    cdef int i
    for i in range(2, x):
        if x % i == 0:
            return 0
    return 1

def runme():
    cdef int i
    cdef int n = 0
    for i in range(1, 500000):
        n += testp(i)
    return n

numba：

import numba as nb


@nb.jit(nb.int32(nb.int32), nopython=True)
def testp(x):
    for i in range(2, x):
        if x % i == 0:
            return 0
    return 1


@nb.njit
def runme():
    n = 0
    for i in range(1, 500000):
        n += testp(i)
    return n

Julia：

function testp(x)
    for i=2:x-1
        x % Cint(i) == 0 && return 0
    end
    return 1
end
 
function runme()
    for i = 1:500000
        testp(Cint(i))
    end
 end

C的运行时间是为30.13s
matlab的运行时间是为212.34s
python+cython的运行时间是为37.3s
python+numba的运行时间为43.8s
Julia的运行时间为33.67s
报告完毕！

BOBO · 2018 年8 月 14 日 07:07

如果要match数据类型的话，这种比较就显得很不自然了。另外，compilation的选项不清楚以及MATLAB貌似最新版会自动用多核。这种语言之间的速度比较特别难，往往很有争议。

Julia的好处在于语法很friendly, 而且速度可以是scalable。因为大部分速度对比的例子都是code不超过几十行，所以对实际应用implications不大。我这里的scalable是指code量和运算量同时增加，这个时候Julia就很有优势。

ps: Mathematica不是主流科学计算语言。。。。。

BOBO · 2018 年8 月 14 日 07:23

纯for loop的写法，compiled language一般占很大优势，Cython和Matlab的mex与C的速度比较差不多就是C和C自己比较。

Julia和Matlab的优势在数学语法的舒适度和性能并存，但说性能好多硬件语言速度要快的多。

honghuzi · 2018 年8 月 14 日 09:06

你的 c 版本运行快慢也和编译器有关吧。我用 clang，几乎是秒算完。不过你的 julia 版本我电脑上(v0.6.4)要85秒。