Fortran经过这么多年的优化，而且是纯compiled language为什么Julia速度能超过它，是包的算法的优势还是语言本身的优势？
Julia用的LLVM也是多年优化的编译器。如果给编译器的信息一样，没有理由Julia会慢。Julia会推导出变量的数据类型再转化成LLVM IR（然后交给LLVM去优化并且产生native code）。
比如说
julia> function fib(n)
a, b = 1, 1
for i in 1:n-1
a, b = b, a+b
end
b
end
fib (generic function with 1 method)
julia> @code_typed fib(5)
CodeInfo(
3 1 ── %1 = (Base.sub_int)(n, 1)::Int64 │╻ -
│ %2 = (Base.sle_int)(1, %1)::Bool ││╻╷╷╷ Type
│ (Base.sub_int)(%1, 1)::Int64 │││╻ unitrange_last
│ %4 = (Base.ifelse)(%2, %1, 0)::Int64 ││││
│ %5 = (Base.slt_int)(%4, 1)::Bool ││╻╷╷ isempty
└─── goto #3 if not %5 ││
2 ── goto #4 ││
3 ── goto #4 ││
4 ┄─ %9 = φ (#2 => true, #3 => false)::Bool │
│ %10 = φ (#3 => 1)::Int64 │
│ %11 = (Base.not_int)(%9)::Bool │
└─── goto #10 if not %11 │
5 ┄─ %13 = φ (#4 => 1, #9 => %14)::Int64 │
│ %14 = φ (#4 => 1, #9 => %16)::Int64 │
│ %15 = φ (#4 => %10, #9 => %22)::Int64 │
4 │ %16 = (Base.add_int)(%13, %14)::Int64 │╻ +
│ %17 = (%15 === %4)::Bool ││╻ ==
└─── goto #7 if not %17 ││
6 ── goto #8 ││
7 ── %20 = (Base.add_int)(%15, 1)::Int64 ││╻ +
└─── goto #8 │╻ iterate
8 ┄─ %22 = φ (#7 => %20)::Int64 │
│ %23 = φ (#6 => true, #7 => false)::Bool │
│ %24 = (Base.not_int)(%23)::Bool │
└─── goto #10 if not %24 │
9 ── goto #5 │
6 10 ─ %27 = φ (#8 => %16, #4 => 1)::Int64 │
└─── return %27 │
) => Int64
可以看出Julia把每一步的数据类型都推导出来了，可以说这样的代码和Fortran或者C几乎等价，因为都有着完全的数据类型信息。那么Julia给编译器提供的信息并不比C/Fortran那种编译语言差。
julia> using BenchmarkTools
julia> @btime fib(25)
1.882 ns (0 allocations: 0 bytes)
121393
julia> @code_native fib(25)
.text
addq $-1, %rdi
testq %rdi, %rdi
jle L51
movl $1, %ecx
movl $1, %edx
nopw %cs:(%rax,%rax)
movq %rcx, %rax
addq %rdx, %rax
movq %rdx, %rcx
movq %rax, %rdx
addq $-1, %rdi
jne L32
retq
L51:
movl $1, %eax
retq
nopl (%rax)
;}
3 个赞
算法当然是一样的啦。毕竟显式Runge-Kutta法也没有什么花样。有些细节不太一样，不过并不重要。不过，对于微分方程来说就算所有东西完全一样，Julia这种动态语言会有优势的，而Fortran这种静态编译语言会有劣势。拿常微分方程据例子 u' = f(u, t)，Runge-Kutta法（简写为RK）就调用很多次这个函数 f。如果 f 很小，Julia就会把 f 给inline（基本上就是把 f 给复制粘贴进去）到整个解法里面，那就不用call stack什么的了，所以Julia会更快。MATLAB或者Python在微分方程就会有很大的劣势（如果不JIT编译的话），因为 f 这个函数是被解释运行的。那么RK就会在 f 上浪费非常多的时间。
1 个赞
其实我觉得速度不是Julia的最大优点（虽然是个很大的优点）。Julia因为generic编程，你可以把很多东西拼在一起，而且他们之间不用知道彼此的存在。比如说，你可以把 Unitful.jl 和 DifferentialEquations.jl 一起用，你解微分方程的时候就有了单位了。你还以把它再放到 Flux.jl 里的layer，那么你就把有单位的微分方程用到机器学习了。用Fortran，C++，Python什么的都会很难做到这一点。
你做 finite difference 的话，如果用手写的loop来做operator，在Julia里你只要定义 LinearAlgebra.mul! 那么你就可以用IterativeSolvers.jl的迭代算法了。因为Julia的multiple dispatch，灵活性可是Python都比不上的。
5 个赞
我看到有些人在用Julia重写blas等线性代数库，是否这些库用julia写出来性能会更高。按理说这些库都是有C和asm写的，用Julia写为什么会更好？
C和Asm的BLAS库只能运算 Float64, Float32, Complex{Float64} 和 Complex{Float32} 。如果你用 GitHub - JuliaPhysics/Measurements.jl: Error propagation calculator and library for physical measurements. It supports real and complex numbers with uncertainty, arbitrary precision calculations, operations with arrays, and numerical integration. 这种库就会很慢。
至于绝对性能，应该不会更高。但是如果是完全Julia写的线性代数库的话，可以通过调节block矩阵的大小来进一步优化。所以，BLAS应该不会更快，但是LAPACK如果是Julia的话，会更快。
2 个赞
昨天读了一下他们12年的论文，实际上Julia在设计之初认真考虑了哪些动态性是我们不需要的，哪些是必要的。实际上大部分的程序都不需要那么动态，这也是为什么他们要重新设计一个语言，而不是去给Python做优化。有些优化是被语言本身的动态性给限制了的，技术上不可行靠砸钱可没用。
然后在我们最近的实践里，我们做了一个叫 LuxurySparse.jl 的稀疏矩阵库，针对例如 general permutation matrix等特殊矩阵进行优化，由于多重派发，你可以直接把别的东西拼过来，只写算法（我们矩阵的打印，以及和其它矩阵的交互等等都是拿的现成的，其实我自己也不知道别人都是怎么实现的）然后扩展一个更优的东西。而这个东西是一个external package，而不是要写进core里面的东西。
比起没有generic function的实现，你会发现这个的使用体验非常的统一。然后我们还可以去按照类型参数dispatch，这样例如 IMatrix （单位矩阵）我们可以在计算的时候不需要就什么都不做，而需要的时候（例如 kronecker乘）就手工构造（你可以针对不同的稀疏矩阵格式和dense matrix进行专门的构造方法）然后实现一个generic function的具体方法。
我觉得当我开始使用Julia的时候很大的一个感受就是Julia里面各种各样的数组类型和矩阵类型简直太多了，但是你使用他们毫无门槛，因为interface高度统一。
Julia的这个duck type比较彻底，没有别的给你用也。当然我觉得有时候会导致你发现某个类型虽然是 AbstractMatrix 但是有个什么方法没定义（因为可能作者本人不建议用这个矩阵类型做这种运算），但是不妨碍你用来干别的。
其实Python里面我们有时候虽然说duck type写还是按照C++来（因为只是做prototype），然后可能基类里面就有一堆 NotImplementedError ，如果这个mixin做的不是很细，那不管怎样你都至少要把上面要求实现的方法给弄一遍。
duck type目前我觉得使用体验还是非常好的。当然也不能说这个方案一定适合工程。
1 个赞
绝对性能在小矩阵上更好。（基于目前的benchmark）因为有更好的kernel？
Julia未来的优化空间还有多大？
这个问题太泛了。说实话我也不做编译器优化。不过看到discourse上Stefan说Julia可以使用很多aot和动态语言的技术。然后也有很多人（在discourse）说1.0其实还有很多优化没做（尽管换了全新的optimizer）接下来还会有很多优化出来（甚至1.0的部分benchmark有regression）
不过我目前是比较希望能尽快进一步减少jit overhead（这不是不可能），这样使用和开发体验会好一些，然后这样就能用cli了，现在repl搞这么好很大一个原因是cli不行。
虽然是动态语言，但是可以动态硬编码程序。
1 个赞
不是的，小矩阵性能好是因为Julia可以得到更多的信息（矩阵大小）从而vectorize得更彻底。如
julia> using StaticArrays
julia> m2 = @SMatrix [ 1 3 ;
2 4. ];
julia> m1 = copy(m2);
julia> @code_native m2*m1;
.text
; Function * {
; Location: matrix_multiply.jl:9
; Function _mul; {
; Location: matrix_multiply.jl:75
; Function macro expansion; {
; Location: matrix_multiply.jl:78
; Function mul_unrolled; {
; Location: matrix_multiply.jl:123
; Function macro expansion; {
; Location: matrix_multiply.jl:138
; Function *; {
; Location: matrix_multiply.jl:9
vbroadcastf128 (%rsi), %ymm0 # ymm0 = mem[0,1,0,1]
vmovddup (%rdx), %xmm1 # xmm1 = mem[0,0]
vmovddup 16(%rdx), %xmm2 # xmm2 = mem[0,0]
vinsertf128 $1, %xmm2, %ymm1, %ymm1
vmulpd %ymm1, %ymm0, %ymm0
vbroadcastf128 16(%rsi), %ymm1 # ymm1 = mem[0,1,0,1]
vmovddup 8(%rdx), %xmm2 # xmm2 = mem[0,0]
vmovddup 24(%rdx), %xmm3 # xmm3 = mem[0,0]
vinsertf128 $1, %xmm3, %ymm2, %ymm2
vmulpd %ymm2, %ymm1, %ymm1
;}}
; Function macro expansion; {
; Location: float.jl:395
vaddpd %ymm1, %ymm0, %ymm0
;}}}}
vmovupd %ymm0, (%rdi)
movq %rdi, %rax
vzeroupper
retq
nopw %cs:(%rax,%rax)
;}
但是 StaticArrays 只对非常小的矩阵好。8x8的矩阵就不能像那样优化了。
julia> m2 = @SMatrix rand(8,8); m1 = copy(m2);
julia> @code_native m2*m1;
.text
; Function * {
; Location: matrix_multiply.jl:9
pushq %r14
pushq %rbx
subq $520, %rsp # imm = 0x208
movq %rdi, %rbx
; Function _mul; {
; Location: matrix_multiply.jl:75
; Function macro expansion; {
; Location: matrix_multiply.jl:78
movabsq $julia_mul_unrolled_37105, %rax
leaq 8(%rsp), %r14
movq %r14, %rdi
callq *%rax
movabsq $__memcpy_avx_unaligned_erms, %rax
;}}
movl $512, %edx # imm = 0x200
movq %rbx, %rdi
movq %r14, %rsi
callq *%rax
movq %rbx, %rax
addq $520, %rsp # imm = 0x208
popq %rbx
popq %r14
retq
nopw %cs:(%rax,%rax)
;}
而且对于线性代数来说，我觉得小矩阵应该是像 36x36 这样的大小（如下）。因为高层线性代数，比如说LU分解，QR分解，都是把大矩阵block起来再用Level 3 BLAS的。更多资料请见 http://www.netlib.org/utk/papers/scalapack/node19.html
julia> A = rand(50,50); q = qr(A);
julia> q.T |> size
(36, 50)
这个大小上 StaticArrays 没有任何优势。
julia> m2 = @SMatrix rand(36,36); m1 = copy(m2);
julia> m3 = rand(36,36); m4=copy(m3);
julia> @btime $m3*$m4;
6.155 μs (1 allocation: 10.25 KiB)
julia> @btime $m2*$m1;
7.136 μs (0 allocations: 0 bytes)
1 个赞
我是在这里看到的optimal kernel以及他自己写的东西和open blas以及static array的比较。好像并不是简单的靠dependent type（static array这种思路）？感觉是更好的simd？
我们在写Yao.jl的时候你可以看到我们也有用static array，但是实际上发现它的效果没有在Julia里面手动permutation的好（手写的和generated出来的差不多，不过我们的比较少）。这个工作还在实验中，我觉得还是可以期待一下。
Julia的SIMD没有OpenBLAS或者BLIS手写的好。他也是从矩阵大小上下工夫的。他的库算出比较好的block大小，然后生成的kernel。相比之下BLIS或者OpenBLAS只有4种microkernel。
okay，也就是说他可以生成大量的kernel？所以在他展示的benchmark里面更好？大矩阵应该就没有什么差别了吧？
BLAS在CPU上没有什么优化的空间了。基本上都能在大矩阵上得到几乎CPU的峰值速度的。Julia可以优化的方向就是更加紧密地连接底层Level 3 BLAS和高层线性代数（LAPACK）的联系。毕竟只要调整高层线性代数的block大小，并且让BLAS对那个大小优化，那么可以让Julia的高层线性代数非常的快。
1 个赞
Okay, 不过感觉我们现在的体验还是在利用type优化一些特别的矩阵类型上，也比较好实现。然而其实我一直没弄明白为啥 StaticArray 没我们手写的permutation好。我改天再开个帖子讨论好了。
好的，直接看看生成的LLVM和汇编应该会比较明显。
StaticArrays are built on tuples which have inference limits. As the compiler improves this could improve. One big thing that is missing is the ability to recognize tuples as homogeneous so that full inference doesn’t have to happen.
4 个赞
You can read Chinese???
