python 高次元の外積について

Question

pythonで100×100000の行列自身の外積を求めて100×100の行列にしたいのですが、メモリエラーとなり実行出来ません。データが大きくても実行出来る方法があれば、どなたか教えていただけないでしょうか？

追記
以下の行列の形にしたいと思っています。

nが100000、aが100で
xが2次元で100000列、cが2次元で100列のデータです。

sigma=10**np.linspace(-3.0,1.0,9)
lam=10**np.linspace(-3.0,1.0,9)
kernel_num = 100
def search_sigma_and_lambda(x, y, c, sigma, lam):
    nx = len(x);ny = len(y);n_min = min(nx,ny);kernel_num = len(c)


    score_new = np.inf;sigma_new = 0;lam_new = 0;
    for s in sigma:
        phi_x = np.array([compute_kernel_Gaussian(x,i,s) for i in c])
        phi_y = np.array([compute_kernel_Gaussian(y,i,s) for i in c])
        h = phi_x.sum(1) / nx
        H = np.multiply.outer(phi_y,phi_y.T) / ny←「ここの部分でメモリエラーになります。」
        phi_x = phi_x[np.arange(n_min)]
        phi_y = phi_y[np.arange(n_min)]
        for la in lam:
            B = H + lam * np.identity(kernel_num)
            B_inv = np.inv(B)
            B_inv_X = B_inv.dot(phi_y)
            print B_inv_X;sys.exit()

def compute_kernel_Gaussian(x,centers,sigma):
    x_c = np.square(x-centers).sum(1)
    return np.exp(-0.5*x_c/(sigma*sigma))

参考にしているのは以下のサイトで、
http://qiita.com/kenmatsu4/items/0a7a3ef71d4e8bb53da0
以下のサイトのソースコードをpythonに変えて実行したいと思っています。
https://github.com/hoxo-m/densratio/blob/master/R/uLSIF_search_sigma_and_lambda.R

Comments

[ikedas]

どういうことをしたいのかよくわかりません (行列自身の外積って何？) 式を書いて説明するなりしてもらえませんか。

[tonarino_sazana]

ご意見ありがとうございます。修正いたしましたので、宜しくお願い致します。

[ikedas]

「式」を書いていただきましたが、「説明」がないので式やそこに使われている変数の意味がわからないです……。たとえば、「100×100000」とおっしゃってましたが、そういう数にあたるような変数がどれなのかもわかりませんし。

[tonarino_sazana]

ご意見ありがとうございます。拙い説明で申し訳ございません。これでどうでしょうか。

[ikedas]

x = (x1, ..., x100000) と、c = (c1, ..., c100) という二つのベクトルから100行100列の行列 f を計算するということよろしいでしょうか。もしもそうなら、この程度でメモリが不足するというのはちょっと不思議です。ご自分で書かれたプログラムのソースコードも示していただけますか。

[hiro-k]

倍制度浮動小数の100x100000の行列なら約80MB程度です。最近のコンピュータなら難なく扱える容量です。しかし現実にメモリエラーとなっているということは、プログラムにメモリ使用量に対して優しくない書き方をしている部分があるはずです。同じ入力から同じ結果を出すプログラムであっても、処理の仕方によって必要なメモリ量は簡単に何桁も変化しえます。ですので、「やりたいこと=式」ではなく「どのように処理しているのか=プログラム」を示していただくことがアドバイスを得るためには必要でしょう。

[tonarino_sazana]

ご意見ありがとうございます。一部ですがソースコードを載せました。宜しくお願い致します

Answer 1

numpy の outer の中身(手元のubuntu16.04 の場合 /usr/lib/python2.7/dist-packages/numpy/core/numeric.py ) を見ると、

def outer(a, b, out=None):
    a = asarray(a)
    b = asarray(b)
    return multiply(a.ravel()[:, newaxis], b.ravel()[newaxis,:], out)

となってまして、 a と b がどんな行列でも、それを ravel() で1次元化して、その外積を求めることになってます。

なので、両方の入力に 10M個(100x100000)の要素を持つ行列を入力すると結果は、10Mx10Mの行列ということで 100T個の要素になりますので、メモリーは足りなくなって当然ですね。

提示された式と numpy.outer() の処理内容が異なっているようです。

Comments

[ikedas]

マニュアルしか見てませんでした。すみません。

元のRのソースでは、この部分にcrossprod(X, Y)というのを使ってます。これは X^t * Y (^tは行列の転置、*は行列の積) をあたえます。またRにもouter()があって、外積をあたえるようです (たぶんhiro-kさんご指摘のnumpyのと同様)。Rはよく知らないんですが、このへんが違うのではないかと。(解決後に訂正: 「要素毎の乗算」は「行列の積」の誤記)

[tonarino_sazana]

回答ありがとうございます。自分でも調べなおしてみたところ、外積ではなく内積でした。
以下のサイトでcrossprod(X,Y)がクロス積だと書かれており、勝手に外積だと勘違いしておりました。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/20.html
内積にしてみたところ、メモリエラーにならず計算出来ました。
本当にありがとうございました。

Answer 2

numpy の使い方を誤って、想定に比べて大きなサイズの配列になっているという可能性は有りませんか？
phi_y = np.array([compute_kernel_Gaussian(y,i,s) for i in c])
の後で
print phi_y.shape
として、配列のサイズが適正なのか確認してみるのは如何でしょうか？

Comments

[tonarino_sazana]

回答ありがとうございます。配列のサイズについては確認済みです。ありがとうございます

[hiro-k]

動かしている python は 32bit版でしょうか？それとも64bit版でしょうか？

http://stmind.hatenablog.com/entry/20120407/1333815923

[tonarino_sazana]

動かしているのは64bit版です

[hiro-k]

64bit版でメモリーが足りなくなるって不思議ですね・・・

[hiro-k]

a=np.random.random((100,100000))
b=np.outer(a,a.T)
でもメモリーエラーになりました。numpy にはこのサイズの外積は扱えないのかもしれません。

[tonarino_sazana]

回答ありがとうございます。そうですか、やはりnumpyでは扱えないのかもしれませんね。調べていただいてありがとうございます。

[hiro-k]

「numpy で扱えない」は誤りであったようです。
numpy.outer([100x100000],[100x100000]) の出力は御期待の100x100 ではなく、10Mx10Mになるようです。コメントではソースコード扱いに出来ないので、別回答にて。

[tonarino_sazana]

回答ありがとうございます、調べなおしてみたところ外積ではなく内積でした。お騒がせして申し訳ございません。ありがとうございました。

Answer 3

 H = np.multiply.outer(phi_y,phi_y.T) / ny

ですが、いくつか気になる点を書きます。

• NumPyのマニュアルによればnp.multiply.outer()とnp.outer()とは違った結果になるようですが、意図して前者を使っておられるのでしょうか (すみません、不勉強で、どっちが正しいのかわかりません)。前者だと必要なメモリも計算時間も増えます。
• outer()の2番目の引数をT()で転置する必要はないのではないでしょうか。

つまり、

H = np.outer(phi_y, phi_y) / ny

ということだったりしないでしょうか。

とりあえずは、いったんカーネル数を減らした上でhとHを計算してみて、print()などを使って行列の内容を出力し、思い通りのものが計算できているか確認してはどうでしょうか。

Comments

[tonarino_sazana]

回答ありがとうございます。np.multiply.outer()の方が高次元が計算出来そうな気がしたので使っていましたが、本当はnp.outer()です。np.outer()で計算出来なかったので試していました。

[ikedas]

いや、私もどっちが正しいのかわからないんですけどね。というのは、何の計算をしているのかわかってないので (ガウスカーネルみたいにもみえるけど、違うかなあ)。

ご提示のコードの参考にした書籍、ウェブページなどがあれば、それも質問に追記していただければ、もっと的確な回答をしてくださる方があるかもしれません。

[tonarino_sazana]

ご意見ありがとうございます。参考サイトを載せました。

[ikedas]

こちら確認しました。Rの crossprod() は t(X) %*% X と同等ということなので、np.outer(phi_y, phi_y) でいいです。「cross product」、「外積」、「outer product」といった言葉の意味が曖昧で、混乱しますね。

コードの残りもそのうち見ておきます、気がついたことがあればコメントするかもしれません。

[tonarino_sazana]

回答ありがとうございます。そうなんですね、調べていただいてありがとうございます。