Probabilistic Graphical Models – Semana 5

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.

Os comentários estão fechados.

%d bloggers like this: