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.

Anúncios

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.