C/C++を呼び出すことにより、Pythonで行列ベクトルの乗算とべき乗を高速化することができます

Question

私は現在、機械学習プロジェクトに取り組んでいます - データ行列Zとベクトルrhoを与えられました。 logistic loss functionはrhoです。計算には基本的な行列 - ベクトル乗算とlog/exp演算が含まれ、数字のオーバーフローを回避するためのトリックがあります（previous postで説明）。

私は現在、以下に示すようにNumPyを使用してPythonでこれを行っています（参考として、このコードは0.2秒で実行されます）。これはうまくいきますが、私はコード内で関数を複数回呼び出すので高速化したいと思います（そして、それは私のプロジェクトに関わる計算の90％以上を占めています）。

私はこのコードの実行時間を並列化せずに改善する方法を探しています（つまり1CPUのみ）。 Pythonで一般に公開されているパッケージを使用したり、CやC++を呼び出すことができて満足しています（これは実行時間が大幅に改善されたと聞いています）。データ行列Zの前処理もOKです。ベクトルrhoが（エントリー= 0の周りの50％）通常まばらであると列よりもはるかより多くの行が通常あることをよりよく計算するためのもので悪用される可能性がありますいくつかのもの（n_cols <= 100ほとんどの場合）

---

import time 
import numpy as np 

np.__config__.show() #make sure BLAS/LAPACK is being used 
np.random.seed(seed = 0) 

#initialize data matrix X and label vector Y 
n_rows, n_cols = 1e6, 100 
X = np.random.random(size=(n_rows, n_cols)) 
Y = np.random.randint(low=0, high=2, size=(n_rows, 1)) 
Y[Y==0] = -1 
Z = X*Y # all operations are carried out on Z 

def compute_logistic_loss_value_and_slope(rho, Z): 
    #compute the value and slope of the logistic loss function in a way that is numerically stable 
    #loss_value: (1 x 1) scalar = 1/n_rows * sum(log(1 .+ exp(-Z*rho)) 
    #loss_slope: (n_cols x 1) vector = 1/n_rows * sum(-Z*rho ./ (1+exp(-Z*rho)) 
    #see also: https://stackoverflow.com/questions/20085768/ 

    scores = Z.dot(rho) 
    pos_idx = scores > 0 
    exp_scores_pos = np.exp(-scores[pos_idx]) 
    exp_scores_neg = np.exp(scores[~pos_idx]) 

    #compute loss value 
    loss_value = np.empty_like(scores) 
    loss_value[pos_idx] = np.log(1.0 + exp_scores_pos) 
    loss_value[~pos_idx] = -scores[~pos_idx] + np.log(1.0 + exp_scores_neg) 
    loss_value = loss_value.mean() 

    #compute loss slope 
    phi_slope = np.empty_like(scores) 
    phi_slope[pos_idx] = 1.0/(1.0 + exp_scores_pos) 
    phi_slope[~pos_idx] = exp_scores_neg/(1.0 + exp_scores_neg) 
    loss_slope = Z.T.dot(phi_slope - 1.0)/Z.shape[0] 

    return loss_value, loss_slope 


#initialize a vector of integers where more than half of the entries = 0 
rho_test = np.random.randint(low=-10, high=10, size=(n_cols, 1)) 
set_to_zero = np.random.choice(range(0,n_cols), size =(np.floor(n_cols/2), 1), replace=False) 
rho_test[set_to_zero] = 0.0 

start_time = time.time() 
loss_value, loss_slope = compute_logistic_loss_value_and_slope(rho_test, Z) 
print "total runtime = %1.5f seconds" % (time.time() - start_time) 

Answer 1

BLASファミリのライブラリは、最高のパフォーマンスを得るためにすでに高度にチューニングされています。だから、C/C++のコードにリンクしようとすると、何のメリットもありません。しかし、いくつかのCPUに特別に調整されたものを含め、かなりの数があるので、さまざまなBLAS実装を試すことができます。

もう一つのことは、theano（またはGoogleのtensorflow）などの計算グラフ全体（上記の関数内のすべての演算）を表すライブラリを使用して、グローバル最適化を適用することです。それは、C++を介してそのグラフからCPUコードを生成することができます（また、GPUコードでも簡単なスイッチを反転させることによって）。また、自動的に記号派生派生を計算することもできます。私は機械学習の問題にtheanoを使用しました。それは本当に素晴らしいライブラリですが、学習するのが最も簡単です。

（それはコメントのためにあまりにも長いですので、私は答えとしてこれを掲示しています）

編集：私は実際にtheanoでこれを行くがあったが、結果は実際にはおよそ2倍である

CPUの速度が遅い理由については、以下を参照してください。 numpyのは非常に最適化されている

import theano 

def make_graph(rho, Z): 
    scores = theano.tensor.dot(Z, rho) 

    # this is very inefficient... it calculates everything twice and 
    # then picks one of them depending on scores being positive or not. 
    # not sure how to express this in theano in a more efficient way 
    pos = theano.tensor.log(1 + theano.tensor.exp(-scores)) 
    neg = theano.tensor.log(scores + theano.tensor.exp(scores)) 
    loss_value = theano.tensor.switch(scores > 0, pos, neg) 
    loss_value = loss_value.mean() 

    # however computing the derivative is a real joy now: 
    loss_slope = theano.tensor.grad(loss_value, rho) 

    return loss_value, loss_slope 

sym_rho = theano.tensor.col('rho') 
sym_Z = theano.tensor.matrix('Z') 
sym_loss_value, sym_loss_slope = make_graph(sym_rho, sym_Z) 

compute_logistic_loss_value_and_slope = theano.function(
     inputs=[sym_rho, sym_Z], 
     outputs=[sym_loss_value, sym_loss_slope] 
     ) 

# use function compute_logistic_loss_value_and_slope() as in original code 

Answer 2

（これはオリジナルのポストからのコードとの完全な部分的にしかコード、です）：私はとにかくそれをここに投稿します、多分それは、より良い何かをする誰か他の人のための出発点です。あなたができることは、ランダムに初期化された同じサイズのデータ​​（0に初期化されていない）を持つ他のライブラリを試し、独自のベンチマークを行うことです。

試してみるには、もちろんBLASを試してみてください。 eigenにも試してみてください。個人的に私のアプリケーションのほうが速く見つかっています。