優化過程對比：Julia 能戰勝 Python 和 R 成最終贏家嗎？

編譯自：Julia vs R vs Python: simple optimization

做者：ZJ，數據科學家，全棧工程師，信用風險模型團隊負責人。

在這篇文章中，做者經過一個簡單的似然函數優化（Maximum Likelihood Optimization）問題來對比 Julia，R 和 Python。這是一個比較小的優化問題，性能上的差別表現可能不太明顯，但解決問題的過程能很好地反應三者各自的優劣勢。

做者在撰寫本文時，對這三種語言的熟悉程度以下：

語言	實戰經驗
R	9 年
Julia	6 個月
Python	新手

Julia 佈道者 ChrisRackauckas 曾經說過：

若是你用 Julia 處理一個 10 秒內的問題，它的優點並不能體現出來。 而一旦處理的問題變複雜，須要花費比較長的時間，這時 Julia 的優點就會慢慢體現了。

有人用 Python 和 Julia 作過對比實驗。以 10⁵ 爲界點進行計算，當數值比 10⁵ 更小時 Python 比 Julia 快的。但數值大於 10⁵ 後，Julia 的速度就比 Python 快不少了。

優化問題

觀察序列 Q​1​​,Q​2​​,...,Q​n，咱們須要找到優化該似然函數的參數 μ 和 σ：

一般咱們會嘗試優化對數似然：

在統計學上，這是截斷的正態分佈的最大似然估計（MLE）。

Julia 的測試狀況

如下是做者使用 Julia 進行測試的狀況。使用 Julia 中的 Optim.jl，能夠直接使用特殊符號（symbols）做爲變量名稱，按照使用習慣，此處做者使用了希臘字母 μσ。Julia 還有一個 JuMP.jl 包用於優化問題。但 JuMP.jl 更適合用於更高級的優化問題，用在此處有點小題大作。

Julia 第一次優化

Julia 在執行第一次優化用了 7 秒，比 R 和 Python 都慢。對此，ChrisRackauckas 指出：

若是你須要解決 100 個 10 秒的優化問題，第一次執行須要花費 17 秒，接下來的優化不須要編譯，大約只須要 10 秒。所以，總運行時常爲 1007 秒。因此，當用 Julia 處理一個 10⁵ 秒的問題時，這 7 秒基本能夠忽略不記；但若是用 Julia 處理 5 秒甚至更小的問題時，這 7 秒的差別就特別明顯。

做者在下方硬編碼了在 MLE 估計中使用的 Q_t 的值：

using Distributions, Optim

# hard coded data\observations
odr=[0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12]
Q_t = quantile.(Normal(0,1), odr)

# return a function that accepts `[mu, sigma]` as parameter
function neglik_tn(Q_t)
    maxx = maximum(Q_t)
    f(μσ) = -sum(logpdf.(Truncated(Normal(μσ[1],μσ[2]), -Inf, maxx), Q_t))
    f
end

neglikfn = neglik_tn(Q_t)

# optimize!
# start searching 
@time res = optimize(neglikfn, [mean(Q_t), std(Q_t)]) # 7.5 seconds
@time res = optimize(neglikfn, [mean(Q_t), std(Q_t)]) # 0.000137 seconds

# the \mu and \sigma estimates
Optim.minimizer(res) # [-1.0733250637041452,0.2537450497038758] # or

# use `fieldnames(res)` to see the list of field names that can be referenced via . (dot)
res.minimizer # [-1.0733250637041452,0.2537450497038758]

輸出效果以下，排版看起來很舒服，也支持數學公示顯示：

Results of Optimization Algorithm
 * Algorithm: Nelder-Mead
 * Starting Point: [-1.1300664159893685,0.22269345618402703]
 * Minimizer: [-1.0733250637041452,0.2537450497038758]
 * Minimum: -1.893080e+00
 * Iterations: 28
 * Convergence: true
   *  √(Σ(yᵢ-ȳ)²)/n < 1.0e-08: true
   * Reached Maximum Number of Iterations: false
 * Objective Calls: 59

由此看出 Julia 的優點：

指數	描述
	Truncated(DN, lower, upper) 是定義截斷分佈的很是簡單的方法
	logpdf 函數適用於任何分佈式函數
	輸出結果條理清晰，可讀性強

Julia 的不足：

指數	描述
	若是隻是處理 10 秒內的簡單問題，7.5 秒的編譯時間會很煩人

R 的測試狀況

R 有一個 truncnorm 用於處理截斷正態

odr=c(0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12)
x = qnorm(odr)

library(truncnorm)
neglik_tn = function(x) {
  maxx = max(x)
  resfn = function(musigma) {
    -sum(log(dtruncnorm(x, a = -Inf, b= maxx, musigma[1], musigma[2])))
  }
  
  resfn
}

neglikfn = neglik_tn(x)

system.time(res <- optim(c(mean(x), sd(x)), neglikfn))
res

結果將輸出：

$par
[1] -1.0733354  0.2537339

$value
[1] -1.89308

$counts
function gradient 
      55       NA 

$convergence
[1] 0

$message
NULL

R 的優點：

指數	描述
	又處理截斷正太的專用包
	立刻輸出結果，編譯比 Julia 快

R 的不足：

指數	描述
	截斷正態沒有對數密度; 沒有簡單的方法來定義任意分佈的截斷分佈; 稀疏輸出

Python 的測試狀況

做者利用已有的 Python 學習經驗想出以下方案，輸入代碼：

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

# generate the data
odr=[0.10,0.20,0.15,0.22,0.15,0.10,0.08,0.09,0.12]
Q_t = norm.ppf(odr)
maxQ_t = max(Q_t)

# define a function that will return a return to optimize based on the input data
def neglik_trunc_tn(Q_t):
    maxQ_t = max(Q_t)
    def neglik_trunc_fn(musigma):
        return -sum(norm.logpdf(Q_t, musigma[0], musigma[1])) + len(Q_t)*norm.logcdf(maxQ_t, musigma[0], musigma[1])
    return neglik_trunc_fn

# the likelihood function to optimize
neglik = neglik_trunc_tn(Q_t)

# optimize!
res = minimize(neglik, [np.mean(Q_t), np.std(Q_t)])
res

輸出結果：

fun: -1.8930804441641282
 hess_inv: array([[ 0.01759589,  0.00818596],
       [ 0.00818596,  0.00937868]])
      jac: array([ -3.87430191e-07,   3.33786011e-06])
  message: 'Optimization terminated successfully.'
     nfev: 40
      nit: 6
     njev: 10
   status: 0
  success: True
        x: array([-1.07334252,  0.25373624])

Python 的優點：

指數	描述
	易於學習，各類支持很是好
	能很快輸出結果，比 Julia 編譯快

Python 的不足：

指數	描述
	輸出的可讀性有待提升

 綜上所述，三種的綜合對好比下：

語言	優點	不足
Julia	易於使用；完美支持截斷正太分佈；可讀性強	第一次運行編譯時間長
R	易於使用	可讀性對比 Julia 較差
Python	易於使用	可讀性對比 Julia 較差