なぜ（非常に単純な）ベクトル化されたコードのオーダーは、Numpyよりも遅いのですか？

Question

ベクトル化などに関するJuliaの哲学をよく知っています.Numpyよりも10倍も遅いコードを受け入れることはできますが、なぜこのコードは、多くはがNumpyよりも遅いのですか？おそらくコードに間違いがあります。問題はループではなくベクトル化を使用することに関連しているとは思いません。

私の要件が強くないため、私はベクトル化を使用しています。さらにメモリ割り当てはそれほど信じられないようです（そしてそれはNumpyと非常に速いことが分かります）。コードも読みやすいです。

次のコードはPythonで書かれています。それは、複素平面の一部にわたって一般化された連続部分を計算する。連続部分は、2つの異なる関数によって定義されます。私はこのブログで画像をプロットするためにこのコードを使用します：https://kettenreihen.wordpress.com/

Pythonコードは以下の通りです。次ジュリアコードが同じ

def K(a, b, C): 
    s = C.shape 
    P1 = np.ones(s, dtype=np.complex) 
    P2 = np.zeros(s, dtype=np.complex) 
    Q1 = np.zeros(s, dtype=np.complex) 
    Q2 = np.ones(s, dtype=np.complex) 
    for n in range(1, 65): 
     A, B = a(C, n), b(C, n) 
     P = A*P2 + B*P1 
     Q = A*Q2 + B*Q1 
     P1, P2 = P2, P 
     Q1, Q2 = Q2, Q 
    return P2/Q2 

行う必要がありますが、それは同じことを計算するために2〜3分かかります。これは、約2秒で640x640の値を計算します。

function K(a, b, C) 
    s = size(C) 
    P1 = ones(Complex, s) 
    P2 = zeros(Complex, s) 
    Q1 = zeros(Complex, s) 
    Q2 = ones(Complex, s) 
    for n = 1:64 
     println(n) 
     A, B = a(C, n), b(C, n) 
     P = A.*P2 + B.*P1 
     Q = A.*Q2 + B.*Q1 
     P1, P2 = P2, P 
     Q1, Q2 = Q2, Q 
    end 
    return P2./Q2 
end 

Answer 1

あなたは抽象要素型で行列を割り当てている：Complexタイプは、すべての特定Complex{T}のタイプの抽象的なスーパータイプです。ここでは、Complex128 == Complex{Float64}やComplex64 == Complex{Float32}のような具体的な要素タイプの配列が必要です。おそらくNumPyのdtype=np.complexは特定の複合タイプを指しており、おそらくはComplex128に相当します。

C行列のさまざまな種類に対応する一般的なコードを書く場合は、Cを複素行列とし、同じ要素の型と形状の1と0の行列を作成し、あなただけの右の要素型と形状で行列を取得するためにCにonesとzeros関数を呼び出すことができます。

function K(a, b, C) 
    P1 = ones(C) 
    P2 = zeros(C) 
    Q1 = zeros(C) 
    Q2 = ones(C) 
    for n = 1:64 
     println(n) 
     A, B = a(C, n), b(C, n) 
     P = A.*P2 + B.*P1 
     Q = A.*Q2 + B.*Q1 
     P1, P2 = P2, P 
     Q1, Q2 = Q2, Q 
    end 
    return P2./Q2 
end 

がうまくいけば、それはパフォーマンスが向上します。 3つの行列を事前に割り振り、行列を回転しながらその場で操作を行うことで、さらにパフォーマンスを向上させることができます。しかし、便利な構文サポートがまだ安定版になっていないので、Julia 0.5ではまだ少し冗長ですが、ベクトル化されたバージョンよりもパフォーマンスが向上します。