配列同士の全組み合わせに対して、関数で定義した計算を通す

Question

以下のaとbのように、それぞれ配列（ベクトル）が複数格納されたデータがあり、この配列同士の全組み合わせに対し、別の関数で定義した配列演算をしたいと思っています（この場合、決定係数という統計指標の算出を使っています）。

この時、for文を二つ回さずに実行する方法を教えていただけないでしょうか？
np.meshgridを使おうとしたのですが、なかなかうまくいかず、止まっています。

import numpy as np

def calc(x,y):
    return 1 - sum((x - y) ** 2) / sum((x - y.mean()) ** 2)

a = np.array([[0, 0, 0, 0, ], [5, 5, 5, 10], [1, 2, 10, 12], [-3, 0, -2, 12]])
b = np.array([[1, 2, 3, 4, ], [8, 0, 10, 22], [5, 6, 7, 0], [-1, 0, 10, -12]])

# 以下は計算可能だが、aとbの全配列組み合わせで実施したい(a[0]〜a[3]×b[0]〜b[3]のイメージ)
calc(a[0], b[1]) 

Comments

[ppaul]

なぜ、for文をふたつ回さずに実行したいのですか。
他の書き方を使っても、性能的に早くなることはありません。
性能以外の原因があるのでしょうか？

[MagMag]

ありがとうございます。センサーからあがってくるデータに対して、この計算をリアルタイムで逐次処理しなければならないため、性能が必要です。

以下のリンクにあるように、for文を使わないと早いという文言を見て「for文使わない方がいいだろうな」と考え、meshgirdでの実装を考えていました。ご指摘のとおり、他の方法で早くならなければ、目的は達成できません。

https://sekailab.com/wp/2018/06/11/numpy-combinatorial-calculation-in-array/

＞他の書き方を使っても、性能的に早くなることはありません。
このリンクではfor文を回さず、meshgirdで実行した方が早いとありますが、そうならない、ということでしょうか？for文で処理せず、numpyで可能な限り一括計算に持ち込んで処理することで、処理性能は向上するものばかりと思っていましたが、、、。

[ppaul]

pythonのfor文が非常に遅いのは間違いありません。長さ4の二重ループを使うと、イテレータの呼び出しが25回、割り込みが5回起こるので、それだけでも膨大な時間がかかります。
これを回避するために、多次元配列の計算部分をC言語やアセンブラで書かれた高速なライブラリを呼び出すのがnumpyの役目です。
meshgridが使えるのであれば良いのですが、meshgridを使うのはほぼ確実に無理です。calc関数自体を変更せずにほかの制御方法を使った場合には性能は変わりません。
meg_の回答のように、calc関数をnumpyを使って書き直すのはかなり良い方法です。明示的な関数呼び出しだけでなく、**演算や/演算もC言語やアセンブラで書かれた関数の呼び出しになるので、かなりの高速化が期待できます。
性能が不足する場合は、numpyにリンクされているblasライブラリを有償のものに取り換えると、SIMD命令を有効利用することで性能が上がります。
それでも性能が不足するようでしたら、calc関数をpythonインタフェースを使ったC言語で作ることになりますが、ご自身で行うことはお勧めしません。

[MagMag]

承知しました。ありがとうございます！

Answer 1

次元を意識してブロードキャストさせていきましょう。
またそのためにsum()ではなくnp.sum()を用いましょう。

In [11]: a = np.array([[0, 0, 0, 0, ], [5, 5, 5, 10], [1, 2, 10, 12], [-3, 0, -2, 12]])
       : b = np.array([[1, 2, 3, 4, ], [8, 0, 10, 22], [5, 6, 7, 0], [-1, 0, 10, -12]])

In [12]: x = a[:, None]
       : y = b[None, :]

In [13]: 1 - np.sum((x - y) ** 2, -1) / np.sum((x - y.mean(-1, keepdims=True)) ** 2, -1)
Out[13]:
array([[-2.00000000e-01, -6.20000000e-01, -3.58024691e-01, -1.07888889e+02],
       [ 1.33333333e-01, -1.70666667e+00, -2.38709677e+00, -1.65424913e+00],
       [ 2.41610738e-01, -2.68456376e-02, -7.61904762e-01, -1.02251082e+00],
       [ 2.58503401e-01,  1.24700240e-01, -8.57142857e-01, -3.26509573e+00]])

Comments

[MagMag]

ありがとうございました！ただ、すみません。ブロードキャストで3次元にして計算するというのが、自分の実力不足でまだ理解できないです、、。

[meg_]

こんなにシンプルに書けるんですね！

[MagMag]

どれをベストアンサーにするのか迷ったのですが、一番高評価だったkiraraさんの案にさせていただきました。

[kirara0048]

aのどの次元とbのどの次元が組み合わさるか、絵を描いてみるなどすると感覚がつかめるかもしれません（難しければ次元や要素数を減らして試行）。

[MagMag]

わざわざアドバイスもいだたき、ありがとうございます！やってみます。

Answer 2

import numpy as np

def calc(x,y):
    X = np.tile(x, (1, y.shape[0])).reshape(-1, x.shape[1])
    Y = np.tile(y, (x.shape[0], 1))
    return 1- np.sum((X - Y) ** 2, axis=1) / np.sum((X - (np.mean(Y, axis=1)).reshape(np.mean(Y, axis=1).shape[0],-1)) ** 2, axis=1)

Comments

[MagMag]

なるほど。tileを使って中で組み合わせを作ってあげるのですね。ありがとうございました！

Answer 3

itertools.productを使ってみてはどうでしょうか。

from itertools import product
import pandas as pd
import numpy as np


def calc(x,y):
    return 1 - sum((x - y) ** 2) / sum((x - y.mean()) ** 2)


a = np.array([[0, 0, 0, 0, ], [5, 5, 5, 10], [1, 2, 10, 12], [-3, 0, -2, 12]])
b = np.array([[1, 2, 3, 4, ], [8, 0, 10, 22], [5, 6, 7, 0], [-1, 0, 10, -12]])

result = map(lambda p:calc(*p), product(a, b))

print(list(result))

Comments

[MagMag]

なるほど、itertoolを使う方法はわかりやすいですね。ありがとうございます。

[sfdust]

自分の回答は速度のことを全く考慮していないので、あまり参考にはなりません。
計算を走らせた時点で、内部的には2重forを回すのと同じ操作が発生するため激遅です。

現時点ではkirara0048さん、meg_さんのコードが圧倒的に速いと思います。

[MagMag]

わざわざ補足いただき、ありがとうございます。