Pythonより遅いJuliaの高速化

現在Juliaにてガウス過程のハイパーパラメーター推定(対数尤度関数の最大化)を分散的に行うプログラムを作成しています。pythonで同様のコードを書いていましたが高速化を図ろうとJuliaに変えてみました。しかし、pythonで書いた場合に比べてなぜか10倍ほど遅いプログラムになってしまいました。＠timeにていろいろと見てみたところ関数"Grad_L()"が遅いようです。調べていく過程で以下の記事を見つけエージェントを一つの構造体にまとめてみましたがそれでも早くなりませんでした。
どうにか高速化する方法はないでしょうか。構造体と関数をpythonのクラスのようにまとめた方が高速に動くのでしょうか？

ご回答よろしくお願いいたします。

https://qiita.com/Keyskey/items/a1dbaef973879dd7c3ef

Juliaでのコード

est_hp.jl
1
2using LinearAlgebra
3using DelimitedFiles
4using Distributions
5using Plots
6using Random
7
8mutable struct Agents_Struct
9    dx::Array{Float64,2}
10    dy::Array{Float64,2}
11    theta::Array{Float64,2}
12    y::Array{Float64,2}
13    z::Array{Float64,2}
14    pi_til::Array{Float64,2}
15    y_R::Array{Float64,3}
16    z_R::Array{Float64,3}
17    grad_L::Array{Float64,2}
18end
19
20function Grad_L(theta, dx, dy)#各エージェントのカーネル計算から勾配までの計算
21    grad_L = zeros(Float64,10,3)
22    grad_K3 = Matrix(1.0I, 100, 100)
23    for i in 1:10
24        x_diff = dx[i,:] * ones(1, 100) - (dx[i,:] * ones(1, 100))'#dxの差（100×100の対象行列）
25        K = theta[i,1]  .* exp.((-(x_diff) .^ 2) ./ theta[i,2]) + Matrix(theta[i,3]I, 100, 100)
26        A = -(x_diff) .^ 2 ./ theta[i,2]
27        K_Inv = inv(K)
28        grad_K1 = exp.(A)
29        grad_K2 = theta[i,1] .* grad_K1 .* ((x_diff) ./ theta[i,2]) .^ 2
30        M = K_Inv * dy[i,:]
31        grad_L[i,1] = tr(K_Inv * grad_K1) - M' * grad_K1 * M
32        grad_L[i,2] = tr(K_Inv * grad_K2) - M' * grad_K2 * M
33        grad_L[i,3] = tr(K_Inv * grad_K3) - M' * grad_K3 * M
34    end
35    return grad_L
36end
37
38function Step1(theta, grad_L, pi_til, k)
39    alpha::Float64 = 1/(k+1000)
40    tau::Int64 = 20
41    z=zeros(Float64,10,3)
42    for i in 1:10
43        z[i,:] = theta[i,:]  + alpha .* (clamp.(theta[i,:] - (grad_L[i,:] + pi_til[i,:]) ./ tau, 1e-5, 100) - theta[i,:])
44    end
45    return z
46end
47
48function Step3(z, grad_L, pi_til, dx, dy, z_R, y_R,w, N=10)
49    theta = similar(z)
50    y = similar(z)
51    pi_til = similar(z)
52    new_grad_L = zeros(Float64,10,3)
53    for i in 1:10
54        theta[i,:] = z[i,:] - z_R[i,:,i] + z_R[i,:,:] * w[i,:]
55        new_grad_L! = Grad_L(theta,dx,dy)
56        y[i,:] = new_grad_L[i,:] - grad_L[i,:] + y_R[i,:,:] * w[i,:]
57        pi_til[i,:] = 10.0 .* y[i,:] - new_grad_L[i,:]
58    end
59    theta, new_grad_L, y, pi_til
60end
61
62function ite(Agents)
63    iteration::Int64 = 10 #反復回数を設定(本来なら30万回ほど)
64    send_number::Int64 = 0
65    weight = [
66        0.4 0.2 0.2 0.2 0 0 0 0 0 0
67        0.2 0.4 0.2 0 0.2 0 0 0 0 0
68        0.2 0.2 0.2 0.2 0.2 0 0 0 0 0
69        0.2 0 0.2 0.4 0 0.2 0 0 0 0
70        0 0.2 0.2 0 0.4 0 0.2 0 0 0
71        0 0 0 0.2 0 0.4 0 0.2 0.2 0
72        0 0 0 0 0.2 0 0.4 0.2 0 0.2
73        0 0 0 0 0 0.2 0.2 0.2 0.2 0.2
74        0 0 0 0 0 0.2 0 0.2 0.4 0.2
75        0 0 0 0 0 0 0.2 0.2 0.2 0.4
76    ]
77    for k in 1:iteration
78        """step1:代替関数の最小化"""
79        Agents.z = Step1(Agents.theta, Agents.grad_L, Agents.pi_til, k )
80        """Step2:情報交換"""
81        for i = 1:10
82            #if norm(Agents.z-Agents.z_R::Real=2)>=E_z || norm(Agents.y-Agents.y_R::Real=2)>=E_y || k==0
83            for j = 1:10
84                if weight[j, i] != 0
85                    Agents.z_R[j,:, i] = Agents.z[i,:]
86                    Agents.y_R[j,:, i] = Agents.y[i,:]
87                    if i != j
88                        send_number += 1
89                    end
90                end
91            end
92        end
93        """Step3:更新"""
94        Agents.theta, Agents.grad_L, Agents.y, Agents.pi_til = Step3(
95            Agents.z,
96            Agents.grad_L,
97            Agents.pi_til,
98            Agents.dx,
99            Agents.dy,
100            Agents.z_R,
101            Agents.y_R,
102            weight
103            )
104    end
105end
106
107
108"""main関数"""
109function Est()
110    init_theta = [
111        1.0 1.0 0.01
112        0.6 0.1 0.013
113        0.6 1.3 0.005
114        1.5 0.4 0.012
115        0.74 0.8 0.006
116        2.0 0.5 0.009
117        1.2 0.71 0.015
118        1.21 0.6 0.011
119        2.0 1.0 0.008
120        0.7 0.9 0.0095
121    ]
122    Random.seed!(1)
123    D_x = randn(Float64, (10, 100))
124    D_y = sin.(pi .* D_x) + rand(Normal(0, 0.1), Base.size(D_x))
125    Agents = Agents_Struct(
126            D_x,
127            D_y,
128            init_theta,
129            zeros(Float64, 10,3),
130            zeros(Float64, 10,3),
131            zeros(Float64, 10,3),
132            zeros(Float64, 10,3 ,10),
133            zeros(Float64, 10,3 ,10),
134            zeros(Float64, 10,3),
135        )
136    Agents.grad_L = Grad_L(Agents.theta, Agents.dx, Agents.dy)
137    Agents.y, Agents.pi_til = Agents.grad_L, 9.0 .* Agents.grad_L
138#反復
139    ite(Agents)
140    println(mean(Agents.theta,dims=1))
141end
142
143Est()
144

Pythonでのコード(一部抜粋)

est_hp.py
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
3import math
4import numpy as np
5import matplotlib.pyplot as plt
6import random
7import time
8
9#-----------------------------------------------------------#
10class agent:
11    """エージェントのクラス"""
12    def __init__(self,name,dx,dy,theta_,weight):
13        self.name = name
14        self.dx,self.dy = dx,dy
15        self.theta = theta_
16        self.weight = weight
17        self.y_R = np.zeros((10,3))
18        self.z_R = np.zeros((10,3))
19        self.grad_K = np.zeros((len(self.theta), len(self.dx), len(self.dx)))
20        self.grad_L = np.zeros(len(self.theta))
21    def kernel(self,x1,x2):
22        return self.theta[0]*np.exp(-(x1-x2)**2/self.theta[1]) + self.theta[2]*(x1==x2)
23    def kgrad(self, xi, xj, d):
24        """kの勾配"""
25        if d == 0:
26            return np.exp(-((xi-xj)**2)/self.theta[1])
27        elif d == 1:
28            return self.theta[0]*np.exp(-((xi-xj)**2)/self.theta[1])*((xi-xj)/self.theta[1])**2
29        elif d == 2:
30            return (xi==xj)*1.0
31    def L_grad(self):
32        self.K = self.kernel(*np.meshgrid(self.dx,self.dx))
33        self.K_Inv = np.linalg.inv(self.K)
34        """Kの勾配"""
35        self.grad_K = np.zeros((len(self.theta), len(self.dx), len(self.dx)))
36        for i in range(3):
37            self.grad_K[i,:,:] = self.kgrad(*np.meshgrid(self.dx,self.dx), i)
38        """Lの勾配"""
39        self.grad_L = np.zeros(len(self.theta))
40        M = self.K_Inv @ self.dy
41        for d in range(3):
42            self.grad_L[d] = np.trace(self.K_Inv @ self.grad_K[d,:,:]) - M.T @ self.grad_K[d,:,:] @ M
43len(self.dx)* math.log(2*math.pi)
44    def init_y_pi(self,N):
45        self.y = self.grad_L
46        self.pi = (N-1)*self.y
47    def Step1(self,alpha,tau,N,k):
48        """近似関数を解く"""
49        self.til_theta = np.clip(self.theta - (self.grad_L+self.pi)/tau,1e-5,10000)
50        self.z = self.theta + alpha*(self.til_theta - self.theta)
51#Step2
52    def receive(self,yr,zr,i):
53        self.y_R[i] = yr
54        self.z_R[i] = zr
55    def send(self):
56        return self.y,self.z
57#Step3
58    def Step3(self,N):
59        self.y_tmp = self.grad_L
60        self.theta = self.z - self.z_R[self.name,:] + self.weight @ self.z_R
61        self.L_grad()
62        self.y = self.grad_L - self.y_tmp + self.weight @ self.y_R
63        self.pi = N * self.y - self.grad_L
64
65#-----------------------------------------------------------#
66if __name__ == '__main__':
67#-----------------------------------------------------------#
68    #(省略) 
69    #---------------------------------------#
70    for m in range(6):#閾値を変えるためのループ
71        """algorithm"""
72        agents.clear()
73        send_number = 0
74　　　　#エージェントの作成
75        for i in range(N):
76            agents.append(agent(i,D_x[i*S:(i+1)*S],D_y[i*S:(i+1)*S],theta[i,:],w[i,:]))
77        sum_L = np.zeros(3)
78        for i in range(N):
79            agents[i].L_grad()
80            sum_L += agents[i].grad_L
81        a = np.clip(mean - sum_L,1e-5,10000)
82        for i in range(N):
83            #agents[i].theta = theta_list[i]
84            agents[i].L_grad()
85            agents[i].init_y_pi(N)
86        """反復"""
87        t1 = time.time()
88        for k in range(iteration):
89            if m == 0:
90                tau=10
91                alpha=1/(k+1000)
92            elif m == 1:
93                alpha=1/(k+500)
94            elif m == 2:
95                alpha=1/(k+1000)**0.9
96            elif m == 3:
97                alpha=1/(k+500)**0.9
98            elif m == 4:
99                tau=20
100                alpha=1/(k+1000)
101            elif m == 5:
102                alpha=1/(k+1000)**0.9
103            """step1:代替関数の最小化"""
104
105            for i in range(N):
106                agents[i].Step1(alpha,tau,N,k)
107
108            """step2:近傍との通信"""
109            for i in range(N):
110                for j in range(N):
111                    if w[i][j]!=0:
112                        agents[j].receive(agents[i].send()[0],agents[i].send()[1],i)
113                        if i != j:
114                            send_number += 1
115            """step3 :更新式"""
116            for i in range(N):
117                agents[i].Step3(N)
118　  #平均
119            theta_list = np.array([agents[i].theta for i in range(N)])
120            mean = np.mean(theta_list,axis=0)
121#以下plot

(簡単のため、一部コードを改変しております。)

補足情報（FW/ツールのバージョンなど）

v.1.5

行動規範の内容に同意します

回答1件

ベストアンサー

Juliaを全く知らずに30分ググっただけの知識で回答します。

まず、Juliaが速いのは、(1)型推論によりPythonのオーバーヘッドを無くすこと、および(2)簡単に並列処理ができることが理由のようです。(1)については初回実行時に型推論のオーバーヘッドがかかるため、何度も実行する関数やループ回数の多いfor文に効果があるようです。

実際、質問者様が示したリンク先のQiitaでは、society.population = 10000 回のforループをまわすコードに対して、Pythonの10倍以上の速度を達成しています。

ところが、(1)に関して質問者様のコードGrad_Lはループ10回しかまわしておらず、Juliaの高速性を生かせない構造になっています。さらには、Pythonコードにおいては、Pythonとは名ばかりのCython（c言語のPython拡張ライブラリ）でがりがりチューンしたnumpyライブラリをフル活用しています。よって、実質的にはnumpy vs ループ10回分しか最適化効果が無いJulia という戦いになっており、numpyに軍配が上がるのは、納得できる結果です。

また(2)に関しても、私は使ったことはありませんがnumpyそのものを並列動作にする方法もあるようです。
参考: NumPyをマルチスレッドで計算させる
さらにGPUを使う方法もあります。
参考: NumPy 互換 GPU 計算ライブラリ cupy
そもそも並列動作は、スレッド間の独立性を極力保ちつつスレッド間通信を設計する必要があり、一般的に難易度が高いです。ライブラリを活用して、そのあたりの難易度を少し軽減するのは有力な方法です。

以上により、質問者様のやりたいことにはJuliaは向かないか、Juliaで高速性を発揮するにはかなりJuliaの特性を知り尽くしてそれにあわせたコードに全面的に書き換えていくか、だと考えます。

投稿2020/11/15 00:18