Julia：高效易用的數值計算/優化編程語言
一、Julia中的數值計算入門
首先，Julia的LLVM-based just-in-time (JIT)編譯器在提供較低難度的編寫環境（和C相比）的情況下，仍能在數值計算計算時間的表現上和C達到類似的運算速度（見下圖，詳見Julia Micro-Benchmarks）。測試的operation包括 
例子1：簡單隨機游走（Simple Random Walk）的hitting time
walk(L)函數定義了一個從 
function walk(L=50) t = 2 x = 2 while x != 1 if rand() < 0.5 x += 1 else x -= 1 end if x > L x = L end t += 1 end return tendfunction run(L, N) times = Int64[] for i in 1:N push!(times, walk(L)) end return timesendfunction repeat(L, N, M) [mean(run(L, N)) for i in 1:M]end
我們的蒙特卡洛模擬返回的平均hitting time是99.90814250000001，并輸出了相應的pmf的histogram。
data = repeat(50, 1000, 10000);using Plotsmean(data)histogram(data)
那么我們可以看到pmf是比較類似高斯分布的，且mean集中在 
因為這個random walk也可以看成一個馬爾可夫鏈，所以我們也可以通過狀態轉移矩陣的乘積去計算它的stationary distribution，好吧這里主要是為了展示更多Julia的命令和比較一下array comprehension和（稀疏）矩陣乘法的速度。我們利用BenchmarkTools的package比較這三種算法的速度。讓 
# Array comprehensionsfunction evolve(p) L = length(p) p_new = vec(zeros(L,1)) p_new[1] = p[1] + p[2]/2 p_new[2] = p[3] / 2 for i in 3:L-1 p_new[i] = (p[i-1] + p[i+1]) / 2 end p_new[L] = (p[L-1] + p[L]) / 2 return p_newend@benchmark begin p = [i == 2 ? 1//1 : 0//1 for i in 1:L]; P = zeros(t_max,L) P[1,:] = evolve(p) for i = 2:t_max P[i,:] = evolve(P[i-1,:]) endend
# Dense matrix@benchmark begin p = [i == 2 ? 1//1 : 0//1 for i in 1:L]; P = zeros(t_max,L) Trans = zeros(L,L) Trans[1,1] = 1 Trans[1,2] = 0.5 Trans[2,3] = 0.5 for i = 3:L-1 Trans[i,i-1] = 0.5 Trans[i,i+1] = 0.5 end Trans[L,L-1] = 0.5 Trans[L,L] = 0.5 P[1,:] = Trans * vec(p) for i = 2:t_max P[i,:] = Trans * vec(P[i-1,:]) endend
# Sparse matrix@benchmark begin t_max = 1000 p = [i == 2 ? 1//1 : 0//1 for i in 1:L]; P = zeros(t_max,L) dl = 0.5 * vec(ones(L-1,1)) dl[1] = 0 d = vec(zeros(L,1)) d[1] = 1 d[L] = 0.5 du = 0.5 * vec(ones(L-1,1)) Trans = Tridiagonal(dl,d,du) P[1,:] = Trans * vec(p) for i = 2:t_max P[i,:] = Trans * vec(P[i-1,:]) endend
注意到在稀疏矩陣那里我們利用了三對角矩陣(Tridiagonal)的數據結構，而benchmarkTools都做了幾百次的sample然后輸出該代碼塊運算時間的各種統計量。我們于是知道稀疏矩陣快于array comprehension快于naive的稠密矩陣。在數值計算中，使用最有效的數據結構和代碼命令永遠是很重要的。至于如果想驗證（sanity check）這三個代碼塊是否給出了一致的輸出結果，可以使用@test_approx_eq命令，這里不多贅述。
例子2：稀疏矩陣，自定義數據結構，分布式數據結構（初步）
我們在上例中看到了一類特殊的稀疏矩陣——三對角矩陣類型的應用。這里我們指出更一般的稀疏矩陣類型在Julia中名為sparse()。一個簡單例子如下：
G = sparse(vec([1;2;3;1;1;2;3]),vec([1;2;3;2;3;1;1]),vec([1.;3;5;2;2;2;2]))
我們看到這種類型的矩陣便只會存儲非0的entry和其對應的值，這和Matlab等都是很類似的。接下來我們嘗試自定義新的矩陣類型（數據結構），我們叫他們SymArrow矩陣，一類對稱矩陣，除了對角線和第一行/列其它entry都是0的稀疏矩陣。
type SymArrowFloat
Diag::Array{Float64,1}
Row::Array{Float64,1}endimport Base: full Base.full(A::SymArrowFloat) = full(sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ))import Base: promote_rule, convert, showBase.show(io::IO, A::SymArrowFloat)=print(io,sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ) )import Base.++(A::SymArrowFloat,B::SymArrowFloat) = SymArrowFloat(A.Diag+B.Diag, A.Row+B.Row)import Base.**(A::SymArrowFloat,v::Array{Float64,1}) = [dot(vec([A.Diag[1];A.Row]),vec(v));[v[1]*A.Row[i-1]+v[i]*A.Diag[i] for i=2:size(vec(v),1)]]
我們給這類新的數據結構自定義了full()和show()函數，前者是將一個稀疏矩陣變成稠密矩陣，后者則是顯示該數據類型。同樣，我們也定義了該數據類型的加減和乘法。
# Function full worksG = SymArrowFloat([1.;3;5],[2;2.])full(G)# Function show worksG# + worksG+G# * worksG*rand(3)
這里不貼運行結果了，大家可以自行嘗試都是work的。這邊我們指出，如果想specify有理數，復數為值域，也可以單獨定義（其實主要只是為了展示Julia的數域類型233）。
type SymArrow{T}
Diag::Array{T,1}
Row::Array{T,1} endimport Base: full Base.full{T}(A::SymArrow{T}) = full(sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ))import Base.++{T}(A::SymArrow{T},B::SymArrow{T}) = SymArrow{T}(A.Diag+B.Diag, A.Row+B.Row)import Base.**{T}(A::SymArrow{T},v::Array{T,1}) = [dot(vec([A.Diag[1];A.Row]),vec(v));[v[1]*A.Row[i-1]+v[i]*A.Diag[i] for i=2:size(vec(v),1)]]import Base: promote_rule, convert, showBase.show{T}(io::IO, A::SymArrow{T})=print(io,sparse(round(Int64,[1:size(A.Diag,1);vec(ones(size(A.Row)));2:size(A.Diag,1)]),round(Int64,[1:size(A.Diag,1);2:size(A.Diag,1);vec(ones(size(A.Row)))]),[A.Diag;A.Row;A.Row] ) )
我們看到和前面相比我們就是特地把{T}加到了和Array相關的Type定義里。這里T可以是整數，有理數，復數，BigFloat等各種數域類型。
G = SymArrow{Rational}([1.;3;5],[2;2.])
G = SymArrow{Complex}([1.+0im;3;5],[2;2.-0im])
G = SymArrow{BigFloat}([1.;3;5],[2;2.])G*[BigFloat(1);2;3]
如果想省事的話，我們也可以用AbstractMatrix類型，不過這種情況下的運算速度會稍微不如之前的定義方式（因為這時候無法利用稀疏性。大家可以自行驗證，之前的那些操作在這里都可以自動work）。
type SymArrow2{T} <: AbstractMatrix{T}
Diag::Array{T,1}
Row::Array{T,1} endimport Base: size, getindexsize{T}(A::SymArrow2{T}) = (length(A.Diag), length(A.Diag))function getindex{T}(A::SymArrow2{T}, i, j)
if i == j
return A.Diag[i]
elseif i == 1
return A.Row[j-1]
elseif j == 1
return A.Row[i-1]
else
return zero(T) # otherwise return zero of type T
endend
這個例子最后我們看一下Julia支持的分布式數據結構（這里只簡單介紹array），和相關的一些分布式運算。至于真正要用Julia做并行計算不在這篇文章里重點談，可能以后可以再專門科普。我們這里只涉及DistributedArrays.jl pacakge。例子是對一個array分布式地累加求和（cumulative sum）。我們用4個核來分布式地完成操作。
using DistributedArraysnprocs() == 1 && addprocs(4)xd = @DArray [i for i = 0:39][DistributedArrays.chunk(xd, i) for i = 1:4]
我們可以看到array被自動partition成4部分分配到了相應的區域里。為了并行運算，我們需要定義prefix（）函數，主要是map-reduce框架里map的部分。
@everywhere function prefix!(op, v, v0 = 0) v[1] += v0 for i = 2:length(v) v[i] = op(v[i-1], v[i]) end return vend@everywhere function prefixlag!(op, v, v0 = 0) vi = v[1] v[1] = v0 for i = 2:length(v) tmp = op(v[i - 1], vi) vi = v[i] v[i] = tmp end return vend
我們于是便可以分布式地得到cumsum。
yd1 = map_localparts(t -> [sum(t)], xd)yd2 = prefixlag!(+, Array(yd1))yd3 = map_localparts((t,s) -> prefix!(+, t, s[1]), xd, distribute(yd2))
例子3：奇異值分解（SVD）
這個例子是我這個回答里的源碼。覃含章：SVD 降維體現在什么地方？
using Images, Colors, ImageCore, ImageViewimg = load("Albert Einstein-rare-pics56.jpg")img_RGB = float(Array(channelview(img)))R_pixel = img_RGB[1,:,:];G_pixel = img_RGB[2,:,:];B_pixel = img_RGB[3,:,:];U1,S1,V1 = svd(R_pixel);U2,S2,V2 = svd(G_pixel);U3,S3,V3 = svd(B_pixel);U1[:,1:100]*diagm(S1[1:100])*V1[1:100,:]# Full SVD: perfect reconstructionImages.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))# 25% SVD reconstruction: near perfectS1[101:400] .= 0S2[101:400] .= 0S3[101:400] .= 0Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))# 2.5% SVD reconstruction: obviously blurredS1[11:400] .= 0S2[11:400] .= 0S3[11:400] .= 0Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))# 1% SVD reconstruction: heavily blurredS1[5:400] .= 0S2[5:400] .= 0S3[5:400] .= 0Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))# only 2 nonzero entries: stripesS1[3:400] .= 0S2[3:400] .= 0S3[3:400] .= 0Images.colorim(cat(3,(U1*diagm(S1)*V1')',(U2*diagm(S2)*V2')',(U3*diagm(S3)*V3')'))
二、JuMP入門：線性規劃，錐規劃，混合整數規劃
除了數值計算，Julia語言對于優化工作者很方便的一點是，類似于Matlab中的YALMIP(https://yalmip.github.io/)，提供了high level的優化模型的輸入接口，并調用不同的solver進行優化模型的求解。在Julia中，這主要是通過JuMP這個package來實現的。
目前支持的solver一覽：（詳見Optimization packages for the Julia language）
在建模速度上，嚴格好于其它市面上的AML(algebraic modeling language)接口，且不輸于商業AML！
例子1：支持向量機（Support Vector Machine, SVM）【線性規劃】
這個例子里，我給出SVM利用線性規劃（linear programming）在Julia中的實現方式，linear kernel和RBF kernel兩種kernel function。具體來說，這里考慮的是soft-margin的SVM的對偶表示,相應對slack variable的cost是
using PyPlot, JuMPfunction linear_K(X)
K = zeros(size(X,1),size(X,1))
for i = 1:size(X,1)
for j = 1:size(X,1)
K[i,j] = dot(vec(X[i,:]), vec(X[j,:]))
end
end
return Kendfunction RBF_K(X,γ)
K = zeros(size(X,1),size(X,1))
for i = 1:size(X,1)
for j = 1:size(X,1)
K[i,j] = exp( -γ*sum(( X[i,:]-X[j,:] ).^2) )
end
end
return Kendfunction SVM_K(X,Y,C,K_type,γ=1)
n = length(Y)
M_dual = Model(solver = GurobiSolver(OutputFlag=0))
@variable(M_dual, 0<=α[1:n]<=C)
@constraint(M_dual, sum(Y.*α)==0)
if K_type == "linear"
@objective(M_dual, Max, -0.5*sum(α'*linear_K(X)*α) + sum(α))
elseif K_type == "RBF"
@objective(M_dual, Max, -0.5*sum(α'*RBF_K(X,γ)*α) + sum(α))
else
error("Wrong kernel type!")
end
solve(M_dual)
# return optimal α and margin
return getvalue(α), 1/(getobjectivevalue(M_dual)-sum(getvalue(α)))*(-2)end
注意到在這里的JuMP model里我調用了Gurobi作為solver。
例子2：計算圖的stable number【半正定規劃】
給定一個無向圖G=(V,E)，一個stable set指的是節點V的一個子集，其中沒有任何兩個節點是相連的。stable number指的就是G上所有stable set里最大的節點數量。這個問題可以relax成一個半正定規劃(semidefinite programming)問題：
Julia代碼如下。核心就是利用@SDconstraint這個macro，以及注意這里我調用了Mosek作為solver。這里有一些具體算例：Copositive Programming簡介
using JuMP, Mosekfunction StableSDP(A) n = size(A,1) StableDD = Model(solver = MosekSolver(MSK_IPAR_LOG = 0)) @variable(StableDD, X[1:n,1:n]) @SDconstraint(StableDD, X>=0) @constraint(StableDD, X.>=0) @constraint(StableDD, sum(sum(A.*X)) == 0 ) @constraint(StableDD, sum(X[i,i] for i=1:n) == 1) @objective(StableDD, Max, sum(sum(X))) solve(StableDD) return getobjectivevalue(StableDD)end
例子3：Gurobi進階【混合整數規劃】
在這里，我們指出使用JuMP的接口我們仍能具有比較自由地和Gurobi進行互動的自由度。比如說，我們可以加入lazy constraint（即一開始只specify一部分constraint，隨著求解過程慢慢加入其它constraint）或者user cut，這都可以穿插在branch and bound過程中。以及Gurobi也支持warmstart，即user可以自己提供一個初始值（比如用某些heuristic得到的，加速branch and bound）；這方面的詳細教程可看Solver Callbacks。
同樣，也可以利用Gurobi做column generation（和lazy constraint類似，只不過一開始只加入一部分variable,隨著求解過程慢慢加入其它variable）。然而不幸的是，Gurobi一直還沒有支持branch and price（即在整數規劃中做column generation），因此目前CG只能用在連續優化問題當中。這方面，Chiwei Yan- JuliaTutorial 的cutting stock problem教程是很好的學習資料。
三、交互性編程
我們指出，類似于IPython，Julia中的IJulia package也可以讓Julia的所有編譯過程在Jupyter notebook里進行。我個人是很喜歡在這個環境里進行Julia編程的（詳見：JuliaLang/IJulia.jl）。當然，我知道也有不少人喜歡Juno的：Juno, the Interactive Development Environment 可能這個更有碼農的感覺吧hhh
進一步的，我們可以有很多交互式的操作，這在Julia中主要通過Interact package實現。比如，我們可以自定義slider，按鈕等對一個參數曲線進行互動。一個例子見如下視頻，或者JuliaGizmos/Interact.jl。
至于各種畫圖，我傾向于使用Plots這個package。入門可以見：Plots - powerful convenience for visualization in Julia.
四、寫在最后
自然，本文給出的只是很少的一些例子和對Julia這門編程語言的最基本的介紹。無論你只是希望有個方便的語言調用solver，或者做數值計算居多，還是比較高級的優化算法專家，多實踐，一邊"get hands dirty"一邊學習我覺得總是最有效率的。
首當其沖的是Julia官網上提供的大量學習資源：包括視頻，具體的算例，和各種pdf教程。
作為Julia cofounder之一的Prof. Edelman的數值計算課程的個人主頁也是非常不錯的參考資料：（主要是那個Github鏈接）
Modern Numerical Computing with Julia courses.csail.mit.edu
工程師必備
- 項目客服
- 培訓客服
- 平臺客服
TOP
