1

私のpythonコードは次のとおりです...永遠にかかります。私が使用できる派手なトリックがいくつかあるに違いありませんか？私が分析している写真は小さく、グレースケールです...

def gaussian_probability(x,mean,standard_dev):
        termA = 1.0 / (standard_dev*np.sqrt(2.0*np.pi))
        termB = np.exp(-((x - mean)**2.0)/(2.0*(standard_dev**2.0)))
        g = (termA*termB)
        return g
def sum_of_gaussians(x):
        return sum([self.mixing_coefficients[i] * 
                    gaussian_probability(x, self.means[i], self.variances[i]**0.5) 
                    for i in range(self.num_components)])
def expectation():
        dim = self.image_matrix.shape
        rows, cols = dim[0], dim[1]
        responsibilities = []
        for i in range(self.num_components):
            gamma_k = np.zeros([rows, cols])
            for j in range(rows):
                for k in range(cols):
                    p = (self.mixing_coefficients[i] *      
                         gaussian_probability(self.image_matrix[j,k],
                                              self.means[i], 
                                              self.variances[i]**0.5))
                    gamma_k[j,k] = p / sum_of_gaussians(self.image_matrix[j,k])
            responsibilities.append(gamma_k)
        return responsibilities

最大化ステップは行列の責任配列のすべての要素をループしますが、比較的速く進むように見えるため、期待値ステップのみを含めました (ボトルネックはすべての gaussian_probability 計算ですか?)

4

1 に答える 1

4

次の 2 つのことを行うことで、計算を大幅に高速化できます。

  • 各ループ内で正規化を計算しないでください! 現在書かれているように、M 個のコンポーネントを持つ NxN 画像の場合、関連する各計算時間を計算して、アルゴリズムN * N * Mに導きます! O[N^4 M^2]代わりに、すべての要素を 1 回計算してから合計で割ると、 になりますO[N^2 M]。

  • 明示的なループではなく、派手なベクトル化を使用します。これは、コードをセットアップした方法で非常に簡単に実行できます。

基本的に、expectation関数は次のようになります。

def expectation(self):
    responsibilities = (self.mixing_coefficients[:, None, None] *
                        gaussian_probability(self.image_matrix,
                                             self.means[:, None, None],
                                             self.variances[:, None, None] ** 0.5))
    return responsibilities / responsibilities.sum(0)

完全な例を提供していないため、これを確認してベンチマークするために少し即興で行う必要がありましたが、簡単に説明します。

import numpy as np

def gaussian_probability(x,mean,standard_dev):
    termA = 1.0 / (standard_dev*np.sqrt(2.0*np.pi))
    termB = np.exp(-((x - mean)**2.0)/(2.0*(standard_dev**2.0)))
    return termA * termB

class EM(object):
    def __init__(self, N=5):
        self.image_matrix = np.random.rand(20, 20)
        self.num_components = N
        self.mixing_coefficients = 1 + np.random.rand(N)
        self.means = 10 * np.random.rand(N)
        self.variances = np.ones(N)

    def sum_of_gaussians(self, x):
        return sum([self.mixing_coefficients[i] * 
                    gaussian_probability(x, self.means[i], self.variances[i]**0.5) 
                    for i in range(self.num_components)])

    def expectation(self):
        dim = self.image_matrix.shape
        rows, cols = dim[0], dim[1]
        responsibilities = []
        for i in range(self.num_components):
            gamma_k = np.zeros([rows, cols])
            for j in range(rows):
                for k in range(cols):
                    p = (self.mixing_coefficients[i] *      
                         gaussian_probability(self.image_matrix[j,k],
                                              self.means[i], 
                                              self.variances[i]**0.5))
                    gamma_k[j,k] = p / self.sum_of_gaussians(self.image_matrix[j,k])
            responsibilities.append(gamma_k)
        return responsibilities

    def expectation_fast(self):
        responsibilities = (self.mixing_coefficients[:, None, None] *
                            gaussian_probability(self.image_matrix,
                                                 self.means[:, None, None],
                                                 self.variances[:, None, None] ** 0.5))
        return responsibilities / responsibilities.sum(0)

これで、オブジェクトをインスタンス化し、期待値ステップの 2 つの実装を比較できます。

em = EM(5)
np.allclose(em.expectation(),
            em.expectation_fast())
# True

タイミングを見ると、5 つのコンポーネントを持つ 20x20 の画像の場合、約 1000 倍高速です。

%timeit em.expectation()
10 loops, best of 3: 65.9 ms per loop

%timeit em.expectation_fast()
10000 loops, best of 3: 74.5 µs per loop

この改善は、画像のサイズとコンポーネントの数が増えるにつれて大きくなります。頑張ってください！

于 2015-11-14T13:37:33.003 に答える