BrBRCEEn0101-74382002000100003

Otimização de entropia: implementação computacional dos princípios MaxEnt e MinxEnt

1.  Introdução O conceito de entropia foi introduzido na Ciência há mais de 150 anos, mas somente a partir de meados do Século XX é que difundiram-se suas aplicações por diversas áreas do conhecimento. Na raiz deste movimento, estiveram os trabalhos de Shannon (1948), que introduziu um conceito de entropia em teoria da informação e uma medida para quantificá-la, e os estudos de Jaynes (1957a,b) e Kullback (1959), que propuseram princípios de otimização da entropia segundo formulações distintas. Atualmente, diferentes áreas, como termodinâmica, probabilidade, estatística, pesquisa operacional, reconhecimento de padrões, economia, finanças, marketing, planejamento urbano e de transportes, dentre outras, vêm usando e desenvolvendo princípios de otimização da entropia (para diversos exemplos de aplicação em várias áreas, ver os livros de Kapur & Kesavan, 1992;Golan Judge & Miller, 1996; eFang, Rajasekera & Tsao, 1997).

O conceito de entropia de Shannon refere-se à incerteza de uma distribuição de probabilidade e a medida que propôs destinava-se a quantificar essa incerteza.

Formalmente, o princípio de Jaynes envolve a busca pela distribuição de probabilidade que maximiza a medida de Shannon, dado um conjunto de restrições lineares. Estas restrições informam características da distribuição procurada, como, por exemplo, sua média e sua variância. O princípio de Kullback, por sua vez, envolve a busca pela distribuição de probabilidade mais próxima de uma outra distribuição a priori, através da minimização de uma medida de divergência entre ambas, dado o mesmo conjunto de restrições. Tanto a medida de Shannon como a de Kullback são funções intrinsecamente não-lineares de probabilidades. Assim, os princípios de Jaynes e Kullback reduzem-se a problemas de programação não-linear cuja solução demanda um algoritmo de busca iterativa.

Em livro recente, Fang, Rajasekera & Tsao (1997, p. ix) apontam que, como esses (e outros) princípios de otimização de entropia foram repetidamente usados em várias áreas, muitos métodos para sua implementação (solução) foram sugeridos e utilizados. Entretanto, esses métodos carecem de uma formulação matemática mais rigorosa, que possa prover suporte aos praticantes interessados tanto no desenvolvimento de implementações computacionais eficientes como na adequada utilização dos princípios em aplicações. Por outro lado, embora esses mesmos autores apresentem um algoritmo para implementação dos princípios de Jaynes e Kullback, eles o fazem para o caso generalizado da medida de Kullback em que as distribuições são de freqüências e não de probabilidades (ver Fang, Rajasekera & Tsao (1997), pp. 51-56). No último caso, as restrições de não negatividade e de soma um para as probabilidades restringem de modo particular o conjunto de soluções viáveis, ao imporem limites inferiores e superiores para os coeficientes lineares das equações de restrição.

Este artigo objetiva apresentar em detalhe, com o auxílio de dois fluxogramas, a construção de um algoritmo para implementação computacional daqueles dois princípios de otimização da entropia no caso de distribuições de probabilidade discretas. O algoritmo aqui apresentado foi desenvolvido seguindo abordagem semelhante à proposta por Agmon, Alhassid & Levine (1979). No entanto, as descrições de ambos os princípios, do algoritmo e a caracterização das condições para existência de solução são feitas aqui de forma mais clara e direta, facilitando a implementação da metodologia através de diferentes linguagens computacionais. Como um subproduto do artigo, os autores disponibilizam (mediante solicitação por e-mail) funções escritas por eles em linguagem MatLabâque implementam os dois princípios e são de fácil uso.

Pretende-se, com isso, estimular a utilização de técnicas de otimização da entropia no contexto brasileiro. O público-alvo são pesquisadores e profissionais que utilizam intensivamente métodos quantitativos em suas áreas de atuação, em particular aqueles que trabalham com metodologias de P.O.

O artigo está organizado da seguinte forma. Na seção 2, o conceito de entropia na teoria da informação é revisto e breves considerações são tecidas sobre sua extensão a outras áreas onde o interesse se centra sobre distribuições de proporção e não de probabilidade. Nas seções 3 a 5, são discutidas as medidas de entropia de Shanon e de Kullback, bem como o formalismo dos princípios de Jaynes e Kullback. Na seção 6, a implementação computacional propriamente dita é descrita em detalhe e também são discutidas as condições para existência de solução factível dos problemas. Na Seção 7, um algoritmo de busca iterativa baseado no método de Newton é sugerido. Na Seção 8, são apresentadas as funções em linguagem MatLabâdesenvolvidas pelos autores para implementar o algoritmo proposto. Na seção 9, é apresentado um exemplo de aplicação do algoritmo em um problema de determinação do número de viagens, típico de planejamento de transportes. Finalmente, na seção 10, são tecidos alguns comentários conclusivos. Há também dois apêndices: o primeiro contém um fluxograma de todo o processamento do algoritmo proposto para implementar os dois princípios de otimização da entropia; o segundo descreve, com a ajuda de outro fluxograma, a implementação da Fase 1 do Método Simplex, usada para testar a existência de solução factível dos problemas de otimização considerados.

2.  Entropia Originária de estudos de termodinâmica, onde foi introduzida para caracterizar o grau de desordem de um sistema, a noção de entropia já foi objeto de muitas controvérsias e distintas formulações. O conceito de entropia adotado por Shannon (1948) foi responsável por aplicações de relevo em diversos campos de investigação científica, embora seu trabalho tenha se destacado mais pela medida de quantificação de entropia que propôs, cujas propriedades despertaram o interesse em outras áreas. Nesta seção, é feita uma breve apresentação do conceito de entropia na teoria da informação, visando preparar o terreno para, nas duas próximas seções, falar-se de medidas para sua quantificação.

Segundo Kapur & Kesavan (1992), o conceito de Shannon poderia ser chamado de entropia na teoria da informação e refere-se à incerteza de uma distribuição de probabilidade. Na verdade, o conceito de incerteza é mais geral, podendo-se falar, basicamente, em três tipos de incerteza: a incerteza determinística, em que não são conhecidos os estados que um sistema pode assumir; a incerteza entrópica, em que são conhecidos os estados possíveis, mas não as chances de ocorrência de cada um deles; e a incerteza probabilística, em que são conhecidos não só os estados possíveis mas também a distribuição de probabilidade para eles (todavia, não se pode determinar qual irá ocorrer com certeza).

A entropia na teoria da informação corresponde à incerteza probabilística associada a uma distribuição de probabilidade. Cada distribuição reflete um certo grau de incerteza e diferentes graus de incerteza estão associados a diferentes distribuições (embora diferentes distribuições possam refletir o mesmo grau de incerteza). De um modo geral, quanto mais "espalhada" a distribuição de probabilidade, maior incerteza ela irá refletir. Por exemplo, se alguém lança um dado de seis faces, sem saber se ele é viciado ou não, a probabilidade mais razoável a ser atribuída a cada resultado possível é 1/6, ou seja, representar a incerteza usando a distribuição uniforme. Esta atitude segue o conhecido princípio da razão insuficiente de Laplace, onde atribuir chances iguais aos eventos possíveis é a maneira mais razoável de alguém refletir sua ignorância (e sua incerteza) quanto às chances de ocorrência de cada evento. Por outro lado, provendo-se a informação de que o dado é viciado e que ele dá números maiores (menores) que a média (=3,5, no caso uniforme) mais freqüentemente, então a pessoa naturalmente irá assumir uma distribuição alternativa à uniforme para expressar sua incerteza. A Figura_1 ilustra graficamente essa situação no caso de distribuições contínuas de probabilidade.

Uma importante característica da entropia na teoria da informação, ou incerteza probabilística, é que ela está diretamente associada ao grau de similaridade entre as probabilidades de uma distribuição. Segundo Kapur & Kesavan (1992), este aspecto confere uma importante versatilidade à essa noção de entropia que lhe permite ser estendida e adaptada, enquanto conceito, à várias outras disciplinas. Entretanto, esta extensão/adaptação já foi questionada na literatura (Georgescu-Roegen, 1971) por não estar em consonância com a noção original de entropia em termodinâmica e nem com a própria noção de entropia na teoria da informação.

Sem pretender aprofundar essa discussão, o fato é que a medida introduzida por Shannon para quantificar entropia em teoria da informação também se presta a quantificar diversos conceitos de interesse em outras disciplinas. Se, ao invés de distribuição de probabilidades, trata-se de distribuição de proporções, como a distribuição intersetorial do produto industrial ou a distribuição espacial da ocupação residencial, é possível utilizar-se de modo interessante as medidas de entropia desenvolvidas em teoria da informação. Sob esta perspectiva, elas servem para medir igualdade, espalhamento, similaridade, diversidade, complexidade de sistemas e outros conceitos que aparecem em diversas áreas do conhecimento, ainda que tais conceitos não tenham uma relação direta com alguma noção clássica de entropia.

3. Medida de Entropia Shannon (1948) derivou uma medida para quantificar o grau de incerteza de uma distribuição de probabilidade. Denominando S a medida de entropia de Shannon, sua expressão formal para distribuições discretas de probabilidade é dada por:

onde pT = [p1, , pN] é a distribuição de probabilidade (o sobrescrito "T" representa transposição matricial). Esta medida é sempre não- negativa e assume seu valor máximo S(p) = ln N quando pT = [1/N, , 1/N] (i.e., a distribuição uniforme). Qualquer outra distribuição faz S ser menor do que ln N. Seu mínimo ocorre em S(p) = 0, situação onde há ausência de incerteza, quando então p é degenerada em uma das pis (i.e., uma pi = 1 e as demais iguais a zero).

Usando um método axiomático, Shannon derivou essa medida de modo que ela refletisse certas características desejadas. Posteriormente, outros matemáticos demonstraram que ela atende a outras propriedades de interesse adicional (Kapur & Kesavan, 1992, pp. 23-35). As propriedades de S mais relevantes para os fins deste trabalho são: S1. S(p1, p2, , pn) é uma função duas vezes diferenciável de p1, p2, , pN.

S2. S(p1, p2, , pn) é simétrica em relação à permutação de p1, p2, , pN.

S3. S(1/N, 1/N, , 1/N) é uma função monotonamente crescente de N.

S4. S é uma função estritamente côncava de p1, p2, , pN.

S1 é importante por permitir a aplicação de técnicas para maximização de funções diferenciáveis. S2 significa que as pis podem ter sua ordem invertida no cômputo de S que esta não se altera. S3 significa que a entropia da distribuição uniforme (máxima entropia possível) cresce quanto maior for o número de resultados possíveis N. Por último, S4 é de especial relevância, como se verá adiante, pois garante que S tenha um único máximo (global), mesmo quando sujeita a restrições lineares.

As propriedades de S permitem que ela também seja aplicada em diversos outros contextos. Quando, ao invés de probabilidades, as pis representarem proporções, isto é: , onde vi = valor da i-ésima parcela não negativa de uma soma, a medida de Shannon também pode ser aplicada, o que viabiliza aplicações em outras disciplinas. Por exemplo, S pode ser usada para medir o grau de igualdade (ou desigualdade) da distribuição de renda entre várias classes sociais, ou o grau de espalhamento das ocupações residenciais dentro de uma cidade.

4.  Medida de Entropia Cruzada Kullback (1959) introduziu outra importante medida em teoria da informação. A medida de entropia cruzada de Kullback é um caso particular de medidas de divergência direcionada e serve para medir a diferença entre duas distribuições de probabilidade. Sejam pT = [p1, , pN] e qT = [q1, , qN] duas distribuições quaisquer, e seja K a medida de Kullback. No caso de distribuições discretas de probabilidade, K é definida como:

Embora K não seja uma medida de incerteza de uma distribuição de probabilidade, ela serve aos mesmos propósitos que a medida de Shannon. A expressão (2) indica que K é uma medida de divergência ou diferença entre p e q. É fácil verificar que K(p:q) ¹ K(q:p), daí ela ser uma medida de divergência direcionada. Quanto maior a diferença/divergência entre p e q, maior será o valor de K. Dado o valor de N, seu máximo é atingido quando p é degenerada e q é a distribuição uniforme, situação em que K = ln N (ver expressão (3) abaixo). Quanto mais parecidas forem p e q, menor será K; no limite, se p = q, então K = 0.

Quando q é a distribuição uniforme, isto é, qT = uT = [1/N, , 1/N], então K também pode ser usada para medir incerteza ou entropia, pois neste caso:

Sendo ln N a entropia da distribuição uniforme (constante para um dado N), os graus de entropia de diferentes distribuições podem ser medidos e comparados entre si com base em suas divergências (medidas com a expressão (3)), em relação à distribuição uniforme. Adicionalmente, K também serve para indicar o grau de similaridade entre as entropias de duas distribuições, ainda que nenhuma delas seja a uniforme. Neste último caso, a medida de Kullback se presta a comparar diferentes distribuições com uma distribuição fixa qualquer.

K também apresenta várias propriedades atraentes. Dentre elas, destacam-se: K1. K(p:q) é uma função duas vezes diferenciável de p1, p2, , pN; K2. K(p:q) é simétrica em relação à permutação dos pares (p1, q1), , (pN, qN); K3. K é uma função estritamente convexa de p1, p2, , pN.

K4. K(p:q) ³ 0 (não negatividade); K5. K(p:q) = 0 se e apenas se p = q; As propriedades K1, K2 e K3 possuem implicações para K análogas às que S1, S2 e S4 apresentaram, respectivamente, para a medida de Shannon. As propriedades K4 e K5 são duas características de distâncias métricas (no entanto, K não é uma métrica pois não atende às propriedades de simetria e de desigualdade do triângulo que toda medida de distância tem de apresentar para ser uma métrica).

Da mesma forma que a medida de Shannon, a medida de Kullback se presta a estudos de entropia em que as distribuições se refiram a proporções e não a probabilidades.

5.  Otimização da Entropia Em teoria da informação, maximizar a entropia significa determinar a distribuição de probabilidade que represente o máximo de incerteza, dadas certas restrições. Ou seja, significa determinar a distribuição com maior grau de similaridade entre suas probabilidades, ou que seja mais parecida com a uniforme e diferindo dela apenas devido às restrições. Estas, por sua vez, refletem algum tipo de informação prévia sobre o fenômeno probabilístico de interesse, como, por exemplo, a média e a variância da distribuição que se quer determinar.

O princípio de maximizar a entropia (MaxEnt) através da medida de Shannon, dado um conjunto de restrições, foi introduzido por Jaynes (1957a,b).

Posteriormente, Kullback (1959) introduziu o princípio de minimização da entropia cruzada (MinxEnt), através do qual se procura minimizar a medida K, de divergência direcionada entre duas distribuições p e q, também sujeito a um conjunto de restrições. O princípio MinxEnt de Kullback generaliza o MaxEnt de Jaynes, pois permite que se incorpore, através de q, alguma informação a priori sobre a forma da distribuição de probabilidade procurada ao se otimizar a entropia. Quando qé a distribuição uniforme, o princípio MinxEnt se reduz ao princípio MaxEnt. Quando q é uma outra distribuição qualquer, o princípio MinxEnt envolve encontrar a distribuição p mais parecida com a priori q, ou a distribuição cuja entropia é a mais próxima da de q.

Nas subseções 5.1 e 5.2, o formalismo característico dos princípios MaxEnt e MinxEnt são introduzidos, buscando-se salientar as implicações para a implementação computacional de ambos que será apresentada na seção 6.

5.1  MaxEnt A aplicação do princípio MaxEnt pressupõe expressá-lo formalmente como um problema de otimização (doravante problema MaxEnt), da seguinte forma:

(onde s.a.. significa "sujeito a"). As funções gri (xi), r = 1, , M, são funções dos resultados possíveis xi, i = 1, , N. Note-se que o conjunto de restrições é formado por M + 1 restrições lineares e N restrições de não- negatividade, constituindo um típico problema de programação não-linear. Porém, a presença do termo ln pi em S implica que esta medida não está definida para valores negativos das pis, de modo que as N restrições de não-negatividade são não operantes (embora ln0 não seja definido, a medida S, no entanto, está definida para valores nulos das pis porque, quando x®0, lim (x.lnx)= 0). Isto simplifica o problema, permitindo que se aplique diretamente o método dos multiplicadores de Lagrange para otimização de funções não lineares com restrições de igualdade apenas.

A primeira das restrições lineares <formula/> é chamada de restrição natural, porque reflete a necessidade de que toda distribuição de probabilidade some um. As M restrições são denominadas de restrições de consistência. Nas aplicações em probabilidade, cada ar geralmente representa o momento de ordem r (o que implica fazer gri(xi) = xir ou gri(xi) = (xi ' m )r, com m representando a média da distribuição) ou então um momento característico da distribuição de probabilidade (sobre momentos característicos, ver por exemplo, Kapur & Kesavan, 1992, p. 359). Em várias aplicações onde as pis são tratadas como proporções, ar, xi e gri representam outro tipo de informação conhecida sobre o fenômeno de interesse (ver o exemplo da Seção 9, e também diversos outros nos livros de Kapur & Kesavan, 1992, e Fang & Tsao, 1997).

Usando-se o método do multiplicador de Lagrange, o problema MaxEnt (4) pode ser posto na seguinte forma irrestrita:

onde (l0' 1) e o vetorzT = [l1, ,lM] representam os M + 1 multiplicadores de Lagrange associados às M + 1 restrições. O multiplicadorl0 ' 1, ao invés de simplesmentel0, foi usado com a primeira restrição por conveniência matemática (Kapur & Kesavan, 1992, pp. 43-44), uma vez que permite simplificar as expressões apresentadas adiante. Aplicando a condição de primeira ordem para um extremo local, dada por ÑLs(p,l 0,z) = 0 (i.e., gradiente nulo ou conjunto das derivadas parciais iguais a zero), e manipulando algebricamente o sistema resultante, são obtidas as seguintes expressões:

A expressão (6) caracteriza a chamada distribuição de probabilidade MaxEnt. O sistema de M + N equações em M + Nincógnitas formado por (6) e (7) apresenta uma relação intrinsecamente não-linear entre as probabilidades pi e os multiplicadores de Lagrange lr, de modo que não é possível derivar uma solução analítica para pi e lr, simultaneamente, em função apenas dos elementos conhecidos ar e gri. Logo, a solução do sistema tem de ser obtida usando-se um algoritmo de busca iterativa. Note-se que um dos multiplicadores de Lagrange, l0, foi eliminado na manipulação algébrica (e, logo, do sistema de equações (6) e (7)), mas é simples verificar que ele pode ser obtido a partir dos demais multiplicadores segundo <formula/> Além disso, o problema (5) pode ser colocado ainda de uma outra forma.

Substituindo os pis do Lagrangeano Ls pelas expressões em (6) e realizando uma pequena manipulação algébrica, obtém-se:

Esta nova expressão apresenta como argumentos apenas os multiplicadores de Lagrange em zT = [l1, ,l M], embora seja intrinsecamente não-linear em relação a eles. Isto permite uma formulação dual irrestrita para o problema MaxEnt que consiste em Min Ls*(z). É possível mostrar que Ls* é uma função estritamente convexa dos multiplicadores de Lagrange z (Golan, Judge & Miller, 1996, pp. 25-26), o que assegura que o problema dual apresenta uma única solução (se houver solução). Posto deste modo, o problema tem sua dimensão reduzida, pois agora o sistema a ser resolvido é composto apenas por (7). Ao invés de se procurar iterativamente as N+ M variáveis em (pT,zT), basta a procura dos M multiplicadores de Lagrange emz. Uma vez determinada a solução ótima para estes, automaticamente ficam determinados os N valores para as pis, através da relação (6).

5.2  MinxEnt Analogamente, o princípio MinxEnt tem de ser formalizado como um problema de otimização (problema MinxEnt), da seguinte forma:

O que muda em relação ao problema anterior é que agora busca-se minimizar a Medida de Kullback. Como mencionado antes, esta mede a distância entre a distribuição definida pelas pis em relação a uma distribuição a priori definida pelas qis. O sistema de restrições lineares é o mesmo e tem o mesmo papel que no problema MaxEnt. Escrevendo <formula/>, verifica-se que esta medida também não está definida para valores negativos das pis (e das qis), de modo que a restrição pi ³ 0 também é não operante aqui. Se qi = 1/N (distribuição uniforme), (9) reduz-se ao problema MaxEnt, conforme se pode perceber pela expressão para K dada em (3).

Usando-se o método dos multiplicadores de Lagrange, o problema MinxEnt também pode ser posto na seguinte forma irrestrita:

onde, novamente, l0 ' 1 e z representam multiplicadores de Lagrange. De forma análoga, aplicando a condição de primeira ordem para um extremo local ÑLs (p,l0,z) = 0 (gradiente nulo ou conjunto das derivadas parciais iguais a zero), e manipulando algebricamente o sistema resultante, são obtidas as seguintes expressões:

A expressão (11) caracteriza a distribuição de probabilidade MinxEnt. O sistema (11)'(12) difere do sistema (6)'(7) apenas pela presença do termo qi multiplicando o expoente nos numeradores e denominadores das expressões para pi e ar. Aqui, da mesma forma, cai-se em um sistema de equações intrinsecamente não-lineares em pi e lr, de modo que é impossível derivar soluções analíticas para estas incógnitas em função dos elementos conhecidos ar e gri. Substituindo a expressão para pi em (11) no Lagrangeano Lk em (10) e manipulando, obtém-se:

Assim, o problema (9) admite uma formulação dual irrestrita, que consiste em Max Lk*(z). É possível mostrar que Lk* é uma função estritamente côncava dos multiplicadores de Lagrange z(Kapur & Kesavan, 1992, pp. 167-168), o que garante que o problema MinxEnt dual também tem uma única solução (se houver solução). Aqui também fica simplificado o problema ao se reduzir a dimensão M + N do espaço de busca para M. O sistema de equações não lineares a ser resolvido no problema dual é formado apenas por (12), que envolve apenas M equações nos M multiplicadores de Lagrange z. Achando-se estes, automaticamente determinam-se as N probabilidades pis por (11).

Os problemas duais para MaxEnt e MinxEnt admitem formulações análogas quando se redefine o último de modo que se torne também um problema de minimização, isto é:

A vantagem disso é que é possível desenvolver um único tipo de programa, voltado para um problema de minimização, que pode ser usado em ambos os casos MaxEnt e MinxEnt. O restante deste artigo apresenta a implementação computacional de um algoritmo para encontrar as soluções de ambos os problemas.

6.  Implementação Computacional Conforme visto na seção anterior, encontrar a distribuição MaxEnt envolve resolver um sistema de M equações intrinsecamente não-lineares em M incógnitas.

Para se achar a solução desse sistema, é necessário algum método de procura iterativa pela solução, o que requer implementação computacional mesmo para pequenas dimensões do problema. Antes, porém, de se discutir a montagem de um tal algoritmo, é preciso examinar as condições para existência de soluções factíveis dos problemas MaxEnt e MinxEnt. Embora os sistemas alvos sejam (7) e (12), que se apresentam não-lineares em z, o ponto de partida será o sistema de restrições lineares nas pis, porque estas é que são as incógnitas de interesse e porque é mais fácil derivar as condições factíveis para um sistema linear.

6.1  O Sistema de Restrições Considere-se, inicialmente, o sistema de restrições comum aos princípios MaxEnt e MinxEnt, apresentado no Quadro_6.1 Ele é composto por uma restrição natural de soma unitária das probabilidades e M restrições de consistência, formando um conjunto de M + 1 equações lineares. Existe um total de N incógnitas, constituídas pelas N pis. As restrições de não negatividade p1 ³ 0, p2 ³ 0,  , pN ³ 0 não aparecem no Quadro_6.1, pois são não operantes, como discutido anteriormente.

Conforme a teoria de sistemas lineares, o sistema do Quadro_6.1 admite as seguintes possibilidades: * se M +1 ³ N, então: ou existe um número infinito de soluções, ou existe uma única solução, ou não existe solução factível; * se M +1<N, então: ou existe um número infinito de soluções, ou não existe solução factível.

Se não existir solução factível, não é possível aplicar o princípio de Jaynes nem o de Kullback e não existe uma distribuição MaxEnt/MinxEnt. Se, por outro lado, existir solução factível, então pode-se afirmar que: a) quando M + 1 ³ N, com pelo menos N restrições linearmente independentes, maximizar a entropia é sinônimo de resolver o sistema de equações lineares formado pelas restrições (se M + 1 = N) ou por um subgrupo das restrições (se M + 1 > N); b) quando M + 1 < N, maximizar a entropia envolve escolher uma única solução dentre um número infinito de soluções.

A afirmação a) significa que não é necessário maximizar a medida de Shannon ou minimizar a de Kullback se o sistema de restrições lineares admite no máximo uma única solução factível. Neste caso, podem ser ignoradas S ou K e, simplesmente, obter-se a distribuição MaxEnt/MinxEnt usando-se procedimentos padronizados para resolução de sistemas de equações lineares.

A afirmação b), por sua vez, significa que encontrar a distribuição MaxEnt/ MinxEnt envolve escolher uma distribuição de probabilidade entre infinitas distribuições alternativas que atendem ao sistema de restrições. Otimizar a medida de Shannon ou a de Kullback, portanto, só faz sentido neste último caso.

6.2  Condições para existência de solução factível Como otimizar a entropia só faz sentido se M + 1 < N, este é o caso relevante a ser considerado a partir daqui. Entretanto, esta situação admite também a possibilidade de não existir solução factível para o sistema de restrições. Se for este o caso, não faz sentido aplicar um algoritmo de busca. É importante, então, estabelecer quais as condições precisas para a existência de solução do sistema de restrições lineares quando M + 1 < N.

Consideremos a seguinte representação matricial do sistema de restrições:

onde l = vetor de dimensão N com todos os elementos iguais a 1, G= {gri(xi)} é a matriz de coeficientes do subsistema de restrições de consistência, pT = [p1, ,pN] e aT = [a1, ,aM]. As condições para existência de solução factível para o sistema (14) quando M + 1 < N são as seguintes: a) a matriz [1 GT] T tem de apresentar posto cheio (= M + 1), o que implica que as restrições têm de ser linearmente independentes entre si; b) Se M = 1, isto é, se há uma única restrição de consistência, de modo que G = [g11, , g1N], então:

A condição a) é óbvia. A condição b) tem uma justificativa intuitiva. Como o coeficiente a1 é uma média ponderada (pelas probabilidades pis) das g1is, ele não pode ser inferior à menor nem superior à maior delas.

Os sinais de desigualdade estrita em (15) são relevantes porque: * se a1 ® min{G}, entãol1®¥ e a distribuição MaxEnt/MinxEnt degenera-se em pT = (1,0, ,0) com S = 0 (não há incerteza com um único resultado possível); * sea1 ® max{G}, entãol 1®- ¥ e a distribuição MaxEnt/MinxEnt degenera-se em pT = (0, ,1)com S = 0; * sea1 = min{G} = max{G}, então G é linearmente dependente de 1T, ou seja, a única restrição de consistência é linearmente dependente da restrição natural. Logo, a condição b) inclui a condição a).

c) Se M > 1, então a condição (15) tem de ser generalizada, o que é obtido pré-multiplicando-se ambos os lados do subsistema de restrições de consistência em (14) por um vetor de constantes cT = [c1, ,cM], tal que cTGp = cTa. Usando-se o mesmo argumento intuitivo, tem-se a condição generalizada:

onde cTG é um vetor-linha (1´ N) e cTa é um escalar (para uma prova formal de (16), ver Alhassid, Agmon & Levine, 1978).

Note-se aqui que, fazendo G(r) a r-ésima linha de G, a generalização (16) não envolve simplesmente estipular que min{G(r)} < ar < max{G(r)} para todo r=1, ,M, porque embora esta condição seja necessária, ela não é suficiente. Uma razão para isso pode ser pensada quando se tem, por exemplo, duas restrições de momento em que <formula/> (média) e [/img/fbpe/pope/ v22n1/a03img22.gif] (2º momento não centrado). Se o problema for definido de modo que a1 esteja próximo de min i {g1i} e a2, em contraposição, esteja próximo de max i {g2i}, então é de se esperar que o sistema não seja factível.

Nesta situação, tende a haver uma inconsistência, pois a média é muito baixa e o 2º momento muito alto. Conseqüentemente, a faixa admissível de valores para os ars tem de ser mais estreita que a sugerida por (15).

Da mesma forma que antes, os sinais de desigualdade estrita em (16) são relevantes porque: * se cTa ® min{cTG}, entãol r ® ¥ (r = 0,1, ,M), a distribuição MaxEnt/MinxEnt é degenerada em p = (1,0, ,0) com S = 0; * se cTa ® max{cTG}, entãol r ® - ¥ (r = 0,1, ,M), a distribuição MaxEnt/MinxEnt é degenerada em pT = (0, ,1) comS = 0; * se cTa = min{cTG} = max{cTG}, então G é linearmente dependente de 1T, ou seja, as restrições de consistência são linearmente dependentes da restrição natural. A condição (16), portanto, inclui as duas anteriores.

A questão agora é: como verificar se o sistema (14) satisfaz (16)? Agmon, Alhassid & Levine (1979) sugerem que se divida o problema em 2 partes: i) verificar se [1 GT] T tem posto cheio (= M + 1) usando-se o procedimento Gram-Schmidt de ortogonalização de matrizes e; ii) verificar se (14) tem solução factível usando-se o algoritmo da Fase 1 do Método Simplex de programação linear.

Na implementação computacional aqui desenvolvida em linguagem MatLabâ, i) pode ser facilmente verificada usando-se a função RANK, que indica o posto de uma matriz. Para verificar ii), entretanto, foi necessário o desenvolvimento de uma função específica e esse é o assunto da próxima subseção.

6.3  A Fase 1 do Método Simplex O objetivo do Método Simplex é encontrar uma solução ótima para o seguinte problema padrão de programação linear:

onde: d = vetor de coeficientes lineares da função objetivo;

x = vetor das N variáveis de interesse;

s= vetor das M variáveis de folga; A= matriz (M´N) dos coeficientes de x nas M restrições e I= matriz (M´M) de identidade.

Assume-se no problema padrão que b ³ 0. O sistema de restrições Ax ' Is = b possui, portanto, M restrições e N + M incógnitas e representa o conjunto de soluções factíveis para o PLP.

Para ser inicializado, o algoritmo do Método Simplex precisa dispor de uma solução básica inicial (SBI) do sistema de restrições. Uma SBI é um determinado valor de (xT,sT)no qual existem Nelementos nulos. A chamada Fase 1 do Método Simplex corresponde a essa etapa de busca de uma SBI para se inicializar o algoritmo. Uma importante característica da Fase 1 é que ela ou acha uma solução básica inicial ou detecta que o PLP não tem solução factível (Solow, 1984). Daí ser útil para indicar se há ou não solução para os problemas MaxEnt e MinxEnt.

Operacionalmente, a Fase 1 usa o mesmo Método Simplex, mas para resolver um PLP modificado (PLPM), no qual são usadas variáveis artificiais y, como segue:

Na prática, as restrições lineares do problema MaxEnt e MinxEnt são de igualdade, o que implica que s= 0. Então, definindo-se [/img/fbpe/pope/v22n1/ a03img25.gif] como a solução ótima do PLPM, isto é, da Fase 1, podem haver duas situações: a)<formula/>= 0, logo existe uma SBI para o PLP original ou b)<formula/>¹ 0, logo não existe solução factível para o PLP original ou para o seu sistema de restrições.

Uma vez detectado que não há solução factível para o sistema de restrições, o programa de computador tem de ser abortado. O algoritmo detalhado para implementação da Fase 1 do Método Simplex está apresentado no Apêndice_II. Na próxima seção, será apresentado o algoritmo desenvolvido para buscar iterativamente a solução dos problemas MaxEnt e MinxEnt.

7.  Método de Newton A obtenção de uma solução para os problemas MaxEnt e MinxEnt duais demanda a utilização de um algoritmo de procura iterativa. Ambos são problemas de otimização não-linear com restrições de igualdade apenas, para os quais existem vários métodos de procura disponíveis. Na presente implementação, foi adotado o conhecido Método de Newton, devido às suas vantagens de convergência rápida e certa (quando a função-objetivo é estritamente convexa), além de ser de fácil programação.

Na implementação deste método, é realizado o cálculo do gradiente, do Hessiano e da inversa do Hessiano da função objetivo a cada iteração. O cálculo da inversa constitui uma desvantagem em problemas maiores, mas como a dimensão do Hessiano nos problemas MaxEnt e MinxEnt é determinada pelo número de restrições M, raramente este número será grande. Obviamente, melhoras podem ser introduzidas com outros enfoques, como por exemplo os métodos de quase-Newton (Bertsekas, 1995).

Considere-se, inicialmente, as representações primal e dual para o problema MaxEnt, escritas em forma vetorial e matricial no Quadro_7.1:

onde <formula/> e zT = [l1, , lM]. A solução do primal é dada por:

A implementação do Método de Newton para o problema dual envolve a busca iterativa de z* = solução ótima do problema através do seguinte algoritmo:

onde k = índice de iteração, dk = direção de Newton, hk = ÑLk* = gradiente de Lk* e Hk = Ñ2Lk* = matriz hessiana de Lk*. A direção de Newton dk usa informações referentes às primeiras e segundas derivadas em relação ao vetor de multiplicadores de Lagrange z, que estão embutidas no gradiente e na matriz hessiana. Como ambos são conhecidos analiticamente, podem ser calculados diretamente, isto é, sem o uso de aproximações, da seguinte forma:

onde pk é calculado a cada iteração segundo (19). Note-se que a matriz hessiana Hk nada mais é do que a matriz M´M de variância-covariância das colunas de Gcalculada sob a distribuição pk, sendo, portanto, definida positiva sempre.

O algoritmo completo de busca pela solução dual de MaxEnt é formado pelas expressões em (20) e (21). Sua implementação demanda ainda resolver a questão da inicialização e a definição do critério de parada. Em ambos os casos, existem procedimentos alternativos que podem ser utilizados. Na inicialização, é preciso definir-se o vetor inicial de multiplicadores de Lagrange z0T = [l10, ,lM0]. Um cuidado a ser tomado é o de que, se z0 for definido na região assintótica, então isto pode provocar uma quase singularidade da matriz H1. Um procedimento mais seguro pode ser definir-se z0 = 0. De fato, nos diversos testes realizados com a função MAXENT escrita em MatLabâ, este valor de inicialização z0 foi o que apresentou os resultados corretos com maior freqüência. Por sua vez, o critério de parada também pode ser definido de diversas formas. A que foi adotada na função MAXENT é baseada na norma do gradiente e em um número máximo de iterações, de modo que o algoritmo pára quando || hk || < 10-5 ou quando o número de iterações exceder 30.