在期望最大化算法中的奇异协方差矩阵

Question

在期望最大化算法中的奇异协方差矩阵

3

我尝试使用Python编写EM算法，用于对不同类型的图像进行聚类。我对EM算法的理解如下：

EM Algo

根据这个理解，我在Python中编写了相应代码。以下是我的代码：

import numpy as np
import sys
from scipy import stats

# μ: mean vector ∈ R^(self.m x self.n)
# Σ: covariance matrix ∈ R^(self.m x self.n x self.n)
# π: probabilities of each component

class gmm:
    def __init__(self,n_components):
        self.m = n_components
        self.π = np.random.random((self.m,))
        self.x = None
        self.Σ = None
        self.μ = None
        self.r = None
        self.n = None # length of data

    @classmethod
    def N(cls, x, μ, Σ):
        #print(Σ)
        #print('\n')
        return stats.multivariate_normal(μ, Σ).pdf(x)
    
    def E_step(self):
        for i,x_i in enumerate(self.x):
            den = 0
            for c in range(self.m):
                #print(self.Σ[c].shape)
                #print(self.μ[c].shape)
                #sys.exit()
                den+= self.π[c]*gmm.N(x_i,self.μ[c],self.Σ[c])
            #print(f'{den=}')
            for c in range(self.m):
                num = self.π[c]*gmm.N(x_i,self.μ[c],self.Σ[c])
                self.r[i,c] = num/den
            print(f'{self.r}\n')

    def M_step(self):
        m_c = np.sum(self.r, axis = 0)
        self.π = m_c/self.m
        for c in range(self.m):
            s1 = 0
            s2 = 0
            for i in range(self.n):
                s1+= self.r[i,c]*self.x[i]
                s2+= self.r[i,c]*(self.x[i]-self.μ[c]).dot(self.x[i]-self.μ[c])
            self.μ[c] = s1/m_c[c]
            self.Σ[c] = s2/m_c[c]

    def fit(self,x, iterations = 10):
        self.x = x
        self.n = x.shape[0]
        self.r = np.empty((self.n, self.m))
        self.μ = np.random.random((self.m, self.n))
        Sigma = [np.random.random((self.n, self.n)) for i in range(self.m)]
        Sigma = [0.5*(s + s.T)+5*np.eye(s.shape[0]) for s in Sigma] # A symmetric diagonally dominant matrix is PD
        #print([np.linalg.det(s) for s in Sigma])
        self.Σ = np.asarray(Sigma)
        for i in range(iterations):
            self.E_step()
            self.M_step()
    
    def params():
        return self.π, self.μ, self.Σ

if __name__ == '__main__':
    data_dim = 5   # No. of data dimensions
    data = [[]]*6
    data[0] = np.random.normal(0.5,2, size = (300,))
    data[1] = np.random.normal(3,2, size = (300,))
    data[2] = np.random.normal(-1, 0.1, size = (300,))
    data[3] = np.random.normal(2,3.14159, size = (300,))
    data[4] = np.random.normal(0,1, size = (300,))
    data[5] = np.random.normal(3.5,5, size = (300,))

    p = [0.1, 0.15, 0.22, 0.3, 0.2, 0.03]
    vals = [0,1,2,3,4,5]
    combined = []
    for i in range(data_dim):
        choice = np.random.choice(vals, p = p)
        combined.append(np.random.choice(data[choice]))

    combined = np.array(combined)

    G = gmm(n_components = 6)
    G.fit(combined)
    pi, mu, sigma = G.params()
    print(pi)
    print(mu)
    print(sigma)

现在出现了问题。在运行代码时，协方差矩阵Σ在一些迭代后变得奇异。具体而言，在特定迭代中，Σ的所有条目突然变为相同的数值。

我尝试在此情况下向Σ添加一些随机噪声，但这似乎只是拖延了不可避免的结果。任何帮助/评论都将不胜感激。先行致谢 :)

- Hara

3

请勿使用字面上的希腊字母作为变量名。虽然在 Python 3 采用了 Unicode 后这是可能的，但它被认为是不良实践，因为它使得代码在其他计算机上运行更加困难，因为它们可能没有相关的 Unicode。相反，只需将它们命名为Sigma、mu、pi等（尽管pi通常表示数字3.1415...）。 - Adriaan

@Adriaan 我理解你的观点，但这样做可以让代码更加优美、易读和易懂。 - Hara

3

除了πℼ以外，所有的Unicode字符都是不同的。这使得阅读更加容易，但维护和使用/编码则非常困难。此外，大多数人没有一种键盘可以轻松地输入这些字符... - Ander Biguri

@AnderBiguri，那么你是支持我还是反对使用Unicode字符？ - Hara

4

我觉得非常明显，很抱歉。使用它们是一个非常糟糕的主意。 - Ander Biguri

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Andrei · Accepted Answer

为了避免协方差矩阵成为奇异矩阵，您可以在矩阵的对角线上添加任意值，即：

val * np.identity(size)

这样做可以确保协方差矩阵始终是正定且具有逆矩阵。例如，sklearn 默认使用 1e-6 作为正则化的默认值。