Probabilistic Graphical Models – Semana 9

Maio 21, 2012

Este post apresenta as notas de aula da semana 9, a última do curso online Probabilistic Graphical Models.

Aprendizado com dados incompletos

Visão geral

Na prática as amostras com as quais trabalhamos muitas vezes são incompletas, podendo ser oriundas de falhas na aquisição dos dados ou mesmo a presença de variáveis latentes (ocultas) no modelo.

Relembrando, utilizamos variáveis latentes no modelo para simplificá-lo. Considere o exemplo abaixo:

Figura 1: Uso de variáveis latentes simplifica o modelo

No exemplo podemos ver que a introdução da variável H reduz o número de arestas no modelo e consequentemente o número de parâmetros.

Lidando com dados incompletos

Podemos considerar dois casos de falta de informação. Se ela ocorre de maneira independente do valor observado, então basta ignoramos as amostras com dados incompletos. Por outro lado, se existe uma dependência entre a falta de informação e o valor observado, então não podemos simplesmente ignorar amostras defeituosas. Nesse caso precisamos considerar o mecanismo que modela a perda desses dados.

Perda ao acaso ou MAR (Missing At Random) é quando a probabilidade da ausência de um dado é independente de seu valor. Um exemplo simples é quando obtemos amostras de lançamentos de uma moeda. Se em ocasiões aleatórias deixamos de registrar o resultado de um lançamento, então temos um MAR.

Por outro lado, se por algum motivo só deixamos de anotar lançamentos que resultaram em cara, então a probabilidade de ausência do dado é influenciada pelo valor da amostra e nesse caso não temos um MAR.

Com a presença de dados restantes, a função de verossimilhança possui múltiplos ótimos locais. Uma intuição da razão disso é que substituindo o valor de uma amostra está faltando, o valor da verossimilhança continua o mesmo, embora o ponto no espaço que gera esse valor seja diferente. Essa característica é denominada identificabilidade.

Métodos para a otimização de verossimilhança

Uma vez que a função de verossimilhança pode ter vários ótimos locais, podemos usar métodos baseados em gradientes para estimar parâmetros que resultam em bons valores de verossimilhança.

Subida de gradiente

Podemos obter o gradiente usando o resultado do seguinte teorema:

\displaystyle \dfrac{\partial \log P(D| \Theta)}{\partial \theta_{x_i, u_i}} = \dfrac{1}{\theta_{x_i|u_i}} \sum_m P(x_i, u_i | d[m], \Theta)

Onde \theta_{x_i, u_i} representa o parâmetro para a atribuição X_i = x_i e Y_i = y_i no CPD de X_i. Temos que P(x_i, u_i | D, \Theta) representa a probabilidade da atribuição X_i = x_i e Y_i = y_i dadas as amostras e os parâmetros \Theta.

Para cada cálculo do gradiente, é preciso computar essa probabilidade para cada atribuição e para cada amostra d[m]. Usando inferência em uma árvore-clique apenas uma vez para um dado d[m], conseguimos obter P(x_i, u_i | D, \Theta).

As vantagens desse método é que ele pode ser generalizado para CPDs que não são na forma tabular. As desvantagens é que é preciso garantir que os parâmetros encontrados são de fato CPDs e para acelerar a convergência pode ser preciso combinar com métodos mais avançados.

Maximização de Esperança

Como alternativa ao método usando gradientes, existe um método específico para maximizar a verossimilhança que é chamado de maximização de esperança ou EM (Expectation Maximization).

A ideia por trás do método é que estimar parâmetros com dados completos é fácil. Por outro lado, computar a distribuição dos dados incompletos é fácil (via inferência) dados os parâmetros.

O algoritmo pode ser descrito pelos seguintes passos:

  • Escolha um ponto inicial para os parâmetros
  • Itere
    • Passo-E (Esperança): completar os dados restantes fazendo amostragem com os parâmetros atuais
    • Passo-M (Maximização): estimar novos parâmetros usando os dados completos

Vamos detalhar mais os Passo-E e Passo-M:

Esperança (Passo-E). Para cada conjunto de dados d[m] calcule P(X,U | d[m], \theta^t) para todas as atribuições de (X|U). Compute as estatísticas suficientes esperadas ou ESS (Expected Sufficient Statistics) para cada atribuição x, u, dadas por

\displaystyle \bar M_{\theta^t} [x, u] = \sum_{m = 1}^M P(x, u| d[m], \theta^t)

Seja \bar x[m] (\bar u[m]) o conjunto de variáveis de x (u) que estão faltando em d[m]. Podemos concluir que

\displaystyle P(x, u| d[m], \theta^t) = P(\bar x[m], \bar u[m] | d[m], \theta)

Maximização (Passo-M). Calcule a MLE usando as estatísticas suficientes esperadas como se fossem as estatísticas suficientes reais, ou seja,

\theta_{x|u}^{t+1} = \dfrac{\bar M_{\theta^t}[x,u]}{\bar M_{\theta^t}[u]}

A vantagem deste método é que podemos aproveitar algoritmos para calcular a MLE e na prática este método obtém bons resultados rapidamente, nas primeiras iterações. A desvantagem é que nas iterações finais a convergência é lenta.

É possível mostrar que a verossimilhança L(\theta: D) melhora a cada iteração nesse método.

Aspectos práticos da maximização da esperança

Algumas outras considerações práticas do algoritmo EM:

  1. Convergência da verossimilhança não necessariamente implica na convergência dos parâmetros
  2. Executar muitas iterações do algoritmo sobre um conjunto de treinamento pode levar a overfitting
  3. O chute inicial de parâmetros pode afetar bastante a qualidade das soluções obtidas

Referências

[1] Wikipedia – Expectation Maximization

Anúncios

Probabilistic Graphical Models – Semana 8

Maio 13, 2012

Este post apresenta as notas de aula da semana 8 do curso online Probabilistic Graphical Models.

Aprendizado estruturado

Na aula passada vimos maneiras de aprender as distribuições de probabilidade a partir de dados, onde a estrutura da rede já era dada. Agora, estamos interessados em aprender a estrutura da rede.

Os método que discutiremos a seguir são baseados em pontuação (score). Definimos uma função de pontuação que avalia quão bem a estrutura representa com os dados. O objetivo é encontrar uma estrutura que maximiza essa função.

Pontuação usando verossimilhança

Podemos definir uma função de pontuação denotada por \mbox{score}_L({\cal G}: {\cal D}), como a log-verossimilhança,

\mbox{score}_L({\cal G}: {\cal D}) = \ell((\hat \theta, {\cal G}): {\cal D})

Onde \hat \theta representa a MLE para \ell({\cal G}: {\cal D})

Definimos a distribuição de probabilidade empírica, denotada por \hat P(x, y) como a frequência de um dado resultado no conjunto de dados {\cal D} = \{x[m], y[m]\}_{m=1}^M, como

\hat P(x, y) = \dfrac{M[x,y]}{M}

Onde M[x,y] é o número de instâncias em {\cal D} tal que x[m]=x, y[m]=y. Para o caso da distribuição multinomial, temos que $\hat P$ é a MLE.

Definimos agora a informação mútua de uma distribuição \hat P, denotada por \hat I_{\hat P}(X; Y), como

\hat I_{\hat P}(X; Y) = \sum_{x, y} \hat P(x, y) \log \dfrac{\hat P(x,y)}{\hat P(x) \hat P(y)}

Temos que \hat I_{\hat P}(X; Y) \ge 0, sendo 0 se e somente se X e Y são independentes.

Definimos a entropia de uma distribuição \hat P, denotada por \hat H_P(X), como

H_P(X) = - \sum_{x} P(x) \log P(x)

Podemos então escrever a função de pontuação em função dessas novas definições

\mbox{score}_L({\cal G} : {\cal D}) = M\sum_{i=1}^n I_{\hat P}(X_i; Pa_{X_i}^{\cal G}) - M\sum_{i}H_{\hat P}(X_i)

Onde Pa_{X_i}^{\cal G} são os pais de X_i na rede \cal G.

Uma limitação dessa função é que para a maioria dos casos, a rede que a maximiza é o grafo completo. Para evitar problemas com overfitting, costuma-se inserir restrições como limitar o número de pais que um vértice pode ter ou o número de parâmetros da rede.

Há também a alternativas mais flexíveis como penalizar redes complexas.

Pontuação BIC e consistência assintótica

Uma maneira para penalizar a pontuação usando log-verossimilhança é através da pontuação BIC (Bayesian Information Criterion), que é dada por

\mbox{score}_{BIC}({\cal G}: {\cal D}) = \ell(\hat \theta_{\cal G}: {\cal D}) - \dfrac{\log M}{2}\mbox{Dim}[{\cal G}]

O último termo é o fator de penalização. \mbox{Dim}[{\cal G}] é equivalente ao número de parâmetros independentes de \cal G.

Essa penalização é um tradeoff entre o fitting dos dados e a complexidade do modelo.

Comportamento assintótico. Detalhando os termos da pontuação BIC teremos

= M\sum_{i=1}^n I_{\hat P}(X_i; Pa_{X_i}^{\cal G}) - M\sum_{i}H_{\hat P}(X_i) - \dfrac{\log M}{2}\mbox{Dim}[{\cal G}]

Note que a entropia \hat H não depende de G. Além disso, a informação mútua cresce linearmente com o número de amostras, enquanto a penalidade cresce logaritmicamente. Conforme M aumenta, maior peso é dado ao fitting e menos à simplicidade do modelo.

Consistência. Seja G^* a rede que representa verdadeiramente as amostras. Pode se mostrar que se M tende a infinito, $G^*$ maximiza essa função de pontuação.

Intuitivamente, arestas sobressalentes não contribuem com a informação mútua, mas aumentam a complexidade do modelo e portanto são penalizadas. Por outro lado, arestas necessárias são adicionadas por causa do crescimento linear comparado com a penalização logarítmica.

Pontuação Bayesiana

Essa função de pontuação é derivada do teorema de Bayes, que no nosso contexto é dado por

P({\cal G} | {\cal D}) = \dfrac{P({\cal D}|{\cal G}) P({\cal G})}{P({\cal D})}

Notando que a distribuição marginal dos dados P({\cal D}) é independente de \cal G, podemos escrever a pontuação Bayesiana como o logarítmo do numerador acima, resultando em

(1) \mbox{score}_B({\cal G} : {\cal D}) = \log P({\cal D}|{\cal G}) + \log(P({\cal G}))

Vamos analisar cada um dos termos com mais detalhes. Primeiramente, definimos a verossimilhança marginal P({\cal D}|{\cal G}) como

(2) P({\cal D}|{\cal G}) = \int P({\cal D}|{\cal G},\theta_{\cal G}) P(\theta_{\cal G}|{\cal G})d\theta_{\cal G}

O primeiro termo dessa integral é a verossimilhança, enquanto o segundo termo é a distribuição a priori dos parâmetros.

Para o caso de redes Bayesianas multinomiais, a expressão acima pode ser reescrita de uma maneira mais fácil de se computar. Para simplificar, vamos definir

\displaystyle A_i^{\cal G} = \prod_{u_i \in Val(Pa_{X_i}^{\cal G})} \dfrac{\Gamma(\alpha_{X_i|u_i})}{\Gamma(\alpha_{X_i|u_i} + M[u_i])}

onde Val(Pa_{X_i}^{\cal G}) é o conjunto de todas as atribuições possíveis dos pais de X_i, \alpha_{X_i|u_i} são os parâmetros da distribuição de Dirichlet, M[u_i] é a contagem das instâncias com o valor u_i e \Gamma é a função gamma, que é uma generalização da função de Fibonacci para o domínio contínuo.

Também definimos

\displaystyle B_i^{\cal G} = \prod_{x_i^j \in Val(X_i)} \dfrac{\Gamma(\alpha_{x_i^j|u_i} + M[x_i^j, u_i])}{\Gamma(\alpha_{x_i^j|u_i})}

Onde Val(X_i) é o conjunto de todos os possíveis valores de X_i. Com isso, podemos escrever a verossimilhança marginal como

P({\cal D}|{\cal G}) = \prod_i (A_i^{\cal G} B_i^{\cal G})

Aplicando o logaritmo na expressão acima, podemos decompô-la em famílias de pontuação,

\log P({\cal D}|{\cal G}) = \sum_i \mbox{FamScore}_B(X_i|Pa_{X_i}^{\cal G} : {\cal D})

Analisada a verossimilhança marginal, vamos prosseguir com o segundo termo de (1), P({\cal G}), denominada estrutura a priori.

Podemos escolher a distribuição inicial de \cal G como sendo uniforme ou penalizando a complexidade da rede como por exemplo o número de arestas ou o número de parâmetros.

A constante normalizadora da estrutura a priori é a mesma para qualquer rede. Como aplicamos o logaritmo, essa constante é somada na função de pontuação independentemente da rede e portanto podemos ignorá-la.

Parâmetros a priori. Voltando a falar sobre a verossimilhança marginal, temos a distribuição a priori dos parâmetros, P(\theta_{\cal G} | {\cal G}). Da mesma forma que na estimativa de parâmetros, vamos utilizar a distribuição de Dirichlet (BDe prior).

Lembrando, temos o tamanho do conjunto equivalente, denotado por \alpha, que representa o número de amostras que foram hipoteticamente observadas. Temos também a rede B_0, contendo as distribuições a priori.

O número de amostras hipoteticamente observadas para cada combinação de valores é dada por \alpha[x_i, Pa_{x_i}^{\cal G}] = \alpha P(x_i, Pa_{x_i}^{\cal G} | B_0). Note que não necessariamente Pa_{x_i}^{\cal G} = Pa_{x_i}^{B_0}

A vantagem da rede Bayesiana a priori B_0 é que podemos usá-la para prover as distribuições a priori para qualquer rede \cal G.

O uso dessa distribuição para a distribuição a priori dos parâmetros tem a propriedade de que rede I-equivalentes possuem a mesma pontuação Bayesiana.

Além disso, conforme o número de instâncias amostradas M tende a infinito, uma rede com distribuição a priori de Dirichlet satisfaz a pontuação de BIC mais uma constante

\mbox{score}_{BIC}({\cal G} : {\cal D}) + c

Aprendizado de árvores

Agora que definimos algumas possibilidades de se avaliar a qualidade de uma rede através de funções de pontuação, vamos focar no problema de encontrar a rede que maximiza uma dada função de pontuação, que é essencialmente um problema de otimização.

Inicialmente, vamos nos restringir a redes que são árvores. Para esse caso, temos que cada vértice tem no máximo um pai. Seja p(i) > 0 o pai de X_i ou p(i) = 0 se X_i não possui pai. Assim, podemos re-escrever uma dada função de pontuação como

\mbox{score}({\cal G} : {\cal D}) = \sum_{i} \mbox{score}(X_i|Pa_{X_i}^{\cal G} : {\cal D})

\displaystyle = \sum_{i:p(i)>0} \mbox{score}(X_i | X_{p(i)} : {\cal D}) + \sum_{i:p(i)=0} \mbox{score}(X_i : {\cal D})

Somando \mbox{score}(X_i : {\cal D}) para i : p(i) > 0 no segundo somatório e subtraindo do primeiro temos

\displaystyle = \sum_{i:p(i)>0} (\mbox{score}(X_i | X_{p(i)} : {\cal D}) - \mbox{score}(X_i : {\cal D})) +

\displaystyle \sum_{i} \mbox{score}(X_i : {\cal D})

Agora o segundo somatório é independente da estrutura da árvore. Vamos definir pesos para as arestas

w(i \rightarrow j) = \mbox{score}(X_i | X_{p(i)} : {\cal D}) - \mbox{score}(X_i : {\cal D})

O objetivo é encontrar o conjunto de arestas que maximiza a soma desses pesos.

Se a função de pontuação for a verossimilhança, temos que w(i \rightarrow j)  = M I_{\hat P} (X_i ; X_j) e todos os pesos são não negativos e portanto a solução ótima é uma árvore (ou seja, um grafo conexo).

Para as pontuações BIC e BDe, a penalização pode fazer com que os pesos nas arestas sejam negativos. Dessa forma, a solução ótima é uma floresta (ou seja, não necessariamente conexo).

Dizemos que uma função de pontuação satisfaz a equivalência de pontuação se estruturas I-equivalentes possuem a mesma pontuação. Esse é o caso das pontuações de verossimilhança, BIC e BDe. Nesse caso, podemos mostrar que w(i \rightarrow j) = w(j \rightarrow i). Portanto, para fins de computar a pontuação de uma árvore, podemos ignorar o sentido das arestas.

Com isso, podemos reduzir esse problema para o da árvore geradora mínima, que pode ser resolvido através de algoritmos como o Kruskal e o Prim.

A entrada é um grafo completo com custo dado por w(i \rightarrow j). Uma vez computada a árvore geradora mínima desse grafo, removemos as arestas de custo negativo dessa árvore, obtendo uma floresta (conexa ou não). Essa floresta representa a estrutura que maximiza a função de pontuação.

Aprendizado em grafos gerais

Vimos que é possível determinar a estrutura que maximiza as funções de pontuação baseadas em verossimilhança, BIC e BDe de maneira exata se nos restringirmos a florestas. Entretanto, para grafos gerais, o problema é NP-difícil, como enuncia o tereoma baixo.

Teorema. Encontrar uma estrutura de pontuação máxima onde os vértices têm no máximo k pais é NP-difícil para k > 1.

Podemos então usar heurísticas tais como a busca local. Alguns tipos de vizinhança comumente utilizados são remoção de aresta, adição de aresta e inversão do sentido da aresta.

Como solução inicial podemos usar uma rede vazia (sem arestas), a melhor árvore, uma rede aleatória ou alguma outra estrutura construída a partir de conhecimento a priori.

Aqui também temos os conhecidos problemas de busca local como ótimos locais e platôs. Aliás, para pontuações que respeitam a equivalência de pontuação o problema de platôs é comum porque estruturas I-equivalentes são geralmente vizinhas com relação à reversão de arestas.

Variantes que visam contornar esses efeitos são a busca tabu e o GRASP, por exemplo.

Complexidade

Vamos analisar a complexidade de um passo da busca local. Há n(n-1) arestas que podem sofrer remoção ou inversão de sentido ou podem ser adicionadas, levando a um número O(n^2) de vizinhos. A complexidade total é dada por O(n^3 M).

Podemos usar a propriedade da “decomponibilidade” da função de pontuação para melhorar a complexidade da busca local. Lembrando, temos que,

\mbox{score}({\cal G} : {\cal D}) = \sum_i \mbox{score}(X_i | Pa_{X_i}^{\cal G} : {\cal D})

É possível mostrar que os únicos componentes do somatório que são modificados entre duas soluções vizinhas são correspondes aos vértices incidentes à aresta alterada. Portanto, a ideia é manter um conjunto de deltas representando a diferença entre o custo de uma solução vizinha e a solução atual. Cada vértice tem O(n) deltas que correspondem às possíveis arestas incidentes a ele. Portanto temos que manter O(n^2) deltas.

Ao realizarmos um movimento, no máximo O(n) deltas devem ser atualizados, correspondendo aos vértices incidentes à aresta alterada. Para atualizar um delta, precisamos re-avaliar todas as amostras, ficando com a complexidade O(M).

O conjunto de deltas pode ser mantido em uma fila de prioridade, de modo que podemos escolher o próximo movimento (o de maior delta) em O(1). Os O(n) deltas recomputados devem ser atualizados na fila, o que pode ser feito em O(n\log n).

A complexidade total de um movimento explorando essas ideias é O(nM + n\log n)

Leitura adicional

[1] Wikipedia – Bayesian Information Criterion
[2] Wikipedia – Minimum Description Length
[3] Wolfram – Gamma Function


Probabilistic Graphical Models – Semana 7

Maio 6, 2012

Este post apresenta as notas de aula da semana 7 do curso online Probabilistic Graphical Models.

Estimativa de parâmetros em Redes Bayesianas

Estimativa de verossimilhança máxima

Dada uma amostra {\cal D} = \{x[1], \cdots, x[M] \} obtida de maneira independente de uma distribuição desconhecida, queremos estimar \theta que represente \cal D.

Uma medida da qualidade dessa distribuição é a verossimilhança, denotada por L(\theta:{\cal D}) e definida como a probabilidade de se obter \cal D dada a distribuição \theta, ou seja

L(\theta:{\cal D}) = P({\cal D}|\theta) = \prod_{i=1}^{M} P(x[i]|\theta)

Definimos estatística suficiente como uma estatística que provê a mesma informação de uma amostra, sendo suficiente para calcularmos a verossimilhança. Vamos defini-la em termos de uma função que recebe uma amostra x e retorne uma tupla \mathbb{R}^k.

Um exemplo é a estatística suficiente para uma distribuição multinomial de k valores. Podemos definir a função s(x_i) = (0, \cdots, 0, 1, 0, \cdots, 0) que recebe um valor x_i e retorna uma tupla apenas com a i-ésima posição contendo 1 e o resto 0.

Agora calculamos a verossimilhança definindo \bar M = (M_1, \cdots, M_k) = \sum_{x[i] \in {\cal D}} s(x[i]). Note que M_k representa o número de amostras com o valor k. Logo, a verossimilhança é dada por

L(\bar{\theta}:{\cal D}) = \prod_{i=1}^{k} \theta_i^{M_i}

Outro exemplo é a distribuição Gaussiana com média \mu e variância \sigma^2. A distribuição pode ser dada por

p(X) = \dfrac{1}{\sqrt{2\pi \sigma^2}} \exp \left(-x^2 \dfrac{1}{2\sigma^2} + x \dfrac{\mu}{\sigma^2} - \dfrac{\mu^2}{\sigma^2}\right)

Podemos definir a estatística suficiente como s(x) = (1, x, x^2). Seja

\bar M = (M_1, M_2, M_3) = \sum_{x[i] \in {\cal D}} s(x[i])

podemos escrever a função de verossimilhança como

L(\theta: {\cal D}) = \left(\dfrac{1}{\sqrt{2\pi \sigma^2}}\right)^k \exp \left(-M_3 \dfrac{1}{2\sigma^2} + M_2 \dfrac{\mu}{\sigma^2} - M_1 \dfrac{\mu^2}{\sigma^2}\right)

O problema da estimativa de verossimilhança máxima (Maximum Likelihood Estimation, MLE) consiste em, dado um conjunto de amostras \cal D, encontrar a distribuição \theta que maximiza a verossimilhança L(\theta: D).

No caso da distribuição multinomial, podemos mostrar que o MLE é

\theta^*_i = \frac{M_i}{\sum_{i=1}^m M_i}

No caso da distribuição Gaussiana, temos

\mu^* = \dfrac{1}{M}

\sigma^* = \sqrt{\dfrac{1}{M} \sum_m (x[m] - \mu^*)^2}

Estimativa de verossimilhança máxima em Redes Bayesianas

Podemos generalizar a estimativa de verossimilhança para as distribuições em Redes Bayesianas. Suponha que tenhamos uma rede Bayesiana e queremos estimar as distribuições de probabilidades, nesse caso, os CPDs.

Novamente, dado um conjunto de amostras, \cal D, definimos a verossimilhança

L(\Theta: D) = \prod_{m} P(x[m]:\Theta)

Podemos fatorar essa distribuição em um produto dos CPDs de cada nó i usando a regra da cadeia.

= \prod_{m} \prod_i P(x_i[m] | \pmb{u}_i[m] : \theta_i)

Neste caso \pmb{u}_i representa os pais da variável i na rede. Re-arranjando os fatores, podemos trocar a ordem dos produtórios e obter,

= \prod_{i} \prod_m P(x_i[m] | \pmb{u}_i[m] : \theta_i)

Definindo a verossimilhança local como

L_i({\cal D}: \theta_i) = \prod_m P(x_i[m] | \pmb{u}_i[m] : \theta_i)

Temos finalmente que

L(\Theta: D) = \prod_{i} L_i(D:\theta_i)

Supondo inicialmente que os CPDs não são compartilhados por nenhuma variável, temos um \theta_i para cada variável. Vamos mostrar como determinar esses CPDs de maneira independente.

L_i({\cal D}: \theta_i) = \prod_{m} P(x[m]|\pmb{u}_[m] : \theta) = \prod_{m} P(x[m]|\pmb{u}_[m] : \theta_{X|\pmb{U}})

Re-arranjando os produtos para agrupá-los por valores de x e \pmb{u} (buckets), temos

\displaystyle = \prod_{x, \pmb{u}} \prod_{\substack{m:x[m] = x\\\pmb{u}[m] = \pmb{u}}} P(x[m]|\pmb{u}_[m] : \theta_{X|\pmb{U}})

Observamos que P(x|\pmb{u}) corresponde a uma entrada do CPD dado por \theta_{X|\pmb{U}}, que denotamos por \theta_{x|\pmb{u}}, portanto,

\displaystyle = \prod_{x, \pmb{u}} \prod_{\substack{m:x[m] = x\\\pmb{u}[m] = \pmb{u}}} \theta_{x|\pmb{u}}

Note que o termo \theta_{x|\pmb{u}} será multiplicado pelo número de amostras tais que x[m] = x e \pmb{u}[m] = \pmb{u}, que denotamos por M[x,\pmb{u}],

\displaystyle = \prod_{x, \pmb{u}} \theta_{x|\pmb{u}}^{M[x,\pmb{u}]}

Finalmente, vemos que a distribuição P(X|\pmb{U} = \pmb{u}) corresponde à distribuição multinomial, então podemos calcular o MLE \theta^*_{x|\pmb{u}} como

\theta^*_{x|\pmb{u}} = \dfrac{M[x,\pmb{u}]}{\sum_{x'} M[x',\pmb{u}]} = \dfrac{M[x,\pmb{u}]}{M[\pmb{u}]}

CPDs compartilhados. Note que a fórmula que obtivemos depende apenas dos valores de x e \pmb{u} e do CPD \theta_{x|\pmb{u}}. Assim, se tivermos CPDs compartilhados entre variáveis, como no modelo temporal, a única diferença é que podemos contar as ocorrências em M[x,\pmb{u}] múltiplas vezes para uma mesma amostra.

Observamos finalmente que o número de buckets cresce exponencialmente com o número de pais da variável. Isso faz com que os valores de \frac{M[x,\pmb{u}]}{M[\pmb{u}]} fiquem muito pequenos, levando à ocorrência da fragmentação e requerendo um maior número de amostras.

Estimativa Bayesiana

Uma maneira alternativa de se fazer estimativas das distribuições dos CPDs é usando uma outra rede Bayesiana. Neste caso supomos que cada amostra de \cal D representa um nó da rede e elas dependem de uma variável contínua que representa a distribuição \theta.

Neste caso a variável \theta tem um CPD correspondente a uma distribuição inicial P(\theta). A ideia é usar a rede para determinar a distribuição a posteriori de P(\theta) uma vez que observamos \cal D. Usando a regra de Bayes temos que,

(1) P(\theta|x[1], \cdots, x[M]) = \dfrac{P(x[1], \cdots, x[M]|\theta) P(\theta)}{P(x[1], \cdots, x[M])}

Veja que o primeiro termo do numerador corresponde à verossimilhança. P(\theta) é a distribuição a priori e P(\theta|x[1], \cdots, x[M]) a distribuição a posteriori. Uma vez que temos o conjunto de amostragem, o denominador é uma constante normalizadora.

Portanto, podemos reescrever (1) como

(2) P(\theta|{\cal D}) \propto P({\cal D}|\theta) P(\theta)

Distribuição multinomial

Suponha que \theta é uma distribuição multinomial sobre k valores. Como um valor a priori da distribuição de \theta, em geral utiliza-se a distribuição de Dirichlet, denotada por Dir(\alpha_1, \cdots, \alpha_k), onde \alpha_i é chamado de hiper-parâmetro. Basicamente, a distribuição é proporcional a um produtório,

Dir(\alpha_1, \cdots, \alpha_k) \propto \prod_{i=1}^{k} {\theta_i}^{\alpha_i - 1}

Intuitivamente os hiper-parâmetros correspondem às frequências de variáveis amostradas discriminadas por valor.

Usando a distribuição de Dirichlet como a distribuição a priori e lembrando que para uma distribuição multinomial a verossimilhança é dada por

P({\cal D}|\theta) = \prod_{i=1}^{k} \theta_i^{M_i}

Por (2), concluímos que a distribuição a posteriori é dada por,

P(\theta|{\cal D}) = \prod_{i=1}^{k} \theta_i^{\alpha_i + M_i - 1}

Note que P(\theta|{\cal D}) = Dir(\alpha_1 + M_1, \cdots, \alpha_k + M_k), ou seja, a distribuição a posteriori continuou sendo uma distribuição de Dirichlet. Devido a essa característica, dizemos que a distribuição de Dirichlet é um conjugado a priori para a distribuição multinomial.

Previsão Bayesiana

Considere uma rede bayesiana simples como a da Figura 1. Suponha que a variável \theta tem distribuição Dir(\alpha_1,\cdots,\alpha_k) = \prod_{j=1}^{k} \theta_j^{\alpha_j - 1}.

Figura 1: Variável condicionada por uma distribuição de Dirichlet

Podemos calcular a distribuição em X multiplicando os CPD’s e marginalizando \theta, ou seja,

(3) P(X) = \int_{\theta} P(X|\theta) P(\theta) d\theta

Para um valor específico de X = x_i, temos que P(X = x_i|\theta) = \theta_i. Substituindo também P(\theta) pela definição da eq. de Dirichlet teremos

P(X = x_i) = \frac{1}{Z} \int_{\theta} \theta_i \left(\prod_{j=1}^{k} \theta_j^{\alpha_j - 1}\right)d\theta

Que segundo as aulas é igual a

(4) = \dfrac{\alpha_i}{\sum_{j=1}^{k} \alpha_j}

Agora, suponha que tenhamos um conjunto de M+1 variáveis condicionadas por \theta, que novamente supomos possuir a distribuição Dir(\alpha_1,\cdots,\alpha_k), conforme a Figura 2.

Figura 2: Rede Bayesiana com M+1 variáveis condicionadas a uma distribuição de Dirichlet

Imagine que observamos as M primeiras variáveis (conjunto \cal D) e queremos prever a distribuição da M+1-ésima, temos que

P(X[M+1] | {\cal D}) = \int_{\theta} P(X[M+1]|{\cal D},\theta) P(\theta|{\cal D}) d\theta

Dado \theta, X[M+1] é independente de \cal D, logo

P(X[M+1] | {\cal D}) = \int_{\theta} P(X[M+1]|\theta) P(\theta|{\cal D}) d\theta

Vimos anteriormente que P(\theta|{\cal D}) = Dir(\alpha_1 + M_1, \cdots, \alpha_k + M_k). Note então que o lado direito da equação acima é idêntico ao lado direito da equação (3). Portanto, podemos calcular $P(X[M+1] = x_i |{\cal D})$, mas usando \alpha_i + M_i, levando a

(5) P(X[M+1] = x_i |{\cal D}) = \dfrac{\alpha_i + M_i}{\sum_{j=1}^{k} (\alpha_j + M_j)}

Definimos como o tamanho da amostra equivalente a soma \alpha = \alpha_1 + \cdots + \alpha_k. Observe que se \alpha_j = 0, teríamos o MLE.

Usando redes Bayesianas, damos um “chute” inicial na distribuição, que eventualmente converge para o mesmo valor que o MLE conforme o número de iterações aumenta. A vantagem é que temos uma maior robustez nas iterações iniciais, quando temos poucas amostras. O gráfico da Figura 3 mostra bem esse fenômeno.

Figura 3: Gráfico de estimativa de uma distribuição de Bernoulli, usando diferentes métodos.

Note que a magnitude de \alpha define o quanto a distribuição inicial terá efeito sobre as distribuições obtidas a priori. Veja como a Dir(.5,.5) (verde-claro) quase não tem diferença sobre a MLE, enquanto Dir(10,10) (preto) muda bastante.

Estimativas Bayesianas para Redes Bayesianas

Podemos generalizar as ideias apresentadas anteriormente para estimar as distribuições dos CPDs de toda uma rede Bayesiana. Consideremos uma rede Bayesiana com uma variável X e uma variável Y que depende condicionalmente de X.

Para estimar os CPDs dessas duas variáveis usando um conjunto de M amostras, podemos construir outra rede Bayesiana conforme a Figura 4.

Figura 4: Estimativa de parâmetros para uma rede Bayesiana”

Note que \theta_X corresponde ao CPD de X e \theta_{Y|X} corresponde ao CPD de Y. Se partirmos de distribuições iniciais para \theta_X e \theta_{Y|X}, temos que as amostras (X[m], Y[m]) são independentes.

Dado o conjunto de amostras {\cal D}, temos que \theta_X e \theta_{Y|X} são independentes. Assim, a distribuição a posteriori de \theta pode ser dada por

P(\theta_X, \theta_{Y|X} | {\cal D}) = P(\theta_X|{\cal D}) P(\theta_{Y|X}|{\cal D})

Podemos refinar o modelo e considerar os possíveis valores de X, criando um CPD de Y para cada valor. Por exemplo, se X for uma variável binária teremos a rede Bayesiana da Figura 5

Figura 5: Estimativa refinada de parâmetros para uma rede Bayesiana

Neste caso, dado o conjunto \cal D, não fica evidente pela rede que \theta_{Y|x^0} é independente de \theta_{Y|x^1}, mas aqui temos independência específica por contexto, pois quando o valor de X é conhecido, uma das duas arestas incidentes a cada Y[m] some.

De qualquer maneira, podemos computar as probabilidades a posteriori de \theta independentemente. Se a distribuição a priori de \theta_{X|\pmb{u}} for Dir(\alpha_{x^1|\pmb{u}}, \cdots, \alpha_{x^k|\pmb{u}}), a distribuição a posteriori será Dir(\alpha_{x^1|\pmb{u}} + M[x^1|\pmb{u}], \cdots, \alpha_{x^k|\pmb{u}} + M[x^k|\pmb{u}])

Uma escolha inicial dos hiper-parâmetros das distribuições de Dirichlet, consiste em escolher um tamanho de amostra equivalente \alpha e para cada CPD distribuí-lo pelo número de entradas. Por exemplo, se o CPD for sobre uma única variável binária, teremos Dir(\alpha/2, \alpha/2). Se tivermos uma variável binária dependendo de outra, fazemos a divisão de \alpha por 4.


Estimativa de parâmetros em Redes de Markov

Vamos apresentar métodos para a estimativa de parâmetros em redes de Markov ou MRF (Markov Random Fields) e para redes de Markov condicionais ou CRF (Conditional Random Fields). Ao invés de um produto de fatores, vamos considerar a representação mais genéricas através de modelos log-lineares.

Verossimilhança máxima para modelos log-lineares

Considere uma distribuição conjunta de um modelo log-linear

P(\pmb{X} : \theta) = \dfrac{1}{Z(\theta)} \tilde P(\pmb{X} : \theta) = \dfrac{1}{Z(\theta)} \exp \left\{\sum_{i=1}^{k} \theta_i f_i(\pmb{x}) \right\}

Note que para simplificar a notação, estamos supondo que f_i ignora todas as entradas de \pmb{x} que não pertencem a seu domínimo D_i. Aqui Z(\theta) é a constante normalizadora dada pela soma de \tilde P(\pmb{X}=\pmb{x}) para todos as possíveis atribuições \pmb{x}, ou seja

Z(\theta) = \sum_{\pmb{x}} \exp (\sum_i \theta_i f_i (\pmb{x}))

Podemos definir a log-verossimilhança de um conjunto de amostras \cal D, denotado por \ell(\theta:{\cal D}) como

\ell(\theta:{\cal D}) = \ln (L(\theta:{\cal D})) = \sum_m(\ln P(\pmb{x}[m] : \theta))

Substituindo P(\pmb{x}[m] : \theta) pela distribuição log-linear temos

= \sum_m (\sum_i \theta_i f_i(\pmb{x}[m]) - \ln Z(\theta))

podemos re-arrajar os somatórios e obter:

= (\sum_i \theta_i (\sum_m f_i(\pmb{x}[m]))) - M\ln Z(\theta)

Diferentemente do caso de redes Bayesianas, aqui não podemos decompor em verossimilhanças locais porque o termo Z(\theta) acopla os termos \theta_1, \cdots, \theta_k.

Porém, podemos mostrar que essa função é côncava, podendo ser resolvida numericamente por um método de descida de gradiente, por exemplo. Na prática, costuma-se utilizar o L-BFGS (Low Memory Broyden Fletcher Goldfarb Shanno).

Para computar os gradientes, vamos considerar a derivada da log-verossimilhança.

\dfrac{\partial}{\partial\theta_i} \dfrac{1}{M} \ell(\theta:{\cal D}) =  \dfrac{\partial}{\partial\theta_i} \left((\sum_i \theta_i \dfrac{1}{M} (\sum_m f_i(\pmb{x}[m]))) - \ln Z(\theta)\right)

Note que dividimos a log-verossimilhança pela constante M antes de derivar. Isso é para simplificar o resultado. Como estamos derivando com relacão a \theta_i, o único termo que restará do primeiro somatório é

\dfrac{1}{M} \sum_m f_i(x[m])

que corresponde à média dos valores de f_i considerando todas as amostras. Definimos esse termo como esperança empírica, denotada por E_{\cal D}[f_i(\pmb{X})]. Veja que este valor independe da distribuição \theta.

Já para o segundo termo, aplicamos inicialmente a regra da cadeia para obter

\dfrac{\partial}{\partial\theta_i} \ln Z(\theta)
= \dfrac{1}{Z(\theta)} \sum_{\pmb{x}} \left( \dfrac{\partial}{\partial \theta_i} \exp \{\sum_j \theta_j f_j(\pmb{x}) \} \right)

Podemos aplicar novamente a regra da cadeia para \exp, obtendo

= \dfrac{1}{Z(\theta)} \sum_{\pmb{x}} \left( f_i(\pmb{x}) \exp \{\sum_j \theta_j f_j(\pmb{x}) \} \right)

Rearranjando os termos ficamos com

= \sum_{\pmb{x}} \left( \dfrac{1}{Z(\theta)} \exp \{\sum_j \theta_j f_j(\pmb{x}) \} \right) f_i(\pmb{x})

Veja que o priimeiro termo do somatório é exatamente nossa distribuição original, P(\pmb{X} : \theta) com \pmb{X} = \pmb{x}. Reescrevendo temos:

\dfrac{\partial}{\partial\theta_i} \ln Z(\theta) = \sum_{\pmb{x}} P(\pmb{X} = \pmb{x}) f_i(\pmb{x})

Esse termo representa a esperança de f_i considerando a distribuição \theta e será denotado por E_\theta[f_i]. Veja que para calcular a probabilidade de cada possível atribuição \pmb{x}, podemos usar os métodos de inferência vistos em aulas passadas.

Resumindo, temos o seguinte resultado:

\dfrac{\partial}{\partial\theta_i} \dfrac{1}{M} \ell(\theta:{\cal D}) = E_{\cal D}[f_i(\pmb{X})] - E_\theta[f_i]

Usando essa diferença para construir o gradiente, obtemos finalmente o seguinte teorema:

Teorema. A distribuição \hat \theta resulta na MLE se e somente se E_{\cal D}[f_i(X)] = E_\theta[f_i] para todo i=1,\cdots,k

MLE para CRFs

No caso de CRFs (Conditional Random Fields), temos a distribuição conjunta condicional

P_\theta(\pmb{Y}|\pmb{x}) = \dfrac{1}{Z_{\pmb{x}}(\theta)}\tilde P_\theta(\pmb{x}, \pmb{Y}),

onde Z_{\pmb{x}}(\theta) = \sum_{\pmb{Y}} \tilde P_\theta(\pmb{x}, \pmb{Y})

Considere o conjunto de amostras {\cal D} = \{(x[m], y[m])\}_{m=1}^M. Podemos novamente definir a log-verossimilhança:

\ell_{Y|x}(\theta:{\cal D}) = \sum_{m=1}^M \ln P_\theta(y[m], x[m] : \theta)

= \sum_m (\sum_i \theta_i f_i(\pmb{x}[m]) - \ln Z_{x[m]}(\theta))

Note que neste caso não podemos extrair o termo Z_{x[m]}(\theta) para fora do somatório porque ele depende de x[m].

Para determinar o gradiente, podemos novamente calcular a derivada em relação a \theta_i:

\dfrac{\partial}{\partial \theta_i} \dfrac{1}{M} \ell_{Y|x}(\theta:{\cal D}) = \dfrac{\partial}{\partial \theta_i} \dfrac{1}{M} \sum_m (\sum_i \theta_i f_i(\pmb{x}[m]) - \ln Z_{x[m]}(\theta))

Usando as mesmas ideias do caso de MRF’s, concluiremos que

\dfrac{\partial}{\partial \theta_i} \dfrac{1}{M} \ell_{Y|X}(\theta:{\cal D}) = \dfrac{1}{M} \sum_{m=1}^{M} (f_i(x[m],y[m]) - E_\theta[f_i(x[m], Y)])

Onde E_\theta[f_i(x[m], Y)] =  \sum_Y P_\theta(Y|x[m]) f_j(Y, x[m])

A diferença aqui é que a cada passo do algoritmo, devemos recalcular P_\theta(Y|x[m])] para cada amostra x[m]. Porém, em geral o conjunto Y será menor para CRF’s.

Estimativa MAP para MRFs e CRFs

Da mesma forma que para redes Bayesianas, podemos “chutar” uma distribuição a priori e obter uma distribuição a posteriori. Entretanto, devido ao acoplamento presente na verossimilhança de redes de Markov, não há uma fórmula fechada para a distribuição, como a de Dirichlet para redes Bayesianas.

O que fazemos é manter um conjunto de parâmetros que maximizem a dstribuição a posteriori, ou seja, uma MAP (Maximum A Posteriori).

Uma possível distribuição a pirori que pode ser usada é a distribuição gaussiana com métia 0.

P(\theta : \sigma^2) = \prod_{i=1}^{k} \dfrac{1}{\sqrt{2\pi}\sigma} \exp \left\{ - \dfrac{\theta_i^2}{2\sigma^2} \right\}

Intuitivamente, quanto maior o valor de \sigma, maior a confiança de que o valor de \theta_i esteja próximo de 0.

Outra alternativa é usar a distribuição de Laplace,

P(\theta:\beta) = \prod_{i=1}^{k} \dfrac{1}{2\beta} \exp \left\{-\dfrac{|\theta_i|}{\beta} \right\}

Queremos encontrar a distribuição \theta que maximiza a distribuição conjunta dado o conjunto de amostras, ou seja,

\mbox{argmax}_\theta P({\cal D}, \theta) = \mbox{argmax}_\theta P({\cal D} | \theta) P(\theta)

Neste caso P({\cal D} | \theta) é a verossimilhança e P(\theta) é a distribuição a priori. Como o logaritmo é monotomicamente crescente com sua entrada, podemos aplicar o log ao termo dentro de argmax:

= \mbox{argmax}_\theta(\ell(\theta:{\cal D}) + \log P(\theta))

O plot de \log P(\theta) para a distribuição Gaussiana, é uma função quadrática. Esta curva corresponde a uma penalização da log-verossimilhança e é conhecida como regularização-L2.

De maneira análoga, o plot para a distribuição de Laplace é uma função linear e é equivalente à regularização-L1.


Probabilistic Graphical Models – Semana 6

Abril 29, 2012

Este post apresenta as notas de aula da semana 6 do curso online Probabilistic Graphical Models.

Inferência em modelos temporais

Tracking da crença em tempo real

Um modelo temporal pode ser representado sucintamente através de templates. Porém, a rede expandida (grounded ou unrolled) pode ficar enorme. Ficaria inviável manter uma distribuição conjunta considerando todas as variáveis do modelo.

O tracking da crença em tempo real é uma maneira aproximada de calcular a crença para um conjunto de variáveis em um dado nível t a partir apenas da crença no nível anterior.

Dada a crença \sigma^{(t)}(S^{(t)}) do conjunto de variáveis S^{(t)} no nível t temos

(1) \sigma^{(t)}(S^{(t)}) = P(S^{(t)} | o^{(1:t)})

Podemos deduzir a seguinte relação de inferência da crença em uma iteração dada a crença na iteração anterior.

\sigma^{(t + 1)}(S^{(t+1)}) \propto P(o^{(t+1)}|S^{(t+1)}) \sigma^{(\cdot t + 1)}(S^{(t+1)})

Note que precisamos renormalizar a distribuição \sigma^{(t + 1)}(S^{(t+1)}) para obter o valor correto.

Nas aulas não ficou muito claro como essa recorrência foi obtida. Complementando com respostas no fórum, decidi escrever uma versão mais mastigada como segue.

Prova:

Seja \sigma^{(\cdot t + 1)} (observe o ponto antes do t) a distribuição do estado S^{(t+1)} sem levar em conta a observação em t+1, que é dada por,

(2) \sigma^{(\cdot t + 1)}(S^{(t+1)}) = P(S^{(t+1)} | o^{(1:t)})

Usando a relação

(3) P(A|B) = \frac{P(AB)}{P(B)},

podemos mostrar que

P(S^{(t+1)}, S^{(t)} | o^{(1:t)}) = P(S^{(t+1)} | S^{(t)}, o^{(1:t)}) P( S^{(t)} | o^{(1:t)})

Marginalizando para S^{(t)}, temos que

P(S^{(t+1)} | o^{(1:t)}) = \sum_{S^{(t)}}P(S^{(t+1)}, S^{(t)} | o^{(1:t)})

ou seja,

\sigma^{(\cdot t + 1)}(S^{(t+1)}) = \sum_{S^{(t)}} P(S^{(t+1)} | S^{(t)}, o^{(1:t)}) P(S^{(t)} | o^{(1:t)})

Pela estrutura desse tipo de rede Bayesiana, temos que S^{(t+1)} é independente de tudo dado S^{(t)}, ou seja,

P(S^{(t+1)} | S^{(t)}, o^{(1:t)}) = P(S^{(t+1)} | S^{(t)})

Usando essa relação e (1) concluímos então que,

\sigma^{(\cdot t + 1)}(S^{(t+1)}) = \sum_{S^{(t)}} P(S^{(t+1)} | S^{(t)}) \sigma^{(t)}(S^{(t+1)})

Dada a relação (3), podemos mostrar que dados os conjuntos de variáveis A, B, C,

P(A|BC) = \dfrac{P(C|AB)P(A|B)}{P(C|B)}

Com isso, podemos encontrar uma relação de recorrência para \sigma^{(t + 1)}(S^{(t+1)}). Por definição temos que

\sigma^{(t + 1)}(S^{(t+1)}) = P(S^{(t+1)} | o^{(1:t)}, o^{(t+1)})

Usando (2) temos

= \dfrac{P(o^{(t+1)}|S^{(t+1)}o^{(1:t)})P(S^{(t+1)}|o^{(1:t)})}{P(o^{(t+1)}|o^{(1:t)})}

Observando que o^{(t+1)} é independente de o^{(1:t)} dado S^{(t+1)}, substituindo (2) e observando que o denominador P(o^{(t+1)}|o^{(1:t)}) é uma constante normalizadora, podemos concluir que

\sigma^{(t + 1)}(S^{(t+1)}) \propto P(o^{(t+1)}|S^{(t+1)}) \sigma^{(\cdot t + 1)}(S^{(t+1)})

Emaranhamento

Uma limitação do tracking da esperança em tempo real é que conforme expandimos a rede Bayesiana dinâmica (t aumenta), as independências condicionais para um dado nível (X(t) \perp Y(t) | Z(t)) tendem a sumir, pois cada vez mais existirão caminhos ativos entre elas (através dos ancestrais).

Assim, a crença em tempo real tende a ficar muito correlacionada, se distanciando da crença exata. Note que neste caso que a crença exata deveria considerar praticamente todas as variáveis do modelo para levar em conta as dependências, ficando com um tamanho exponencial.


Tomada de decisão

Utilidade esperada máxima

Considere um problema de decisão \cal D consistindo de um conjunto \pmb{X} de variáveis aleatórias, um conjunto de ações A e uma distribuição dos estados das variáveis condicionados por A, denotada por P(\pmb{X}|A). Considere uma função utilidade denotada por U(\pmb{X}, A), que corresponde a função peso para as atribuições de \pmb{X} e A.

Definimos a utilidade esperada do problema de decisão \cal D como a soma ponderada da utilidade de uma atribuição (\pmb{x},a), ponderada por sua probabilidade, ou seja:

EU[{\cal D}[a]] = \sum_x P(\pmb{x}| a) U(\pmb{x},a)

O problema da utilidade esperada máxima (Maximum Expected Utility (MEU)) consiste em escolher a ação a que maximiza a utilidade esperada, ou seja,

a^* = \mbox{argmax}_a EU[{\cal D}[a]]

Há casos em que a ação é condicionada por variáveis aleatórias, como no exemplo da Figura 1, onde a ação Found (relacionada a fundar ou não fundar uma empresa), está condicionada ao resultado de uma pesquisa de mercado, associada à variável Survey.

Figura 1: Exemplo

Neste caso, ao invés de encontrar uma única ação a ser tomada, devemos encontrar uma tabela onde as entradas são os possíveis valores das variáveis, a qual definimos regra de decisão, denotada por \delta_A e corresponde a essencialmente ao CPD, P(A | \mbox{predecessores}(A)).

Neste caso, a soma ponderada para calcular a utilidade esperada engloba todos os valores de a, ou seja,

EU[{\cal D}[\delta_A]] = \sum_{\pmb{x},a} P_{\delta_A}(\pmb{x}| a) U(\pmb{x}|a)

O objetivo agora é encontrar o CPD \delta_A.

No caso de uma rede Bayesiana, podemos desmembrar a distribuição conjunta P_{\delta_A}(\pmb{x}| a) em um produto de fatores. Definindo \pi_X como os pais da variável X na rede, podemos reescrever EU[{\cal D}[\delta_A]] como,

= \sum_{X_1,\cdots, X_n,A} \left((\prod_{i} P(X_i|\pi_{X_i}) \cdot U(\pi_U) \cdot \delta_A(A|\pi_A) \right)

Definindo W como \{X_1,\cdots, X_n\} \setminus \pi_A, podemos fatorar o somatório e isolar \delta_A,

= \sum_{\pi_A,A} \delta_A(A|\pi_A) \sum_W ((\prod_{i} P(X_i|\pi_{X_i}) \cdot U(\pi_U))

Denotando por \mu(A, \pi_A) o produto dos fatores marginalizando W, temos

= \sum_{\pi_A,A} \delta_A(A|\pi_A) \mu(A, \pi_A)

Com isso é podemos construir a solução ótima \delta^*_A(A|\pi_A) como um CPD determinístico, cujas entradas podem ser obtidos de maneira gulosa:

\delta^*_A(A|\pi_A) =   \begin{cases}  1 & \mbox{Se } a = \mbox{argmax}_A \mu(A, \pi_a) \\  0 & \mbox{Caso contrario}  \end{cases}

Funções de utilidade

A escolha da função utilidade é complicada, principalmente porque a utilidade é um conceito subjetivo.

Além disso, a utilidade não é uma função linear do valor esperado. Por exemplo, para a maior parte das pessoas, ganhar R$500 é melhor do que ganhar R$1000 com 50% de chance.

O gráfico da Figura 2 é uma aproximação da percepção de utilidade em função da quantia de dinheiro. Note que a função cresce menos do que a função linear.

Figura 2: Gráfico utilidade vs. recompensa

Devido a essa característica, quando consideramos risco, como no exemplo de se ganhar R$1000 com 50% de chance, temos que o valor utilidade nesse caso é menor do que R$500 (o valor esperado), mais especificamente é equivalente a ganhar ~R$400 com 100% de chance.

Um caso onde isso fica mais evidente é o paradoxo de São Petersburgo, onde o valor esperado é infinito, mas a utilidade em geral é equivalente a R$2!

Valor da informação perfeita

Uma vez que obter informações que auxiliem no suporte a decisão — como pesquisa de mercado, sensores mais precisos, exames médicos — possuem um custo (equipamentos mais caros, exames mais perigosos), devemos ser capazes de medir qual o ganho no valor utilidade eles trazem.

Considere um diagrama de decisão \cal D. Denotamos por {\cal D}_{X \rightarrow A} como o diagrama \cal D adicionado da aresta X \rightarrow A. Definimos o valor de uma observação X antes de decidir por uma ação de A como a melhoria na utilidade que essa observação traz, ou seja,

\mbox{VPI}(A|X) = \mbox{MEU}({\cal D}_{X \rightarrow A}) - \mbox{MEU}({\cal D})

Teorema. \mbox{VPI}(A|X) \ge 0, sendo que \mbox{VPI}(A|X) = 0 se e somente se, a decisão ótima para \cal D é a também a solução ótima para {\cal D}_{X \rightarrow A}.

Conclusão

Logo no final do ano passado eu havia falado que um dos objetivos ao fazer o curso era criar uma boa base para aprender técnicas de tomada de decisão sob incerteza. Foi uma surpresa boa descobrir que este assunto foi coberto nesse curso!

Em um post anterior, eu já havia indicado este post que apresenta uma discussão bem bacana sobre loterias e funções de utilidade.


Probabilistic Graphical Models – Semana 5

Abril 22, 2012

Este post apresenta as notas de aula da semana 5 do curso online Probabilistic Graphical Models.

Inferência de MAP

Continuando a falar sobre inferência de probabilidade máxima a priori (MAP). Vamos apresentar casos fáceis para os quais existem algoritmos polinomiais e então descrever uma heurística que consiste em decompor instâncias gerais em subproblemas fáceis de resolver.

Casos especiais

Alguns casos especiais onde resolver o MAP é fácil.

Correspondência entre variáveis. No caso em que as entradas do fator representam a correspondência entre dois grupos de variáveis (o primeiro grupo é o domínio do fator, o segundo são os valores), podemos reduzir o MAP para o problema de encontrar o emparelhamento bipartido de custo máximo, que pode ser resolvido eficientemente.

Fatores associativos. Redes arbitrárias, desde que as variáveis usem fatores singletons, ou seja, com domínio igual à própria variável ou potenciais supermodulares entre pares, que são fatores entre duas variáveis binárias cuja tabela é uma matriz 2×2 onde a soma dos elementos da diagonal principal é maior ou igual do que os da anti-diagonal.

Nesse caso em particular, é possível obter soluções exatas usando corte mínimo em grafos (que por sua vez pode ser encontrada através do fluxo máximo).

Fatores de cardinalidade. Fatores sobre um número arbitrário de variáveis binárias X_1, \cdots, X_k, onde o custo de uma atribuição é igual à soma das variáveis.

Fatores com padrões esparsos. Um fator sobre variáveis X_1, \cdots, X_k, com um número pequeno de entradas para as quais o custo está especificado, sendo o custo das entradas restante igual (constante).

Fatores convexos. Variáveis binárias ordenadas X_1, \cdots, X_k. Fator possui domínio X_1, \cdots, X_k, mas só estão definidas as entradas onde as variáveis com 1’s são “conexas”, ou seja se X_i = 1 e X_j = 1 para i \le j, então X_k = 1 para todo i < k < j.

Decomposição dual

Considere um conjunto de fatores divididos entre fatores singletons e não-singletons (grandes), denotados respectivamente por \theta_i(x_i) e \theta_F(\pmb{x}_F).

MAP(\theta) = \max_x \left( \sum_{i = 1}^n \theta_i (x_i) + \sum_{F} \theta_F(\pmb{x}_F) \right)

Seja \lambda_{F,i}(x_i) uma função sobre o valor da variável x_i associada ao fator F. Por exemplo, se x_i é binária, teremos dois possíveis fatores, \lambda_{F,i}(0) ou \lambda_{F,i}(1)

Observamos primeiramente que \sum_{i=1}^n \sum_{F:i\in F}  \lambda_{F,i}(x_i) = \sum_F \sum_{i \in F}  \lambda_{F,i}(x_i)

Pois ambos os somatórios iteram sobre todos as ocorrências das variáveis x_i considerando todos os elementos em F. Com isso temos

MAP(\theta) = \max_x ( \sum_{i = 1}^n (\theta_i (x_i) + \sum_{F:i\in F}  \lambda_{F,i}(x_i))
+ \sum_{F} (\theta_F(\pmb{x}_F) - \sum_{i \in F}  \lambda_{F,i}(x_i)) )

Vamos definir as variáveis

\bar \theta_i(\lambda) = \theta_i (x_i) + \sum_{F:i\in F}  \lambda_{F,i}(x_i)

\bar \theta_F(\lambda) = \theta_F(\pmb{x}_F) - \sum_{i \in F}  \lambda_{F,i}(x_i)

Seja então

L(\lambda) = \sum_{i=1}^n \max \bar \theta_i(\lambda) + \sum_{F} \max \bar \theta_F(\lambda)

Note que para cada fator de entrada temos um subproblema de maximização, que chamamos de escravo. Afirmamos que L(\lambda) \ge MAP para qualquer valor de \lambda. Para ver porquê, considere uma atribuição ótima x^* que maximiza MAP. Temos que,

\bar \theta_i(\lambda) \ge \left( \theta_i (x^*_i) + \sum_{F:i\in F}  \lambda_{F,i}(x^*_i) \right)

pois x^* pode ser usado para obter \bar \theta_i(\lambda). Analogamente, temos que

\bar \theta_F(\lambda) \ge \left( \theta_F(\pmb{x}^*_F) - \sum_{i \in F} \lambda_{F,i}(x^*_i) \right)

Dado que MAP(\theta) = \sum_{i = 1}^n (\theta_i (x^*_i) + \sum_{F:i\in F}  \lambda_{F,i}(x^*_i))
+ \sum_{F} (\theta_F(\pmb{x}^*_F) - \sum_{i \in F}  \lambda_{F,i}(x^*_i))

fica fácil ver que nossa afirmação vale. Observe que L é uma relaxação do MAP original criada a partir de sua decomposição em duas partes e fornecendo um limitante dual, dando razão ao nome decomposição dual.

Note que os escravos não necessariamente precisam ser apenas sobre um fator apenas. A ideia é decompor o problema em subproblemas que possam ser resolvidos de maneira eficiente, como por exemplo os casos especiais descritos acima.

Algoritmo de decomposição dual

A ideia do algoritmo é começar com dois subproblemas independentes, mas usar as variáveis \lambda para amarrar os problemas. Conforme as iterações passam, o algoritmo vai penalizando aquelas cujos valores discordam nos diferentes subproblemas.

  • Inicializa todos os \lambda‘s com 0
  • Repita para t = 1, 2, \cdots
    • Obtenha a solução que otimiza cada subproblema (escravo):
      • \pmb{x}^*_F = \mbox{argmax}_{\pmb{x}_F}\{\bar \theta_i(\lambda)\} \forall F
      • x^*_i = \mbox{argmax}_{x_i}\{\bar \theta_F(\lambda)\} \forall i
    • Obteremos duas soluções ótimas possivelmente distintas
    • Para todo fator F e i \in F:
      • Se x^*_{F_i} \neq x^*_i, então faça
        • \lambda_{F,i}(x^*_{F_i}) =  \lambda_{F,i}(x^*_{F_i}) + \alpha_t
        • \lambda_{F,i}(x^*_{i}) =  \lambda_{F,i}(x^*_{i}) - \alpha_t

Convergência. Uma condição suficiente para a convergência é a escolha do aumento das penalidades \alpha:

\sum_t \alpha_t \rightarrow \infty e \sum_t \alpha^2_t < \infty

Mesmo obtida a convergência, não há garantias de que as soluções obtidas nos subproblemas concordam. Se for este o caso, temos uma atribuição ótima para o MAP. Caso contrário, devemos resolver o problema de decodificação, que consiste em obter uma atribuição viável para o MAP. Existem várias heurísticas como por exemplo usar o valor das variáveis ótimas de um só escravo.

Conclusões

Fiquei interessado no conceito de matrizes supermodulares e como essa característica permite reduzir o MAP para o problema do corte mínimo.

Achei essa heurística de quebrar em subproblemas para resolvê-los de maneira exata muito parecida com a relaxação lagrangeana. Inclusive, o algoritmo usado para resolvê-la é o de otimização de subgradiente, que usei em outro post para resolver um exemplo de relaxação lagrangeana.

Uma curiosidade que me surgiu: dado que essa relaxação retorna limitantes duais, será que não é interessante utilizar um algoritmo do tipo branch-and-bound para resolvê-lo?

Essa parte de inferência exata me agradou bastante, pois tem muita interseção com minha área de concentração que é otimização combinatória.


Métodos de Amostragem

Nesta seção vamos ver uma maneira aproximada de se obter uma distribuição de probabilidade através de amostragem. Vamos apresentar os tradeoffs entre número de amostragem e precisão da estimativa. Também descreveremos um algoritmo para obter amostras de uma rede Bayesiana.

Para o caso mais geral que funciona para redes de Markov, apresentaremos as cadeias de Markov, as cadeias de Gibbs e um algoritmo simples para obter amostras sobre essas cadeias.

Finalmente falaremos de um algoritmo mais complexo com convergência mais rápida, chamado de Metropolis-Hastings.

Estimação baseada em amostragem

Seja {\cal D} um conjunto de amostras de uma variável X, dado por \{X[1], X[2], \cdots, X[M]\} e obtido de uma distribuição P(X) de maneira independente e idêntica ou IID (Independent and Identically Distributed)

Uma estimativa para o valor de P(X = x) é dado por T_{\cal D} = \dfrac{1}{M} \sum_{m=1}^M 1(X[m] = x), onde M é o número rotal de amostras e 1(X[m] = x) é a função indicador, que retorna 1 se X[m] = x e 0 caso contrário. Em palavras, T_{\cal D} é a fração das amostras cujo valor é igual a x

Limitantes. Há dois resultados úteis que ajudam na estimativa do erro das amostras em função do número de amostras. A desigualdade de Hoeffding,

P_{\cal D}(T_{\cal D} \not \in [p - \epsilon, p + \epsilon]) \le 2e^{-2 M \epsilon^2}

Ou seja, dado um erro \epsilon e M amostras, a probabilidade de que a diferença entre T_{\cal D} e p seja no máximo \epsilon é no máximo 2e^{-2 M \epsilon^2}.

Temos também a desigualdade de Chernoff,

P_{\cal D}(T_{\cal D} \not \in [p(1 - \epsilon), p(1 +\epsilon)]) \le 2e^{-2 Mp \epsilon^2 / 3}

O número de amostras para garantir esses limitantes pode ser impraticável. No caso do limitante aditivo, podemos reescrever a desigualdade de Hoeffding como:

M \ge \dfrac{\ln(2/\delta)}{2\epsilon^2}

Como o número de amostras necessários para se obter um erro menor que \epsilon com probabilidade maior que 1 - \delta.

No caso do limitante multiplicativo, a desigualdade de Chernoff diz que

M \ge 3 \dfrac{\ln(2/\delta)}{p\epsilon^2}

Estes limitantes são ruim de se obter para probabilidades pequenas porque \epsilon tende a ser menor ainda e neste caso M deve ser bastante grande.

Amostragem progressiva. (forward sampling) Para obter amostras dados os CPD’s de uma rede Bayesiana, podemos processar as variáveis em ordem topológica. Uma vez amostrado os valores das variáveis ancestrais, fazemos a amostragem dos descendentes considerando a distribuição condicional para os valores observados.

Obtidas várias amostras, podemos estimar a distribuição de probabilidade. No caso em que queremos estimar uma distribuição dado um conjunto observado, podemos usar a rejeição de amostra quando as variáveis setadas na amostra não corresponderem ao valor observado. Note que quanto mais variáveis são observadas, maior o número de rejeições das amostras e maior o número de amostragens necessárias para estimar o valor da distribuição.

Esse método de amostragem não funciona para redes de Markov pois não há a noção de ordem topológica para grafos não direcionados.

Cadeia de Markov

Uma cadeia de Markov é um grafo direcionado onde os vértices correspondem a estados e as arestas formam o modelo de transição, pois a cada aresta (x, x') está associado um valor T(x \rightarrow x') que representa a probabilidade de ir do estado x para um novo estado x'.

Temos que a soma das probabilidades de se sair de um estado x deve ser 1, ou seja, \sum_{x'} T(x \rightarrow x') = 1.

A dinâmica temporal define a relação de recorrência para se determinar a probabilidade de se estar em um estado em uma iteração a partir da iteração anterior.

P^{(t+1)}(X^{(t+1)} = x') = \sum_{x} P^{(t)} (X^{(t)} = x) T(x \rightarrow x')

Se a probabilidade de se estar em um estado estabiliza entre duas iterações consecutivas, temos a distribuição estacionária, ou seja,

P^{(t)}(x') \approx P^{(t+1)}(x') = \sum_{x} P^{(t)} (x) T(x \rightarrow x')

Neste caso ao invés de P^{(t)} usamos a representação

(1) \pi(x') = \sum_{x} \pi(x) T(x \rightarrow x')

O valor de \pi pode ser encontrado analiticamente resolvendo o conjunto de equações dado por (1) para todo estado x, além da restrição de que a soma das probabilidades deve somar 1, ou seja, \sum_{x} \pi(x) = 1.

Não é sempre que uma cadeia de Markov converge para um estado estacionário. Uma condição suficiente para que isso ocorra é que a cadeia seja regular. Uma cadeia de Markov é dita regular se existe um inteiro k tal que no passo k, a probabilidade de sair de x para x' é maior que 0, para todo par de estados x,x'.

Teorema. Uma rede de Markov regular converge para um estado estacionário único independentemente do estado inicial.

Por sua vez, uma condição suficiente para que uma cadeia de Markov seja regular é que quaisquer dois estados x, x' estejam conectados (ou seja, que exista um caminho direcionado de x para x' e vice-versa) e para todo estado, exista uma auto-transição (auto-aresta).

Algoritmo de Monte Carlo para cadeias de Markov

Suponha que queiramos estimar uma distribuição de probabilidade P(X) através de amostragem, mas que não exista uma forma direta de obter amostras. Se pudermos definir uma cadeia de tal forma que a distribuição \pi corresponda a P(X), podemos usar o algoritmo de Monte Carlo para Cadeias de Markov (MCMC).

Primeiramente, obtenha uma amostra da distribuição inicial P^{(0)}. Depois, execute t iterações, em cada uma obtendo uma nova amostra a partir da anterior através de um passeio aleatório (random walk), usando o modelo de transição.

O número de repetições t deve ser tal que a distribuição P^{(t)} se aproxime de \pi, caso em que dizemos que a distribuição está misturada. Note porém, que não temos como garantir que em uma dada iteração t, a distribuição P^{(t)} esteja misturada, mas podemos fazer medições estatísticas que podem nos dizer quando ela não está.

Uma abordagem é executar o algoritmo começando com duas amostras distintas e computar a distribuição para cada uma das execuções usando as k últimas amostras obtidas. Para cada estado da cadeia, plotamos em um gráfico de dispersão um ponto com as coordenadas (p_1, p_2), onde p_1 e p_2 são as probabilidades para o estado em cada uma das execuções.

Se a distribuição estiver misturada, o gráfico de dispersão tende a se concentrar na diagonal. Na Figura abaixo à esquerda, podemos concluir que não houve mistura, enquanto no gráfico da direita pode ter havido mistura.

Figura 1: Gráfico de dispersão

Uma vez que acreditemos que a distribuição tenha misturado, podemos obter mais um conjunto de amostras executando o algoritmo a partir de P^{(t)} para obter a distribuição \pi.

Amostragem de Gibbs

Vimos que podemos obter amostras a partir de uma cadeia de Markov para estimar a distribuição \pi. Mas como construir uma cadeia de forma que \pi corresponda a P?

No caso de uma distribuição de Gibbs P_{\Phi}(X_1, \cdots, X_n), que corresponde a um produto de fatores, podemos construir uma cadeia de Markov para estimar P_\Phi.

Para tanto, considere uma cadeia de Markov onde os estados são todas as possíveis atribuições do conjunto \{X_1, \cdots, X_n\}. O modelo de transição T(\pmb{x} \rightarrow \pmb{x}') de um estado \pmb{x} é dado com o seguinte algoritmo:

  • Para cada i=1,\cdots, n, obtenha o valor de x_i através da distribuição P_\Phi(X_i | \pmb{x}_{-i}), onde \pmb{x}_{-i} corresponde à atribuição \pmb{x} exceto a componente x_i. (Note que o novo valor de x_i é usado em P_\Phi(X_j | \pmb{x}_{-j}) para j > i).
  • Depois de obtido o valor de todas as componentes, fazemos \pmb{x}' = \pmb{x}.

Essa cadeia de Markov é chamada de cadeia de Gibbs.

É possível mostrar que a distribuição estacionária desta cadeia é P_\Phi e portanto podemos usar o algoritmo da seção anterior sobre essa cadeia para estimá-la. Antes porém, precisamos definir como obter a distribuição P_{\Phi}(X_i | \pmb{x}_{-i}).

Temos por defnição:

P_{\Phi}(X_i | \pmb{x}_{-i}) = \dfrac{P_{\Phi}(X_i, \pmb{x}_{-i})}{P_{\Phi}(\pmb{x}_{-i})}

É possível ver que as contantes de normalização são as mesmas para \tilde P_{\Phi}(X_i, \pmb{x}_{-i}) e \tilde P_{\Phi}(\pmb{x}_{-i}) e portanto podemos reescrever

P_{\Phi}(X_i | \pmb{x}_{-i}) = \dfrac{\tilde P_{\Phi}(X_i, \pmb{x}_{-i})}{\tilde P_{\Phi}(\pmb{x}_{-i})}

Temos que \tilde P_{\Phi}(X_i, \pmb{x}_{-i}) é o produto de todos os fatores dado um conjunto de observação \pmb{x}_{-i}, ou seja, \tilde P_{\Phi}(X_i, \pmb{x}_{-i}) = \prod_{\phi_j \in \Phi} \phi_j.

Também podemos obter \tilde P_{\Phi}(\pmb{x}_{-i}) a partir de \tilde P_{\Phi}(X_i, \pmb{x}_{-i}) marginalizando a variável X_i, ou seja,

\tilde P_{\Phi}(\pmb{x}_{-i}) = \sum_{X_i} \tilde P_{\Phi}(X_i, \pmb{x}_{-i})

Veja que é possível colocar pra fora do somatório os fatores que não contêm X_i. Assim, ao fazer a razão entre
\tilde P_{\Phi}(\pmb{x}_{-i}) e \tilde P_{\Phi}(X_i, \pmb{x}_{-i}), vemos que esses fatores poderão ser cancelados. Além disso, observamos que \tilde P_{\Phi}(\pmb{x}_{-i}) é uma constante. Com isso chegamos à seguinte relação

P_{\Phi}(X_i | \pmb{x}_{-i}) \propto \prod_{j: X_i \in \mbox{Escopo}[\phi_j]} \phi_j(X_i, \pmb{x}_{j, -i})

Onde \pmb{x}_{j, -i} corresponde à atribuição das variáveis do escopo de \phi_j exceto X_i.

Propriedades. Uma cadeia de Gibbs não necessariamente é uma cadeia de Markov regular. Um contra-exemplo é a rede bayesiana do ou-exclusivo determinístico, Y = X_1 \oplus X_2 com conjunto de observação Y = 1. Podemos verificar que dado um estado inicial X_1=0, X_2=1, nunca alcançaremos o estado X_1=1, X_2=0.

Uma condição suficiente para que uma cadeia de Gibbs seja regular é que os fatores sejam positivos, ou seja, não tenham entradas negativas ou nulas, condição claramente não respeitada por CPD’s determinísticos.

Uma desvantagem deste método é que geralmente o número de iterações necessários para que a distribuição misture é alto. Além disso, devemos ser capaz de obter amostras dos fatores, o que é complicado para fatores que não são representados na forma tabular, como por exemplo aqueles envolvendo variáveis contínuas ou cujo fatores sejam muito grandes, o que pode tornar a tarefa tão difícil quanto computar P_\Phi explicitamente.

Algoritmo Metropolis-Hastings

Uma cadeia de Markov é dita reversível se para toda aresta da forma T(x \rightarrow x') há uma aresta na direção contrária T(x \rightarrow x') respeitando a condição de balanço detalhado, que é dada por

\pi(x)T(x \rightarrow x') = \pi(x')T(x' \rightarrow x)

A condição de balanço detalhado é uma condição suficiente para que haja uma distriuição estacionária única \pi.

Para descrever o algoritmo de Metropolis-Hastings, definiremos dois conceitos: a distribuição proposta Q(x \leftarrow x') e a probabilidade de aceitação A(x \leftarrow x') as quais explicaremos a seguir.

Em cada estado x, ao invés de usar T, obtemos a amostra x' a partir da distribuição proposta Q(x \rightarrow x'). Então, com uma probabilidade A(x \leftarrow x') aceitamos a proposta. Do contrário permanecemos no estado x.

Podemos escrever T(x \leftarrow x') em função de Q e A como:

Se x' \ne x,

T(x \rightarrow x') = Q(x \rightarrow x') A(x \rightarrow x')

Se x' = x,

T(x \rightarrow x') = Q(x \rightarrow x) \sum_{x'' \neq x} Q(x \rightarrow x'') (1 - A(x \rightarrow x''))

Se quisermos que T respeite a condição de balanço detalhado, devemos ter:

\pi(x)T(x \rightarrow x') = \pi(x')T(x' \rightarrow x)

Substituindo por A e Q,

\pi(x)Q(x \rightarrow x') A(x \rightarrow x') = \pi(x')Q(x \rightarrow x') A(x \rightarrow x')

Supondo \pi(x')Q(x \rightarrow x') \le \pi(x)Q(x \rightarrow x'), ficamos com a seguinte razão

\dfrac{A(x \rightarrow x')}{A(x \rightarrow x')} = \dfrac{\pi(x')Q(x \rightarrow x')}{\pi(x)Q(x \rightarrow x')} = \rho \le 1

Se escolhermos as probabilidades A(x \rightarrow x') = \rho e A(x' \rightarrow x) = 1, satisfazemos a propriedade de balanço detalhado. Mais sucintamente podemos definir A como

A(x \rightarrow x') = \min\left\{1, \dfrac{\pi(x')Q(x \rightarrow x')}{\pi(x)Q(x \rightarrow x')}\right\}

A escolha de Q é uma ciência inexata e segundo a Wikipedia, uma distribuição comumente utilizada é uma Gaussiana com média x.

Conclusão

Achei bem difícil acompanhar essas últimas aulas referentes a amostragem. O algoritmo de Metropolis-Hastings em especial ficou muito vago, mesmo tentando complementar com a Wikipedia. Não consegui entender direito qual a motivação por trás dos conceitos de distribuição proposta e da probabilidade de aceitação.

Intuitivamente, me parece que na cadeia de Gibbs sempre tomamos passos “gulosos” (não é estritamente guloso porque é probabilístico, mas ele tende a ir por arestas de maior probabilidade). No algoritmo de Metropolis-Hastings, parece que podemos ter uma distribuição que permite ir para estados menos prováveis, evitando dar passos muito pequenos e a probabilidade de aceitação é usada para garantir a propriedade de balanço detalhado.


Probabilistic Graphical Models – Semana 4

Abril 15, 2012

Este post apresenta as notas de aula da semana 4 do curso online Probabilistic Graphical Models.

Propriedades da propagação de crença

Para concluir o algoritmo de passagem de mensagens em grafos cluster, vamos apresentar mais algumas definições.

Calibragem. Um cluster grafo é calibrado se para toda aresta (C_i, C_j), ambos têm a mesma crença para as variáveis no sepset entre eles, ou seja, a crença marginalizada de C_i e C_j para S_{ij} é igual.

Matematicamente temos que

\sum_{C_i \setminus S_{ij}} \beta(C_i) = \sum_{C_j \setminus S_{ij}} \beta(C_j)

Podemos mostrar que se houve convergência nas mensagens o grafo cluster está calibrado.

Definimos a crença do sepset como o produto das mensagens trocadas na aresta correspondente, ou seja:

\mu_{ij}(S_{ij}) = \delta_{j \rightarrow i} (S_{ij}) \delta_{i \rightarrow j} (S_{ij})

Reparametrização. é uma propriedade que grafos cluster têm devido à seguinte relação:

\dfrac{\prod_i \beta_i}{\prod_{ij} \mu_{ij}} = \tilde P_{\Phi} (X_1, \cdots, X_n)

O nome da propriedade vem do fato de que esta é uma forma alternativa de parametrizar uma distribuição de probabilidade usando crenças. Essa relação pode ser deduzida da seguinte maneira:

\dfrac{\prod_i \beta_i}{\prod_{ij} \mu_{ij}} = \dfrac{\prod_i (\psi_i \prod_{j \in N_i} \delta_{j \rightarrow i})}{\prod_{ij} \delta_{i \rightarrow j}} =

A observação chave é notar que cada \delta_{i \rightarrow j} tanto do denominador na fração acima, aparece exatamente uma vez para cada aresta. Cancelando os termos, o que sobre é

= \prod_{i} \psi_i = \tilde P_{\Phi} (X_1, \cdots, X_n)

que é definição da distribuição não normalizada.


Passagem de mensagens em Árvore-Clique

No caso particular de um caminho, é fácil mostrar que a crença em um vértice corresponde exatamente à distribuição marginalinada contendo as variáveis desse vértice.

Uma árvore-clique é uma árvore não-direcionada, onde os vértices correspondem a clusteres, a cada aresta está associada um sepset, S_{ij}, que neste caso é exatamente igual a C_i \cap C_j.

Sendo um caso especial de um grafo cluster, temos que satisfazer as propriedades mencionadas no último post, mais especificamente, a preservação de familiar e o caminho de interseções.

Note que a unicidade de caminhos é trivialmente garantido numa árvore, mas para que a propriedade do caminho de interseções ocorra, aindas precisamos garantir que vértices contendo uma mesma variável X estão conectados por um caminho contendo X.

Estas propriedades são suficientes para garantir que a crença obtida ao final do algoritmo corresponda às distribuições marginalizadas para as variáveis no cluster. Para esse tipo especial de grafo, o algoritmo de passagem de mensagem é uma variante do algoritmo de eliminação de variáveis.

Algoritmo de passagem de menssagem

Dizemos que uma mensagem é final se ela não se modificará independente de quantas iterações do algoritmo de passagem de mensagens sejam executadas. Com essa definição, temos a seguinte propriedade:

Propriedade 1. quando um vértice C_i recebe mensagens finais de todos os seus vizinhos exceto C_j, então \delta_{i \rightarrow j} é final.

Vamos apresentar um algoritmo de passagem de mensagem para árvores-cliques que se beneficia dessa propriedade, para obter mensagens finais usando apenas uma iteração. Primeiramente, note que uma folha da árvore satisfaz a Propriedade 1 acima trivialmente para seu único vizinho C_j.

A primeira parte do algoritmo consiste em sempre processar folhas dessa árvore, enviando a mensagem para seu vizinho e removendo-a em seguida. É fácil ver que as folhas do grafo resultante vão enviar mensagens finais para seu único vizinho, já que por hipótese, todos os seus outros vizinhos já lhe enviaram suas mensagens.

Desta maneira, eventualmente processaremos todos os vértices, obtendo as mensagens finais de cada aresta, porém em apenas uma das direções. Como fazer para obter as mensagens na direção oposta?

Considere o último vértice v' a ser processado no algoritmo acima. Temos que todos os seus vizinhos no grafo original já lhe enviaram mensagens na primeira etapa. Assim, todas as mensagens enviadas por v' a partir de agora são finais, devido à Propriedade 1.

Seja v'' o penúltimo vértice removido. Se ele enviou mensagem para v', então agora v' já lhe enviou uma mensagem de volta. Independentemente de ter enviado a mensagem para v' ou não, agora qualquer mensagem enviada por v'' é final.

Esse argumento construtivo nos sugere a segunda parte do algoritmo: processar os vértices na ordem reversa daquela da primeira etapa e cada vértice processado deve enviar uma mensagem para todos os que lhe enviaram mensagem na primeira fase.

Note portanto, que é possível encontrar uma ordem de passagem de mensagens em uma árvore de modo que com apenas uma iteração e todas as mensagens sejam finais. Note que nesta única iteração, foram enviadas exatamente 2(n-1) mensagens, onde n é o número de vértices/clusters.

Consultas

Uma vez calculadas as crenças nos clusters, podemos fazer consultas do tipo P(Y|E=e). Supondo que faremos essa consulta várias vezes para diferentes valores de e, podemos usar o algoritmo acima como um pré-processamento.

Uma vez calculadas as crenças, se existe uma clique contendo Y e E, podemos simplesmente marginalizar a crença para ficar com \tilde P(Y,E) e então remover as entradas E \ne e (reduzir).

Por outro lado, se não existe uma clique contendo Y e E, teremos que considerar duas cliques, uma contendo Y, denotada por C_Y e outra contendo E, denotada por C_E. Primeiramente, reduzimos a crença em C_E, o que ocasionará uma mudança nas mensagens enviadas por esse cluster e depois fazer a consulta no cluster C_Y. Entretanto, nesse caso é possível reaproveitar até n-1 sem recomputá-las.

Para isso, na primeira fase do algoritmo acima, o último vértice a ser processado deve ser C_E (note que é sempre possível escolher uma ordem de processamento em que isso acontece). Nenhuma das n-1 mensagens enviadas nessa primeira etapa partiu de C_E e portanto não dependem de \psi(C_E). Logo, só a segunda etapa do algoritmo precisa ser recomputada para cada consulta de E = e.

Porém, essa ideia é atrelada ao conjunto E. Se a consulta for feita usando outro conjunto de observação, a primeira fase deverá ser computada novamente.

Árvores-clique e independência

Seja T uma árvore-clique e (C_i, C_j) uma aresta de T. Como T satisfaz as propriedades de preservação familiar e o caminho de interseções, podemos particionar o conjunto de variáveis em três conjuntos.

Considere as duas componentes resultantes da remoção da aresta (C_i, C_j), denotando por W_i aquela que contém C_i e W_j aquela que contém C_j.

Seja W_{<(i,j)} o conjunto de variáveis que pertencem a algum cluster em W_i, mas não pertencem a S_{ij}, ou seja, W_{<(i,j)} = (\cup_{C \in W_i} C) \setminus S_{ij}

Analogamente, defnimos W_{<(j,i)} como W_{<(j,i)} = (\cup_{C \in W_j} C) \setminus S_{ij}

Observamos que W_{<(i,j)} e W_{<(j,i)} são disjuntos. Caso contrário, existiria uma variável X em W_{<(i,j)} e em W_{<(j,i)}, mas não em S_{i,j}. Isso implica que não existe um caminho contendo X conectando esses dois clusters em que X aparece, contradizendo a hipótese de que a árvore satisfaz a propriedade do caminho de interseções.

Com essa observação, podemos ver que esses três conjuntos particionam, de fato, o conjunto de variáveis em T. O seguinte teorema associa a ideia de independência ao sepset:

Teorema: Uma árvore T satisfaz a propriedade de caminhos de interseção ou RIP (Running intersection property) se e somente se, para toda aresta (i,j), P_\Phi \models (W_{<(i,j)} \perp W_{<(j,i)} | S_{ij})

Podemos ver que em uma árvore-clique correta, se observarmos as variáveis no sepset, as variáveis nas componentes W_{<(i,j)} e W_{<(j,i)} tornam-se independentes. O nome sepset agora fica sugestivo como conjunto separador.

Árvore-clique e eliminação de variáveis

Podemos construir uma árvore-clique a partir dos fatores intermediários obtidos no algoritmo de eliminação de variáveis.

Recapitulando, a cada passo i escolhemos uma varíavel Z para ser eliminada, obtendo um fator

\Phi'_i = {\phi_k \in \Phi_i | Z \in \mbox{Escopo}(\phi_k)}

\psi_i = \prod_{\phi_k \in \Phi_i'} \phi_k

\lambda_i = \mbox{Escopo}(\psi_i)

\tau_i = \sum_{Z} \psi_i

\Phi_{i+1} \leftarrow  \Phi_i \setminus \Phi' \cup \{\tau_i\}

Nesse caso, construímos uma árvore-clique onde os clusters correspondem às variáveis em \lambda_i os quais denotamos por C_{\lambda_i}. Existe uma aresta entre um cluster C_{\lambda_i} e C_{\lambda_j} se \tau_i foi usado no produtório para calcular \psi_j. Nesse caso, o sepset é igual ao domínio de \tau_i.

Após a construção dessa árvore-clique, pode haver clusters que são subconjuntos de outros cluster, o que é redundante. Assim, uma fase de pós-processamento consiste em removê-los.

Propriedades

O grafo cluster resultante é de fato uma árvore porque um \psi_i é usado em exatamente uma vez em alguma iteração posterior j > i. Sendo n o número de iterações, podemos ver que o grafo tem n-1 arestas e n nós e também não possui ciclos, logo concluímos que é uma árvore.

É possível ver que o grafo resultante respeita a preservação de família pois todo fator usado no algoritmo possui um escopo que é subconjunto de algum \lambda_i.

Temos também que o grafo resultante respeita a propriedade de caminho de independência. Considere uma variável X que aparece em dois clusters diferentes, C_i e C_j. Segundo o algoritmo, existe um caminho contendo X de cada um desses clusters até o cluster onde essa variável é eliminada. Assim, como C_i e C_j têm caminhos com vértice em comum na árvore, podemos concluir que existe um caminho único entre estes dois clusters.

Complexidade

Podemos reduzir a eliminação de variáveis para o algoritmo de passagem de mensagens em uma árvore-clique. Para calcular a marginalização de um conjunto de variáveis usamos aproximadamente o mesmo número de operações em ambos os casos. A vantagem da passagem de mensagens é o re-uso de mensagens no caso de querermos computar a maginalização de outro conjunto de variáveis.

Propagação de Crença na prática

Em um dos vídeos são discutidos aspectos práticos da propagação de crença, como convergência e corretude.

Variantes do BP

Há dois principais meios de passagem de mensagem: síncrono e assíncrono. No caso síncrono as mensagens são enviadas ao mesmo tempo enquanto no assíncrono as mensagens são enviadas em ordem.

A vantagem do caso síncrono é que o algoritmo é facilmente paralelizável. Porém, na prática o assíncrono tende a convergir mais rápido. O problema da abordagem assíncrona é que a convergência e corretude dependem da ordem utilizada para passar as mensagens.

Heurísticas

Há algumas heurísticas para a escolha da passagem das mensagens. Uma delas é o TRP (Tree reparametrization) que consiste em cada iteração escolher árvores geradoras e usar o algoritmo de árvore-clique. Outra abordagem é passar mensagens por arestas entre clusters com a crença marginalizada para o sepset mais distantes.

Outras alternativas consistem em modificar o modo como a mensagem é calculada na tentativa de obter convergência. Uma delas é as mensagens suaves (smoothing messages), onde uma mensagem é construída com uma combinação convexa (suavização) entre seu valor anterior e o novo valor, ou seja,

\delta_{i \rightarrow j} \leftarrow \lambda (\sum_{C_i \setminus S_{ij}} \psi_i \prod_{k \neq j} \delta_{k \rightarrow j}) + (1 - \lambda) \delta_{i \rightarrow j}^{old}

Conclusão

A ordem de passagem de mensagens usada pelo algoritmo em árvores-cluster não foi muito detalhado e aqui apresentei uma versão própria que me pareceu correta. Entretanto para consultas, este algoritmo é atrelado ao conjunto observado E, sendo necessário computar todas as mensagens novamente quando muda E. Porém, acho que é sempre possível aproveitar as mensagens computadas na primeira fase, independente de E.


Inferência de MAP

Vamos agora focar no tipo de consulta de probabilidade a posteriori máxima ou MAP, ou seja,

\mbox{MAP}(Y|E=e) = \mbox{argmax}_y\{P(Y = y | E = e)\}

Algoritmos Max-Soma

Dada uma distribuição P_\Phi(X) \propto \prod_k \phi_k (D_k), queremos calcular a entrada x para a qual P_\Phi(X=x) é máxima. Vimos que isso é equivalente a

\mbox{argmax} \prod_k \phi(D_k)

Agora, vamos transformar o produtório em somatório. Seja \theta_k = \log (\phi_k). Podemos reescrever a equação acima como:

\mbox{argmax} \sum_k \theta(D_k)

A razão para usar log é por questões práticas. Como não estamos interessados no valor da distribuição máxima e sim na atribuição que corresponde a esse valor, não precisamos ficar lidando com multiplicações, que é computacionalmente mais custoso do que somas.

Podemos definir operações análogas ao processo de obtenção de probabilidades marginais. O produto de fatores agora corresponde à soma de log de fatores. Vamos definir um análogo à distribuição de probabilidade conjunta:

\theta(X) = \sum_k \theta(D_k).

Da mesma forma, a marginalização, que é a eliminação de um conjunto de variáveis, pode ser definida para o MAP. Lembrando, dado um fator \phi com domínio Z e um conjunto de variáveis X que queremos remover via marginalização, o que fazemos é somar entradas de \phi com valores de Z \setminus X iguais. Agora, definimos a marginalização-max como o máximo entre as entradas de phi com valores de Z \setminus X iguais.

A ideia dos algoritmos que iremos apresentar é obter um fator com um conjunto menor de variáveis Y sobre o qual se deseja fazer a consulta MAP(Y|E=e). Como estamos trabalhando com somas seguidas da função max, temos os algoritmos max-soma.

Os algoritmos de eliminação de variável e passagem de mensagens pode ser adaptados para esse caso, trocando o produto de fatores pela soma do log de fatores e na marginalização a função soma é trocada pela função max. E no caso da passagem de mensagem, as crenças agora correspondem a distribuição de probabilidade max-marginal.

Obtendo a atribuição MAP

Para este caso existe uma propriedade interessante quando o grafo cluster está calibrado. Considere os clusters adjacentes latex C_i = \{A, B\} e C_j = \{B, C\}. Se o sepset entre C_i e C_j é \{B\}, então a crença marginalizada para B de C_i e C_j são iguais.

Seja b^* o valor para o qual a crença marginalizada para B é máxima e seja tal valor P^*. Isso quer dizer que existe uma entrada na crença \beta_i com valor (a^*, b^*) tal que \beta_i(a^*, b^*) = P^* e em \beta_j existe (c^*, b^*) tal que \beta_j(c^*, b^*) = P^*.

Dessa forma, podemos obter a atribuição de A, B, C que geram MAP, simplesmente encontrando a atribuição em \beta_i e \beta_j separadamente.

Porém, isso só funciona se existir uma única atribuição que gere o valor máximo. Do contrário, podemos usar algumas alternativas como perturbação dos fatores para garantir uma atribuição ótima única; ou ao construir a atribuição ótima incrementalmente, armazenar os valores já atribuídos e ao processar uma nova clique verificar a consistência dos valores.


Probabilistic Graphical Models – Semana 3

Abril 8, 2012

Este post apresenta as notas das aulas da semana 3.

Engenharia do conhecimento

Essa aula discute aspectos práticos da modelagem por redes Bayesianas ou de Markov. Não existe necessariamente um modelo correto, o que caracteriza a tarefa de criação desses modelos como uma arte e não como uma ciência.

Os modelos podem ser classificados em diversos tipos, dependendo de suas características. A seguir são descritas situações em que um tipo é mais indicado que outro, servindo mais como guia, não como regra, notando também que é possível usar modelos híbridos.

Template vs. específico. Se o problema possui variáveis com características compartilhadas, é possível usar um modelo por templates, que tem um número menor de variáveis em comparação com um modelo específico, que possui um número grande de variáveis distintas.

Generativo vs. discriminativo. Um modelo discriminativo é mais apropriado quando se tem uma tarefa particular. Nesse caso definimos características bem precisas, de modo a ter baixa correlação entre elas. Isso em geral garante um bom desempenho do modelo. Por outro lado, um modelo generativo é usado quando a tarefa não está bem especificada, pois ele é mais flexível e mais fácil de treinar.

Tipos de variáveis

Podemos classificar as variáveis de um modelo em 3 tipos:

  • alvo – as que queremos prever (saída)
  • observadas (entrada)
  • latente – escondidas. Não são explicitamente usadas, mas podem ser usadas para simplificar o modelo. Como no exemplo da Figura 1, onde a variável GMT serve apenas para não termos que definir a dependência entre W_1, \cdots W_k explicitamente usando arestas.

Figura 1: Varíavel GMT é latente

Estrutura

Em especial para redes Bayesianas, decidir a causalidade, ou seja, a orientação das arestas, não é uma tarefa simples. Dependendo da escolha, podemos ter um modelo maior ou menor. Por exemplo, se no exemplo da Figura 1 quiséssemos inverter a orientação das arestas, tornaríamos as variáveis W indpendentes, o que exigiria adicionar novas arestas entre elas.

Refinamento iterativo do modelo

Uma abordagem na criação dos modelos é o uso de iterações para refinar o modelo e ajustar parâmetros.

Na minha opinião, a modelagem de redes Bayesianas principalmente, lembra muito a modelagem de classes em engenharia de sofware, onde precisamos decidir o relacionamento entre classes (arestas) e sua hierarquia (orientação). Processos iterativos são bastante utilizados nesses casos também.


Consultas de probabilidade condicional

Uma operação comum em redes Bayesianas e de Markov é determinar a distribuição de probabilidade dado uma obervação. Suponha que queiramos determinar a distribuição de um conjunto de variáveis Y dado uma evidência E = e, ou seja, queremos calcular P(Y|E = e).

Uma alternativa é calcular a distribuição conjunta e depois fazer operações de marginalização e redução. Porém, vimos que o tamanho de um fator é exponencial no número de variáveis do modelo. Será que não há um meio mais eficiente, visto que estamos interessados em apenas uma consulta específica?

Provavelmente não, pois esse problema é NP-difícil. Mesmo o problema de, dada uma rede \Phi, decidir se P_\Phi(X = x) > 0 para um dado x é NP-difícil.

Muitos algoritmos exatos ou heurísticos usados para calcular essa distribuição são chamados de soma-produto. Isso porque a distribuição conjunta é equivalente à mulitiplicação dos fatores e a marginalização corresponde à soma de entradas do fator resultante, ou seja, temos a soma de produtos.

Vamos formalizar essa ideia. Usando a definição de probabilidade condicional, temos

P(Y | E = e) = \dfrac{P(Y, E = e)}{P(E = e)}

A distribuição P(Y, E = e) pode ser obtida a partir da distribuição conjunta removendo-se as variáveis do modelo W que não estejam nem em Y nem em E, ou seja, W = \{X_1, \cdots, X_n\} - Y - E, através de marginalização.

P(Y, E = e) = \sum_W P(Y, W, E=e)
\quad = \sum_W (\prod_{\phi_k \in \Phi} \phi_k(D_k, E = e))

Note que P(E = e) é equivalente a uma constante de marginalização. Como sabemos que P(Y | E = e) é uma distribuição, suas entradas devem somar 1, não precisamos calcular P(E = e) bastando re-normalizar P(Y, E = e).

Probabilidade a posteriori máxima

Outra consulta que podemos querer fazer em uma rede é a probabilidade a posteriori máxima ou MAP (Maximum a Posteriori). Essa consulta basicamente consiste em dizer qual o valor de maior probabilidade dado um conjunto de evidência e é importante para decidir o valor de uma variável.

Matematicamente podemos definir um MAP como:

\mbox{MAP}(Y|E=e) = \mbox{argmax}_y\{P(Y = y | E = e)\}

Nesse caso, temos que Y = \{X_1, \cdots, X_n \} - E.

Infelizmente, dada uma rede com distribuição \Phi, determinar uma atribuição x tal que P_\Phi(X = x) é máximo, também é NP-difícil.

Para esse problema, os algoritmos exatos e heurísticos são baseados em max-produto, pois consistem em encontrar o valor máximo de um produto de fatores. Note que nesse caso não temos soma (marginalização), pois não há nenhuma variável fora de Y e E a ser removida. Já vimos que:

P(Y | E = e) = \dfrac{P(Y, E = e)}{P(E = e)} \propto P(Y | E = e)

Também temos que P(Y, E = e) corresponde ao produto de fatores reduzidos, ou seja, sem as entradas correspondentes E \neq e. Vamos denotar por \phi'_k(D'_k) o fator reduzido \phi_k(D_k)

P(Y, E = e) = \frac{1}{Z} \prod_{k} \phi'_k(D'_k) \propto \prod_{k} \phi'_k(D'_k)

Para encontrar a variável que maximiza uma função, não importa a constante multiplicativa, de modo que

\mbox{argmax}_Y \{P(Y|E=e)\} = \mbox{argmax}_Y \{\prod \phi'_k (D'_k)\}


Eliminação de variáveis

Definidos esses dois problemas NP-difíceis, vamos apresentar algoritmos para resolvê-los. Primeiramente, vamos nos concentar na consulta de probabilidade condicional.

O algoritmo de eliminação de variáveis consiste em determinar a marginalização do produto de fatores, eliminando uma variável de cada vez, ao mesmo tempo em que calcula o produtório. Note que as variáveis podem ser eliminadas em qualquer ordem. Entretanto, veremos que a ordem em que as mesmas são eliminadas pode afetar a complexidade do problema.

Um passo do algoritmo é sucintamente descrito a seguir. Seja Z a variável que queremos eliminar. Denotamos por \Phi' todos os fatores quem têm Z em seu domínio, ou seja,

\Phi' = \{\phi_i \in \Phi : Z \in \mbox{Escopo}[\phi_i]\}

Multiplicamos esses fatores e obtemos um fator \psi:

\psi = \prod_{\phi_i \in \Phi'} \phi_i

Podemos fazer então a marginalização para remover Z, gerando um novo fator \tau:

\tau = \sum_{Z} \psi

Perceba que reduzimos a um novo problema com um conjunto de fatores sem a variável Z. Procedendo dessa maneira, eventualmente chegaremos no fator final.

\Phi := (\Phi \setminus \Phi') \cup \{\tau\}

Observamos que esse algoritmo funciona tanto para BN’s quanto para MN’s.

Complexidade

Podemos provar que a complexidade do algoritmo é O(Nm), onde N é o número de entradas do maior fator, m é o número de fatores. Entretanto, N \in O(d^r), onde d é a cardinalidade da maior variável e r é o tamanho do maior domínio considerando todos os fatores.

O tamanho do maior domínio depende da ordem em que as variáveis são eliminadas.

Encontrando uma ordem de eliminação

Vamos apresentar um algoritmo que visa encontrar uma ordem para melhorar a complexidade do algoritmo de eliminação de variáveis. Para isso, é conveniente definirmos grafos auxiliares.

Rede de Markov induzida. Da mesma forma que definimos uma rede de Markov induzida para a distribuição de Gibbs, podemos defini-la para uma rede Bayesiana se olharmos apenas os fatores. Lembrando, esse grafo induzido tem conjunto de vértices correspondendo às variáveis e uma aresta entre duas variáveis que aparecem no domínio de um mesmo fator.

Um detalhe é que duas variáveis que não estão conectadas no grafo induzido original, podem aparecer no domínio dos fatores intermediários gerados ao longo do algoritmo. Nesse caso, adicionamos novas arestas ao grafo, que serão chamadas arestas de preenchimento (fill edges).

O grafo induzido final é denotado por I_{\Phi,\alpha}, dado um conjunto de fatores \Phi e uma ordem \alpha. É fácil ver que todo fator usado no algoritmo de eliminação de variável corresponde a uma clique em I_{\Phi,\alpha}.

Contrariamente, temos outra observação que não é tão trivial de ver:

Teorema. toda clique maximal em I_{\Phi,\alpha} é um fator produzido durante o algoritmo.

Vamos definir mais alguns termos. A largura do grafo induzido corresponde ao tamanho da maior clique no grafo menos 1. A largura induzida minimal é a menor largura de I_{\Phi,\alpha} considerando todas as ordenações \alpha. Com isso, podemos enunciar o seguinte teorema:

Teorema. Encontrar a ordem de eliminação que gere a largura induzida minimal de I_{\Phi,\alpha} é NP-difícil.

Note porém que mesmo com uma ordem ótima, o problema de decidir a inferência ainda pode ser exponencial.

Heuristicas

Há diversas heurísticas para encontar uma boa ordem de eliminação. Apresentamos a seguir algumas gulosas baseando-se em critérios como:

  • escolha do vértice que produz o fator de menor domínio;
  • escolha do vértice que produz o fator com menor número de entradas;
  • escolha do vértice que produz o menor número de arestas de preenchimento.

Há também uma outra heurística baseada no seguinte teorema:

Teorema. O grafo induzido I_{\Phi,\alpha} é triangulado, i.e., não possui ciclos de tamanho maior que 3.

A heurística consiste então em encontrar uma triangulação com pouca largura do grafo induzido original e então é possível obter uma ordem que gera esse grafo.

Conclusões

Fiquei curioso para saber qual é a redução utilizada para mostrar a NP-completude dos problemas. Se eu tiver um tempo vou pesquisar.

A última observação ficou meio nebulosa. Como obter uma ordenação de um grafo triangulado?


Algoritmo de Passagem de mensagem

Nesta seção é apresentado um algoritmo baseado em passagem de mensagem entre vértices de um grafo. Trata-se de um algoritmo heurístico para obter distribuições marginalizadas de um dado subconjunto de variáveis.

Vamos a algumas definições.

Grafo cluster. Grafo não direcionado, onde o conjunto de vértices são clusters (subconjuntos de variáveis). Associadas às arestas, temos um subconjunto de variáveis chamado sepset, denotado por S_{i,j} \subseteq C_i \cap C_j

Dado um conjunto de fatores \Phi, atribuímos cada \phi_k \in \Phi a um cluster C_{\alpha(k)}, tal que \mbox{escopo}[\phi_k] \subseteq C_{\alpha(k)}, onde \alpha(k) corresponde ao índice do cluster ao qual o fator k foi atribuído.

Definimos o potencial de um cluster C_i, denotado por \psi(C_i) como o produto dos fatores atribuídos a ele, ou seja, \psi(C_i) = \prod_{k: \alpha(k) = i} \phi_k

Definimos a passagem de mensagens entre dois clusters C_i e C_j como um fator dado sucintamente por

\delta_{i \rightarrow j} (S_{ij}) = \sum_{C_i \setminus S_{ij}} \psi(C_i) \prod_{k \in (N(i) \setminus \{C_j\})} \delta_{k \rightarrow i}

lembrando que S_{ij} é o sepset e corresponde ao conjunto de variáveis que se quer transmitir. N(i) é o conjunto de vizinhos de C_i.

Em palavras, a definição de uma mensagem entre C_i e C_j é o produto do potencial de C_i (\psi(C_i)) pelas mensagens recebidas dos vizinhos de C_i exceto C_j, marginalizado para ficar com domínio igual a S_{ij}.

A crença (belief) de um cluster C_i, denotada por \beta(C_i) é definida como:

\beta(C_i) = \phi(C_i) \prod_{k \in N(i)} \delta_{k \rightarrow i}

ou seja, o produto de todas as mensagem recebidas dos vizinhos de C_i e seu potencial.

Com isso, podemos definir o algoritmo de propagação de mensagens. Basicamente, inicializamos cada vértice do grafo com o fator nulo e efetuamos trocas de mensagens por um certo número de informações.

Algoritmo de propagação de crenças.

  • Atribua cada fator \phi_k \in \Phi para cada cluster C_{\alpha(k)}
  • Construa os potenciais \psi(C_i)
  • Inicializa todas as mensagens com 1 (equivalente ao fator nulo)
  • Repita:
    • Selecione uma aresta (i,j) e faça a passagem de mensagem, atualizando \delta_{i \rightarrow j}
  • Calcule a crença de cada cluster \beta(C_i)

Alguns detalhes que o algoritmo não especifica é: quantas vezes iterar? Uma observação é que o algoritmo pode não convergir. Além disso, como escolher a ordem das arestas? Uma sugestão é o uso do Round-Robin.

Propriedades dos grafos cluster

Um grafo cluster deve possuir as seguintes propriedades.

Preservação de família. (family preservation) Para cada fator \phi_k \in \Phi, deve existir um cluster C tal que \mbox{Escopo}[\phi_k] \subseteq C, para que possamos atribuir o fator a algum cluster.

Caminho de Interseções. (tradução livre de Running Intersection) Essa propriedade diz que para cada par de vértices C_i e C_j no grafo e uma variável X \in C_i \cap C_j, deve existir um único caminho entre C_i e C_j tal que X \in S_{u,v}, para toda aresta (u, v) neste caminho.

Intuitivamente, isto serve para garantir que X possa ser “transportado” entre C_i e C_j. Por outro lado, a existência de mais de um caminho implica em ciclos, o que poderia causa uma retro-alimentação indesejável.

Uma forma alternativa de definir essa propriedade é através da decomposição em árvore, que diz que para cada variável X, o conjunto de vértices e arestas que contêm X, (\{C_k | X \in C_k\} e \{(i,j) | X \in S_{ij}\}, respectivamente), devem formar uma árvore.

Construção de um grafo cluster

Grafo cluster Bethe. é um tipo especial de grafo cluster bastante simples de construir e que satisfaz as duas propriedades definidas acima. Ele é composto de dois tipos de vértices:

Clusters grandes correspondendo ao domínio de cada fator, ou seja, C_{\phi_k} = \mbox{Escopo}(\phi_k) e clusters pequenos, formado por apenas uma variável, C_i = \{X_i\}.

Existe uma aresta (C_{\phi_k}, C_{i}) se X_i \in C_k.

É fácil ver que esse grafo satisfaz a propriedade dos caminho de interseções, pois ao isolarmos uma variável, o conjunto de vértices e arestas que a contém é uma estrela.

Conclusões

Não ficou muito claro como podemos usar esse segundo algoritmo para calcular uma distribuição de probabilidade condicionada a uma observação, por exemplo P(Y|E=e). Pelo que entendi, a crença de um cluster representa a distribuição marginalizada das variáveis contidas nele.

Assim, podemos inicialmente reduzir os fatores removendo as entradas tais que E \neq e e adicionar um cluster contendo Y. Executamos o algoritmo e fazemos a consulta da crença neste cluster.