技术▍Julia 解决了 C++/Python/Matlab 的哪些痛点？

其实说Julia解决了什么痛点，回到2012年Julia第一篇论文里也许可以体会到：

论文地址：https://arxiv.org/pdf/1209.5145v1.pdf

过去动态语言大部分都是面向解释器设计的，这使得动态语言非常简单，动态。但是实际上我们真正使用的程序不需要这么动态的特性，Julia是一种新的平衡，它是面向JIT等优化技术设计的动态语言。它去掉了一些不那么重要的动态性，加入了一些限制。这是过去Python，MATLAB的局限。

从来没有觉得0.4入Julia太早，0.5入太早，0.6太早，我是从0.3最后几个月开始用的，对于科学计算领域我一直觉得它非常合适。但是现在1.0已经上线了，我觉得你们在开始关注的同时应该等一个月入，23333

Julia可以调用Python的所有东西（JuliaPy/PyCall.jl），可以调用大部分的R，所以即便有历史负担，也不必太担心，除非你的任务很急，就是一个月要搞出来什么东西。因为虽然Julia学习曲线平滑，但是想用Julia写出性能好，抽象干净的代码是需要一定时间的。Python这么简单是优点，也同时带来了缺点。

此外，Julia社区也从未说放弃Python，因为no silver bullet（没有银弹），只有在科学计算这个领域，Julia目前才是更加合适的解决方案，因为它为科学计算而生，但是在其它领域Julia就几乎没有优势了。所以你还有pyjulia来帮助你在Python里使用Julia。当然我们也许可以期待，在将来我们很多时候不需要调用Python直接使用Julia。我们最近也有一个使用Julia开发的东西，你通过这个工程看看我们相比有一定相似功能的Python库Project Q有什么优势：

Yao.jl 一个可扩展，高效的量子算法设计框架 ：https://zhuanlan.zhihu.com/p/38395675

我们（Julia中文社区，算是Julia社区的一个延伸，Julia有很多本地化组织）暂定在29日在北京中国科学院软件研究所办一次Julia User Meet Up（报名网站:Julia用户见面会（Meetups））也希望感兴趣的朋友可以过来参加。

先讲讲我所接触到的做科学计算的一些场景：

（以下任务对性能时有基本要求的，所以性能上的需求都不再提）

哈密顿量的基态计算：常见的方法有DFT，张量网络，量子蒙卡等都是可以大规模并行的方法。在工程上往往用Fortran，C++，MATLAB来写，往往需要有集群/服务器配合做大规模并行计算。此外由于量子力学本身依赖于线性代数和希尔伯特空间（这也意味着简单的矩阵运算远远不够），这一类任务严重依赖于多维数组。

量子计算的一些数值模拟/硬件交互：数学对象比较复杂，并且几乎每一种不同结构的线路经典模拟方式都不一样，有很多特化的模拟算法。此外和硬件交互需要对C有良好的兼容性。

深度学习：除了对性能上的要求，还需要能够灵活地定义算符，并且在经过PyTorch的实践之后，尽量希望能做成动态图这样更加方便的接口，但是依然需要保证性能（PyTorch选择了自己加JIT）。此外还需要有良好的GPU支持，在目前的方案里都是使用CUDA然后通过Python C API接入。

自动微分：自动微分的应用场景往往是一些比较复杂的表达式需要求导，但是又对性能有一定要求（不然也许就用更简单的有限元/符号求导了）。自动微分有两种方案：代码生成和重载。一般来讲代码生成更快一些。

然后再列几个具体的任务：

1、张量缩并：缩并指对某些角标求和，几乎所有的张量运算都可以用爱因斯坦求和记号，或者所谓的张量网络来标记，而一串张量的缩并顺序将会严重影响你的计算性能。比如下面这个按照某种顺序求和一系列高维数组，具体的求和顺序（先求和谁）是需要提前计算出来的，但是我们不希望之后每次计算都再来一遍。

2. 并行Monte Carlo抽样：Monte Carlo里例如一些基于MCMC的方法需要将大量有一定工作量（复杂度）的相同的任务（比如MCMC的一个chain）在大量进程上并行，一般来讲，我们都使用MPI。

3. 序列化：计算完毕之后我们需要保存结果，这个结果往往保存在CSV，HDFS等文件格式中。当然Python中还有numpy的二进制格式和pickle。

4. 数学对象的处理：往往对某些数学对象的拓展是通过增加少数几种特性（比如新的二元关系，新的对称性等）而获得的新的数学对象，那么这个过程中用class可能就需要mixin来做，而大量的继承就不可避免出现菱形继承（或者也称为钻石继承），这种继承关系本身在数学上可能就不存在（或者是根本不需要用）而是为了使用class人为加上去的，造成了不必要的复杂性。

这里讲的比较抽象，用Python举个例子：

我有一个矩阵类型 Matrix（data我就用numpy的来存了

import numpy as np
class Matrix:
    def __init__(self, m, n):        self.data = np.zeros(m, n)

然后接下来，我有一个对称矩阵Symmetric

class Symmetric(Matrix):
    def is_symetric(self):          return true

然后还有下三角矩阵LowerTriangle和上三角UpperTriangle

class LowerTriangle(Matrix):
     def is_tril(self):          return true
class UpperTriangle(Matrix):
     def is_triu(self):          return true

然后我有厄米矩阵

class Hermitian(Symmetric):
     def is_hermitian(self):          return true

然后还有对角矩阵Diagnal，它是Matrix，是Symmetric，也是下三角和上三角

class Diagnal(Matrix, Symmetric, LowerTriangle, UpperTriangle):
     def is_diagnal(self)           return true

然后我有一个厄米的对角矩阵(可能就是个你要用的哈密顿量呢也许)。。。

class MyMatrix(Hermitian, Diagnal):
      def is_mymatrix(self):            return true

然后请注意的是以上每一种类型都有一堆方法可以针对性地去优化特定类型的性能。所以每一个都至少有一个类似于apply之类的方法来让其可以使用特别的算法去计算。OK，现在请在大脑里解析一下它们之间的各种方法的调用关系....等到了其它场景（比如我们的量子线路模拟器里，这个会变得更复杂）。

诚然，在C++里我们用虚基类和虚方法来定义接口来搞这个事情。也不失是一种可行的方案。但是Julia提供了一种不太一样的面向对象的方式。所以反对有人说Julia没有引入新的编程范式，这是错误的。

Julia是一种从语言设计上尝试为以上场景提供更好的解决方案的一种新尝试，并且在实际应用中确实展现了预期的性能和抽象能力。以至于甚至有人提出用Julia重写“BLAS”。（目的是提供对不同性质矩阵的优化）。

所以这是为什么说：Julia's performance is not magic, its language design is.

Julia的面向对象和多态是基于多重派发的，我们姑且称之为基于多重派发的OO，它有点像是一个结合了模式匹配（pattern matching）的类型系统，但是比模式匹配要弱。Julia是不具有class的，所以我反对把Julia的代码当Python来设计。

关于这一点可以参考我的文章：浅谈Julia语言：Julia的面向对象 https://zhuanlan.zhihu.com/p/38434080

这里讲一点故事吧，其实在起初整个community也不清楚要怎么使用多重派发的，因为从未有人这么试过。所以你能在比较早期的Julia package里（那个时候还是0.3左右吧）经常看到在代码文件开头的一大串注释：某个抽象类的子类型一定要包含xx成员，定义xx方法。这就是从传统class过来的思路，在这个时候这些抽象类和文档里定义的要重载的xx方法，和xx成员就是C++里的虚基类。但是后来就不这么写了，因为Julia被设计成duck type，只要定义好了一套interface和抽象类的层次那么剩下的东西根据具体情况随意写就好。行为（behavior）定义了这个类型是什么。所以现在去读一些新的代码，你也许就会看到某个文件里（常常是core.jl）列了一堆抽象类和它们的接口。甚至连官方文档里都开始出现了interface的介绍（可惜关于怎么在语法里体现interface还没有定论）。

那么这里主要解决了之前的OOP方法里出现菱形继承的问题。然后顺便我们获得了更加轻量级的类型（struct）。配合多重派发我们可以针对不同类型派发不同的方法，于是Julia不仅仅有很快的BLAS和多维数组，还有不同类型的矩阵（比如三角矩阵，对角矩阵等等），他们可以被派发不同的方法，实际上从0.7开始通过对这些类型的重载还有很多方法直接被做成了lazy的。而在C++里这需要更复杂的expression template来实现（这个可以参考李沐在mshadow里的note）。

实际上我们在最近开发的Yao.jl中实现了一个自己的稀疏矩阵库（LuxurySparse，这个名字是为了嘲讽现在的矩阵库没我们快，当然快的代价是类型变多了，逃），针对不同类型的稀疏矩阵（例如identity，general permutation matrix等）进行了优化。而最终的性能在实际计算中的表现也远好于scipy等的实现。而这些都是在纯Julia里实现的，scipy可不是纯Python实现的。

以上我想说，这样的编程范式在数学概念上更自然（你不需要去管继承这件事情，因为也没有，可以直接按照数学上的写法来写东西），方法就是不同对象之间的interaction，而更加简单的类型也带来的更好的性能。

宏是Julia的另外一个特性，虽然继承自lisp但是没有那么多括号（实现上确实就是S-Expr），语言中提供了充分的语法糖来让一个宏在代码里出现的更加干净。值得一提的是虽然Julia继承了Lisp但是这种继承是克制的。Julia的宏的作用主要是标记出那一部分的代码（AST）被修改了，而不是像template或者C的宏一样作为代码的一部分参与编译，例如。

template< typename T>class A {     T *a;
     A();
     // 一些代码     void myfunc(T *c);};

这个T如果不通过特别的命名规则，在后面的代码里就变得不那么明显是一个模板参数了，或者说模板也许并不是设计来进行元编程的（虽然想要实现expression template好像这个是必须要使用的）。而Julia中，宏的作用是标记

@mymacro a + b
@mymacro begin     a + b + c     c[I] + b[i]end

这个mymacro的实现一般来说会直接调用另外某些作用在一个Expr对象上的方法

macro mymacro(ex::Expr)     implementation(ex)end

这是因为概念上宏只是一个入口让你/编译器分清楚谁的代码被修改了（这很重要，explicit is better than implicit）。

我先讲一个比较简单的例子，求一个级数，因为级数往往有一些系数。比如在C里面我们可能会选择定义一个数组（为了简单，我就随便写了哈，想了想泰勒展开还要求导）

double constants[10] = {2, 4, 6, 8, 10, 12, 14, 16, 18, 20};

然后真正计算级数的函数就不会在运行时每次自己去算系数，不过这就意味着你能算的级数是有限的，由于C11之后有了可变参数的宏函数，你也许可以像那个用C的宏实现的Lisp一样自己搞一套，但是对于普通人来讲，这样的hack始终不合适（太复杂）。

double series(x) {    double result = 0.0;    for(int I=0; I< 10; I++)        result += constants[I] * x;    return result;}

而在Julia里面你可以用generated function

@generated function series(::Val{N}, x) where N    ex = :()    for i = 1:N        ex = :(ex + $(2i) * x)    end    exend

这样编译器会在编译时（compile time，是的这是一个有compile time的动态语言）计算系数并且展开这个for循环（当然你也可以不展开，但是对于小的循环for循环是有overhead的，有时候需要手动展开，不过很多时候编译器也会自己发现）。

那么这样的特性（元编程）有什么用呢，我要举TensorOperations.jl这个例子，我们上面讲到，张量的求和顺序，自身的性质等等是需要预处理然后才计算。而TensorOperation.jl就是这样通过利用Julia的宏进行代码生成从而在保证拥有抽象的同时，将大量抽象在真正运行时给inline掉了。就好像另外一个C++项目taco一样（不过没taco做的好，目前，这也是社区小带来的一个缺点，很多好点子搞起来很慢）。

我之前有简单对比过同样的功能在numpy里和TensorOperations.jl的性能，这是一个在凝聚态物理模拟里非常常用的计算过程（你只要知道我想表达这个benchmark在这个场景下是fair的就好了）：

numpy.einsum的实现

import numpy as np
def propagate(LHS, X, Y):    P = np.einsum('ijk, ipq', LHS, X)    Q = np.einsum('jkqp, jvrq', P, Y)    R = np.einsum('kprv, kvm', Q, X)    return R
LHS = np.random.randn(200, 10, 200)X = np.random.randn(200, 2, 200)Y = np.random.randn(10, 2, 10, 2)
%timeit propagate(LHS, X, Y)

numpy.tensordot的实现（调用底层C实现）：

def propagate(LHS, X, Y):    P = np.tensordot(LHS, X, axes=([0, ], [0, ]))    Q = np.tensordot(P, Y, axes=([0, 2], [0, 1]))    R = np.tensordot(Q, X, axes=([0, 3], [0, 1]))    return R
LHS = np.random.randn(200, 10, 200)X = np.random.randn(200, 2, 200)Y = np.random.randn(10, 2, 10, 2)
%timeit propagate(LHS, X, Y)

Julia的TensorOperations.jl实现：

using TensorOperationsusing BenchmarkTools
function propagate(LHS, X, Y)    @tensor R[6,7,8] := LHS[1,2,3]*X[1,4,6]*Y[2,5,7,4]*X[3,5,8]end
BLAS.set_num_threads(1)LHS = randn(200, 10, 200); X = randn(200, 2, 200); Y = randn(10, 2, 10, 2);
@benchmark propagate(LHS, X, Y)

三种实现都采用OpenBLAS作为backend，numpy.einsum是326ms，numpy.tensordot是45.1ms，Julia是24.59ms。为什么更快呢？因为求和顺序之类的预处理工作再编译的时候就给inline掉了。

而且因为可以修改AST，实际上Julia直接使用了爱因斯坦求和记号，更加符合数学表达式。

正是因为Julia拥有强大的宏，使得Julia很适合去开发DSL。例如在Julia社区中有一个很著名的优化库（实际上应该几乎所有做优化的人都知道这个吧）：JuMP.jl 早期有linda hua这样的大牛加入，JuMP本身就是一套使用DSL来方便编写优化任务和调用其它优化库的package。

借助llvm，Julia应该是目前除了CUDA以外唯一一个可以给GPU编程的动态语言。这和调CUDA库不一样。

至于调库大家都是标配，几乎每个语言都有cublas，cudnn这些的wrapper，没什么好说的。

JIT没什么新的，谁都可以用。但是不是谁都好用。觉得numba和Cython特好用的话，可以试试自己用它们写一个稀疏矩阵库看看谁的代码量少呗。但是值得一提的是目前numba一些benchmark要比Julia好一点点（差别不是很大，但也是好）。比起Cython，Julia是一个更加高级的语言，而非简单的底层。

Julia由于语言设计的上的优势，不仅能用JIT的技术，还能用一些AOT的技术来优化。这里就不懂了，从discourse上看到开发者讲的。不过可以给大家看看Julia的编译过程，图来自于CUDAnative的论文，自己重新画的。

Julia可以调用MPI，但是自己也可以被部署在服务器上，它是为云计算设计的。有自带的工具链（比如SharedArray这种）来在集群上进行并行计算，而集群上的并行计算接口和单节点的多进程是类似的。此外还有很多针对并行的语法糖（这个从MATLAB学的），比如这个会将sin函数广播到整个矩阵。

sin.(A)

广播也支持多个参数等

some_operation.(1, rand(10), rand(10))

彻底解决两语言问题是比较困难的，甚至几乎不可能。但是有没有可能在有一定性能需求，和能够容忍语言中加入一些新概念（比如强类型）的前提下，使用一个语言进行开发呢？Julia就是这样一种尝试，Julia能够为你提供接近Python的简单语法（接口），但是也能够提供充分的优化空间（不是什么代码都快，这取决于实现，但是不会说有因为语言本身导致的快不起来，纯Julia实现可以达到/接近C++的速度）。

那么这里就说下Python+ C++的问题吧。一个是两种语言，工程的复杂性增加，然后性能还有可能因为Python层面的逻辑过多而有所下降。此外从C++接入Python即便有了Pybind11之后，也是需要写额外的胶水代码的（用C++写）。那么这个代码量会很大。以至于像PyTorch这样的工程专门实现了一个基于yaml的代码生成器来生成wrapper（新版的PyTorch不是了）。

Julia的所有对象都可以直接dump到HDF5格式的文件里去

（JLD.jl），此外还有依赖更少，性能更好的纯Julia实现： JLD2.jl。然后语言本身的MIME接口也决定了dump这件事情更加容易和清楚。

光听我给您讲没人意思，看看别人怎么说：

这是Julia的开发者之一 Kristoffer Carlsson 的观点：

The thing I wish people understood about Python performance is that the difficulties come from Python’s extremely rich object model, not from anything about its dynamic scopes or dynamic types. The problem is that every operation in Python will typically have multiple points at which the user can override the behavior, and these features are used, often very extensively. Some examples are inspecting the locals of a frame after the frame has exited, mutating functions in-place, or even something as banal as overriding isinstance. "Sure you can make Python fast in some cases where you use a small subset of the language. However, Python as a whole is extremely difficult to make fast because of the incredible freedom you have. The “compiler” can assume very little, so it is difficult to generate efficient code. Just consider the incredible amount of time and money that has been spent on trying to make it faster and the relatively mediocre success it has had (Pyston, Numba, PyPy etc).

这是另外一个开发者Rdeits 的感受：

Another piece of the “why Julia” puzzle is the fact that Julia does not treat a specific set of built-in types differently. That is, I can make my own type or immutable that can do anything that Julia’s built-in types can do, and my type can be just as fast as the one provided with Julia (if I were as skilled at numerical computing as the devs, which I am not). This isn’t just a nice feature; it’s essential to making things like automatic differentiation and other numeric types fast and easy to use. It’s also why we have really nice abstractions like the idea that an image is not a matrix of doubles or bytes but instead a matrix of colors. That works because Julia lets us define a whole family of color-related datatypes which provide useful behavior and abstraction at little or no run-time cost. That’s been the biggest benefit of Julia, in my experience.

以上说了很多优点，也并不是没有缺点。我讲这些缺点，有一些是会被改进的但仍需时间，有一些根本不是Julia的出发点，我希望当你入手这门语言之后不要期待过高，合理地估计自己的需求去选择技术才是正确的。

• 社区小，一个新诞生的语言，没有大公司资金上的扶持。缺乏工业应用很难火起来，这也就意味着你遇到了问题，或者需要什么轮子需要自己解决。但是我们也不能对一个下个月才1.0的语言期望过高。这也是我坚持去搞中文社区的原因之一，社区壮大起来对所有的Julian都有好处。

• Web支持比较烂，这个没什么可洗的，用Julia搞个爬虫或者后端不是不行，但是比起flask之类的还是算了吧。不过顺便一提，Julia在IO上的interface有MIME这种东西，还是很喜欢的。

• module比较烂，0.7之前不能load当前目录的module，所有的module必须和其它package混一起或者手动添加一个搜索路径才行。然后还得写一堆include，手动include代码，不能像Python一样直接把文件当module用。这一点将会在1.1得到改善。

• CLI应用目前完全没有，因为有JIT的overhead，编译器的启动是有一点延迟的（实际上已经从0.3的0.5s左右优化到现在的0.2s左右了），如果使用CLI开发，可能是个噩梦，而我写Python一般都是不用IDE的，直接编辑器+Terminal就可以撸了。所以IDE对于Julia在最近一两年还比较重要。不过等到后面JIT的overhead能在需要的时候被去掉，这个问题会自然解决（我之前问过Julia Team为什么包管理器不提供像npm/pip这样的命令行工具的时候回答我说 ,Will Pkg3.jl support CLI? · Issue #168 · JuliaLang/Pkg.jl ）。我目前推荐Atom的Juno插件作为开发工具，这个目前bug少，更稳定。当然由千里冰封开发的Jetbrain插件也很好用，而且如果你英文不好可以用这个哈。

• 能写出好代码的人少。这个和社区小有关，和写出好的Julia代码的学习难度也有关。好上手和能写出好的工程代码是有区别的，Julia确实上手和Python一样容易，虽然学习曲线平滑，但是要学的应该相对较多。一个是编程方式和Python/C++由很大不同，这对很多人来说就会产生不适应，觉得这个语言设计的不好，遂放弃的想法，于是流行不起来。前段时间和我一起写代码的Leo就经历了从Python到Julia的过程：先是觉得不好，然后三个月后现在比我信仰还坚定...

• IDE支持不够好，由于有JIT的overhead，有明显的启动时间（实际上还有一些precompile）。IDE的问题实在是因为现在用户少，这个没有办法，所以只能到处忽悠别人来用了。JIT的overhead其实一直都在优化，我觉得这里还需要给Julia团队时间。

• 中文资料几乎为0。Julia在中国几乎没有用户（没人serious地用它写自己要用的package，至少GitHub上没怎么见过），这是一个很大的问题，甚至我接触的很多人在我介绍前连这个语言都不知道，不论它成不成熟，一个新的技术至少要让大家知道有才行。这也是我参与中文本地化的建设原因之一。推销一下自己写的一个中文教程：https://zhuanlan.zhihu.com/p/41802723 

没有完美的东西，以上不是唯一的缺点，坑很多，做好心理准备。（Julia有坑，Python/MATLAB/etc.都是有坑的）。

上面讲的已经很简单了，有些地方我也只是简要的介绍了一下，更多的优势还请去社区论坛看看帖子：https://discourse.julialang.org

关注公众号：拾黑（shiheibook）了解更多

[广告]赞助链接：

四季很好，只要有你，文娱排行榜：https://www.yaopaiming.com/
让资讯触达的更精准有趣：https://www.0xu.cn/