Ajustando um modelo de mistura gaussiana usando descida de gradiente estocástico

Estou trabalhando em um modelo de aprendizado de categoria on-line que usa descida de gradiente estocástico para se ajustar a um modelo de mistura gaussiana. O modelo é baseado no modelo de aprendizado on-line usado em Toscano e McMurray (2010).

Embora a descida do gradiente pareça estar funcionando razoavelmente bem para estimar as médias e frequências / probabilidades de mistura das categorias, estou tendo problemas com a estimativa das covariâncias dos componentes da mistura. As derivadas parciais que tenho usado para a atualização da descida do gradiente vêm de Petersen & Pedersen (2008) (p. 44)

Começando com

$p(x) = \sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)$

Petersen & Pedersen fornecem a derivada parcial em relação à matriz de covariância como$\Sigma$

$\frac{\delta \ln p(x)}{\delta \Sigma_j}=\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}\frac{1}{2}[-\Sigma_j^{-1}+\Sigma_j^{-1}(x-\mu_j)(x-\mu_j)^T\Sigma_j^{-1}]$

A etapa de gradiente de descida para cada , como eu tenho isto implementado em Python é (esta é uma ligeira simplificação ea ô Σ para todos os componentes é calculada antes de realizar a atualização):$\Sigma_j$$\Delta\Sigma$

j.sigma += learning_rate*(G(x)/M(x))*0.5*(-inv(j.sigma) + inv(j.sigma).dot((x-j.mu).dot((x-j.mu).transpose())).dot(inv(j.sigma)))

$j$$\frac{\rho_j\mathcal N_x(\mu_j,\Sigma_j)}{\sum _k \rho_k \mathcal N_x(\mu_k,\Sigma_k)}$

Então, eu estou me perguntando se há algo errado com meu código (altamente provável) ou se essa é realmente uma maneira muito ruim de ajustar esse tipo de modelo ao lidar com dados com mais de duas dimensões (consulte Toscano e McMurray para algoritmos univariados). e dados bivariados que definitivamente funcionam).

referências: Toscano, JC e McMurray, B. (2010). Integração de pistas com categorias: ponderação de pistas acústicas na fala usando aprendizado não supervisionado e estatísticas distributivas. Ciência Cognitiva, 34, 434-464.

Petersen & Pederson. The Matrix Cookbook, Versão: 14 de novembro de 2008

mus[d]$\mu_j$j.sigma$\Sigma_j$G(x)/M(x)$j$$x$

• Eu esperaria que o acesso à média, à covariância e ao cálculo da posterior para todos envolvam um jou outro d, qualquer que seja a variável que representa o componente para o qual você deseja calcular o gradiente no seu código. Se você nos o que dizer je drepresentam, poderemos ser capazes de lhe dizer mais.
• Se G(x)/M(x)acessos j.Sigmapara calcular o posterior, seu código pode não calcular o que você acha que ele faz. Talvez seja melhor calcular primeiro todos os gradientes de todos os parâmetros e executar a atualização.

• A descida do gradiente estocástico geralmente não é a primeira escolha para otimizar as misturas de gaussianos. Na maioria das vezes, a maximização de expectativas (EM) é usada (ver, por exemplo, Bishop, 2007). Mesmo que você não use o EM, considere o BFGS ou o L-BFGS (implementado em scipy.optimize) antes de usar o SGD. E mesmo se você seguir o SGD, considere usar vários pontos de dados ("lotes") por vez para estimar o gradiente ou, pelo menos, incluir um termo de momento . Analisando brevemente o artigo de Toscano e McMurray, meu palpite é que eles escolheram usar o SGD porque estavam interessados ​​em modelar a aquisição de fala de uma maneira biologicamente mais plausível, em vez de obter o melhor ajuste possível e fazer isso on-line (ou seja, um dado ponto de cada vez). Se você não precisar disso, meu conselho seria usar o EM.

  (Acabei de perceber que você pediu especificamente para o aprendizado on-line , portanto, a única opção viável para você é adicionar o termo de momento para acelerar um pouco as coisas.)

• A maneira como você escolheu calcular o gradiente é bastante ineficiente, o que desacelerará ainda mais o aprendizado. Você pode não ter visto resultados razoáveis ​​porque leva uma eternidade para o algoritmo convergir para algo interessante. Aqui está uma maneira um pouco melhor para calcular o gradiente:

  sigmaInv = inv(j.sigma)
  dSigma = G(x)/M(x) * 0.5 * (-sigmaInv + numpy.sum(sigmaInv.dot(x - mus[d]) * x))

  $\Sigma_j$