优化过程 PK :Julia 能打败 Python 和 R 笑到最后吗?

简介:

在这篇文章中,作者通过一个简单的似然函数优化(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 快很多了。

优化问题

观察序列 Q1,Q2,...,Qn,我们需要找到优化该似然函数的参数 μ 和 σ:

2b68c46382a00a7db098f11c5b807774de5.jpg

通常我们会尝试优化对数似然:

d033feafed7dd0786f2fcc44ba165b17671.jpg

在统计学上,这是截断的正态分布的最大似然估计(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\observationsodr=[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 parameterfunction 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 estimatesOptim.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.073

输出效果如下,排版看起来很舒服,也支持数学公示显示:

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 的优势:

指数 描述
0cc31a5336d5ecdaf4de9555f5529c42fcb.jpg Truncated(DN, lower, upper) 是定义截断分布的非常简单的方法
5b2fb20e42fe30f3e082a256079d6f1085b.jpg logpdf 函数适用于任何分布式函数
3d8798365c5e7f8cafe90fe173aacfefed7.jpg 输出结果条理清晰,可读性强

Julia 的不足:

指数 描述
5096833d909abb6a79de48578e1cb7eb6bf.jpg 如果只是处理 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 的优势:

指数 描述
e040978b7b8cad6b3375de7db4d5c50e987.jpg 又处理截断正态的专用包
b3d328c6b5c8b08734cfb0e9318d5a72761.jpg 马上输出结果,编译比 Julia 快

R 的不足:

指数 描述
07382ef2c1df957e02b70f4d9cf4b9dfa43.jpg 截断正态没有对数密度; 没有简单的方法来定义任意分布的截断分布; 稀疏输出

Python 的测试情况

作者利用已有的 Python 学习经验想出如下方案,输入代码:

import numpy as npfrom scipy.optimize import minimizefrom scipy.stats import norm# generate the dataodr=[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 datadef 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 optimizeneglik = 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 的优势:

指数 描述
1c9940b0847180b2f6e7c0582c81c641556.jpg 易于学习,各种支持非常好
e2b16b29eab80baa2a99a1150cabb2c3192.jpg 能很快输出结果,比 Julia 编译快

Python 的不足:

指数 描述
e2b16b29eab80baa2a99a1150cabb2c3192.jpg 输出的可读性有待提高

 综上所述,三种的综合对比如下:

语言 优势 不足
Julia

易于使用;完美支持截断正态分布;可读性强

第一次运行编译时间长
R 易于使用 可读性对比 Julia 较差
Python 易于使用 可读性对比 Julia 较差

相关链接

本文来自云栖社区合作伙伴“开源中国”
本文作者:王练
相关文章
|
1天前
|
数据库 开发者 Python
实战指南:用Python协程与异步函数优化高性能Web应用
【7月更文挑战第15天】Python的协程与异步函数优化Web性能,通过非阻塞I/O提升并发处理能力。使用aiohttp库构建异步服务器,示例代码展示如何处理GET请求。异步处理减少资源消耗,提高响应速度和吞吐量,适用于高并发场景。掌握这项技术对提升Web应用性能至关重要。
17 10
|
2天前
|
存储 缓存 算法
如何优化Python代码?
【7月更文挑战第14天】如何优化Python代码?
14 6
|
5天前
|
机器学习/深度学习 数据采集 算法
Python实现ISSA融合反向学习与Levy飞行策略的改进麻雀优化算法优化支持向量机回归模型(SVR算法)项目实战
Python实现ISSA融合反向学习与Levy飞行策略的改进麻雀优化算法优化支持向量机回归模型(SVR算法)项目实战
|
4天前
|
网络协议 Python
python对tcp协议栈进行优化之一
**TCP优化摘要:** - MSS优化涉及调整TCP最大段大小,Python中可使用`socket.getsockopt()`查询MSS。 - Scapy是Python库,用于创建和发送网络包,可用于测试和优化协议栈性能。 - LwIP是轻量级TCP/IP协议栈,适合嵌入式设备,可通过分析和调整提升性能,特别是实时性和资源管理。
|
2天前
|
存储 大数据 数据处理
优化Python中的数据处理效率:使用生成器提升性能
在Python编程中,有效的数据处理是提升性能和效率的关键。本文将探讨如何利用生成器(generator)优化数据处理过程,通过实例展示生成器如何在内存效率和执行速度上带来显著提升。
|
4天前
|
搜索推荐 Python
快速排序的 Python 实践:从原理到优化,打造你的排序利器!
【7月更文挑战第12天】Python的快速排序**以分治策略实现高效排序,平均时间复杂度$O(nlogn)$,优于$O(n^2)$的冒泡排序。基本实现通过选取基准元素分割数组,然后递归排序两部分。优化版使用随机基准避免最坏情况。对比显示优化后排序更稳定,适应不同数据集,提升程序性能。
18 4
|
5天前
|
机器学习/深度学习 数据采集 算法
Python实现WOA智能鲸鱼优化算法优化支持向量机分类模型(SVC算法)项目实战
Python实现WOA智能鲸鱼优化算法优化支持向量机分类模型(SVC算法)项目实战
|
5天前
|
机器学习/深度学习 数据采集 算法
Python实现ISSA融合反向学习与Levy飞行策略的改进麻雀优化算法优化支持向量机分类模型(SVC算法)项目实战
Python实现ISSA融合反向学习与Levy飞行策略的改进麻雀优化算法优化支持向量机分类模型(SVC算法)项目实战
|
1天前
|
数据可视化 算法 数据挖掘
在水利系统工程中,系统工程的方法论同样适用,并可以通过Python等编程语言来实现相关的分析和优化。
在水利系统工程中,系统工程的方法论同样适用,并可以通过Python等编程语言来实现相关的分析和优化。
8 0
|
1天前
|
机器学习/深度学习 开发框架 数据可视化
我们可以从系统工程的角度来讨论如何优化组织架构,并给出一些可能涉及的Python应用领域的示例。
我们可以从系统工程的角度来讨论如何优化组织架构,并给出一些可能涉及的Python应用领域的示例。
17 0