2017-11-17 8 views
0

私はAとBの2つの行列を持っており、ここで与えられたmin-plus積を計算したい:Min-plus matrix multiplication。そのために私は以下を実装しました:PythonでMin-plus行列の乗算を高速化する方法は?

def min_plus_product(A,B): 
    B = np.transpose(B) 
    Y = np.zeros((len(B),len(A))) 
    for i in range(len(B)): 
     Y[i] = (A + B[i]).min(1) 
    return np.transpose(Y) 

これはうまくいきますが、大きなマトリックスでは遅いですが、それを速くする方法はありますか?私は、Cでの実装やGPUの使用が良い選択肢かもしれないと聞いてきました。

+1

あなたの入力はどれくらいですか?あなたの 'min_plus_product'の速度は、' dot 'を使った通常の行列乗算の速度とどうやって比較されますか? – user2357112

+0

* "...大きな行列のために遅い" * @ user2357112と同じ質問:あなたの行列の典型的なサイズは? –

+0

私はA 10000で100、B 100は20000でテストし、鉱山は135、np.dotは1.95を得ました。私の行列はこれよりもさらに大きくなることがあります。 – Rael

答えて

2

中間ディメンションが十分に大きく、エントリが均一に分散されている場合、ビットを節約するアルゴです。これは、最小の合計が典型的に2つの小さな用語からなるという事実を利用する。

import numpy as np 

def min_plus_product(A,B): 
    B = np.transpose(B) 
    Y = np.zeros((len(B),len(A))) 
    for i in range(len(B)): 
     Y[i] = (A + B[i]).min(1) 
    return np.transpose(Y) 


def min_plus_product_opt(A,B, chop=None): 
    if chop is None: 
     # not sure this is optimal 
     chop = int(np.ceil(np.sqrt(A.shape[1]))) 
    B = np.transpose(B) 
    Amin = A.min(1) 
    Y = np.zeros((len(B),len(A))) 
    for i in range(len(B)): 
     o = np.argsort(B[i]) 
     Y[i] = (A[:, o[:chop]] + B[i, o[:chop]]).min(1) 
     if chop < len(o): 
      idx = np.where(Amin + B[i, o[chop]] < Y[i])[0] 
      for j in range(chop, len(o), chop): 
       if len(idx) == 0: 
        break 
       x, y = np.ix_(idx, o[j : j + chop]) 
       slmin = (A[x, y] + B[i, o[j : j + chop]]).min(1) 
       slmin = np.minimum(Y[i, idx], slmin) 
       Y[i, idx] = slmin 
       nidx = np.where(Amin[idx] + B[i, o[j + chop]] < Y[i, idx])[0] 
       idx = idx[nidx] 
    return np.transpose(Y) 

A = np.random.random(size=(1000,1000)) 
B = np.random.random(size=(1000,2000)) 

print(np.allclose(min_plus_product(A,B), min_plus_product_opt(A,B))) 

import time 
t = time.time();min_plus_product(A,B);print('naive {}sec'.format(time.time()-t)) 
t = time.time();min_plus_product_opt(A,B);print('opt {}sec'.format(time.time()-t)) 

サンプル出力:

True 
naive 7.794037580490112sec 
opt 1.65810227394104sec 
+0

@ラエル受け入れてくれてありがとう。あなたがまだCに行くことを考えている場合のための発言です(私の経験では、numpyと痛みを伴わずにインターフェイスし、バージョン番号よりも成熟しています)。その場合、コードを大幅に単純化できます。基本的には、チャンクする必要はなく、マスクや再インデックスもしないで、 'a'の行をループし、各行の短絡を一度に行います。 –

0

可能簡単な経路がnumbaを使用することです。 1000x1000 Aに

from numba import autojit 
import numpy as np 
@autojit(nopython=True) 
def min_plus_product(A,B): 
    n = A.shape[0] 
    C = np.zeros((n,n)) 
    for i in range(n): 
     for j in range(n): 
      minimum = A[i,0]+B[0,j] 
      for k in range(1,n): 
       minimum = min(A[i,k]+B[k,j],minimum) 
      C[i,j] = minimum 
    return C 

タイミング、B行列は、次のとおり

1ループ、3の最良のオリジナルコードのループ当たり4.28秒

1ループ、3の最もよい:ループ当たり2.32秒numbaコード

関連する問題