病みつきエンジニアブログ

機械学習、Python、Scala、JavaScript、などなど

スパースな行列のPearson相関係数

Scipyには、ピアソン相関係数を計算するための関数、scipy.stats.pearson というものがあるのですが、残念ながらスパースな行列(scipy.sparse)には対応していません。

実際、実装を見てみると(stats.py)、

mx = x.mean()
my = y.mean()
xm, ym = x - mx, y - my

という実装があり、単純に実装すると疎行列としては効率が悪くなります。 仮にxが疎行列だったとして、x - mx は要素がほぼ -mx な密行列になってしまうからです。

これは、定義の\sum_{i=1}^{n}{\left( x_i-\bar{x} \right)\left( y_i-\bar{y} \right)} に当たりますが、式変形をすると、陽に平均を引く必要がなくなります。(Wikipediaより引用)

f:id:yamitzky:20150527150315p:plain

ということで、これをコードに直して、

import numpy as np
import scipy.sparse

def pearsons(a, b):
    if not (scipy.sparse.issparse(a) and scipy.sparse.issparse(b)):
        raise ValueError("only sparse arrays are supported")
    if a.shape != b.shape:
        raise ValueError("shape of sparse arrays must be same")
    if a.shape[0] != 1 or b.shape[0] != 1:
        raise ValueError("size of dimention 1 must be one")
    n = a.shape[1]
    a_sum = a.sum()
    b_sum = b.sum()
    nmr = n * a.multiply(b).sum() - a_sum * b_sum
    return (nmr /
            np.sqrt(n * a.multiply(a).sum() - (a_sum) ** 2) /
            np.sqrt(n * b.multiply(b).sum() - (b_sum) ** 2))

実際に検証してみると、以下のように同じ値になります。

import scipy.sparse
from scipy.stats import pearsonr

a = scipy.sparse.csr_matrix([1, 2, 1, 0, 1])
b = scipy.sparse.csr_matrix([1, 0, 1, 4, 2])

print pearsons(a, b) # -0.93250480824031368
print scipy.stats.pearsonr(a.toarray()[0], b.toarray()[0])[0] # -0.93250480824031379