27799


最近数値計算をpythonでやり始めた者です。

例えば

[code]

    A.shape #(3,5,5)
    B.shape #(3,)
    
[/code]

の形をもつA,Bについて

[code]

    numpy.array([ A[i,:,:] * B[i] for i in range(3)]).reshape(5,5)
    
[/code]

を計算したいのですが、この方法では規模が大きくなると遅いような気がしています。  
この操作をより速く動くように書く方法はないものでしょうか？


CC BY-SA 3.0
2016-07-20T13:45:21.647
27765
2016-08-25T10:36:25.630
2016-07-20T14:30:30.147
3068
17234
question
python numpy
numpyでブロードキャストできない掛け算をforを使わず書く方法


Benchmarkerを使って試してみました。サイズ10000*7*7を10000回回した結果、

  * `A*np.array([[B]]).T` が一番早く15秒‌​
  * 次が内包表記の245秒（zipを使った場合が220秒）
  * 最後が `A * np.dstack(B).T` で516秒でした。

とりあえず最も速かった `A*numpy.array([[B]])` を使おうと思います。  
ありがとうございました。


CC BY-SA 3.0
2016-07-21T14:42:41.080
27799
2016-07-21T18:00:23.797
2016-07-21T18:00:23.797
8000
17234
27765
answer


別のやり方を説明を兼ねて。

基本的に質問者様が行いたいのは要素積だと思います。numpy.arrayの要素積はnp.multiplyもしくは*で行えますが、今の場合、配列の次元が違うのが厄介です。  
np.multiplyは高次元の配列Aと一次元配列Bを掛けあわせた場合には、配列Aの最後の軸とBの要素積を行う仕様になっているようです。具体的には、C[i,j,k,l]
= A[i,j,k,l] * B[l]のような積を返します。

今の場合はC[i,j,k] = A[i,j,k] *
B[i]という積を行いたいので、そのためには以下の様な方法が可能だと思います（これ以外にも方法はあるかもしれません）。

  1. 配列Bの次元をAと揃える。具体的には、B[i]に余分な軸を付け足すことでB[i,j,k]のような配列に変形してから要素積をとる。このやり方としては、np.newaxisを使う、Noneを使うなどがあります。前出の回答のnp.array([[B]])も実質的に同じことをしているはずです。（正確には、追加した軸の次元は１なので、AとBの配列は違う次元ですが、その辺はブロードキャストされて計算できるようです）
  2. 転置をとることで、無理やりC[i,j,k] = A[i,j,k] * B[k]の形に落とす。
  3. np.einsumを使う。この関数を使うと、任意のテンソルを任意の形で積をとることができるので、単純に実装するだけでいいなら万能だと思いますが、場合によっては速度の問題が起こるかもしれません。というのも、一般には配列の積はどの順番でループをとるかで（numpyを使う場合には明示的ループは書きませんが）、速度が変わるはずですが、この関数ではループの順番を明示的に指定できないからです。
  4. 内包表記で具体的に書いてしまう。

くらいでしょうか。

[code]

    import numpy as np
    A = np.random.rand(3,5,5)
    B = np.random.rand(3)
    #Method 1
    print A*B[:,np.newaxis,np.newaxis]
    print A*B[:,None,None]
    print A * np.array([[B]]).T
    print A*np.expand_dims(np.expand_dims(B,1),1)
    #Method 2
    print (A.T * B).T
    #Method 3
    print np.einsum('ijk,i->ijk',A,B)
    
[/code]

手元で速度を比べてみましたが、A.shape = (10000, 7, 7)では、方法1,2,3はほぼ同じで、方法4だけ10倍程度遅かったです。


CC BY-SA 3.0
2016-08-25T10:15:47.727
28486
2016-08-25T10:36:25.630
2016-08-25T10:36:25.630
18595
18595
27765
answer