Métodos de regiões de confiança para resolução do problema de quadrados mínimos: implementação e testes numéricos

Métodos de regiões de confiança para resolução do problema de quadrados mínimos: implementação e testes numéricos¹

J.L.C. Gardenghi^I; S.A. Santos^II

^IDepartamento de Ciência de Computação, IME - Instituto de Matemática e Estatística, USP - Universidade de São Paulo, 05508-090 São Paulo, SP, Brasil. [e-mail]

^IIDepartamento de Matemática Aplicada, IMECC - Instituto de Matemática, Estatística e Computação Científica, UNICAMP - Universidade Estadual de Campinas, 13083-859 Campinas, SP, Brasil. [e-mail]

RESUMO

O problema de quadrados mínimos possui várias aplicações no campo de otimização. No presente trabalho, abordamos duas estratégias para sua resolução: Levenberg-Marquardt e Gradientes Conjugados. Cada uma explora características próprias do problema, e ambas usam regiões de confiança para a globalização. Nossa contribuição está na implementação de ambos os métodos no CAS Maxima e na análise comparativa do desempenho desses métodos na resolução de uma família de problemas de quadrados mínimos da literatura.

Palavras-chave: Quadrados mínimos, regiões de confiança, implementação computacional.

ABSTRACT

The least squares problem has many applications in the field of optimization. In the present work, we use two strategies for its resolution: Levenberg-Marquardt and Conjugate Gradients. Each one exploits some problem features, and both are globalized by the trust-region strategy. Our contribution consists in the implementation of both methods using the CAS Maxima and in the comparative analysis of these methods in the resolution of a family of least squares problems from the literature.

Keywords: Least squares, trust region, computational implementation.

1. Introdução

Neste trabalho, consideramos o problema de quadrados mínimos

sendo F : Rⁿ - R^m a função residual, ao menos duas vezes continuamente di-ferenciável, F(x) o resíduo do problema em x e || · || a norma Euclidiana. No campo de otimização, quadrados mínimos tem várias aplicações, destacando-se seu uso para ajustar parâmetros de modelos frente a conjuntos de dados coletados em experimentos, e neste contexto m ^> n. Para um recente panorama do estado da arte deste problema, tanto do ponto de vista teórico quanto prático, veja [14] e suas referências.

Na próxima seção, apresentamos dois métodos para resolução de (1.1): o algoritmo de Levenberg-Marquardt (LMA) e o método de Gradientes Conjugados (MGC). Em seguida, destacamos os detalhes da implementação e os resultados dos experimentos numéricos que a validam e ilustram o desempenho comparativo das duas abordagens distintas. Finalizamos o texto com uma breve conclusão.

2. Métodos de Regiões de Confiança para Quadrados Mínimos

Métodos de regiões de confiança para minimização irrestrita, min f(x) s.a. x ∈ R^n, são estratégias iter ativas baseadas em modelos quadráticos para a função objetivo f em torno de um dado ponto [2]. Em geral, seja x_c ∈ Rⁿ o ponto corrente, define-se um modelo quadrático q_c para f em torno de x_c e o subproblema é formulado como

em que 0 < Δ_c ∈ R é o raio da região de confiança. A solução pc, exata ou aproximada, deste problema é usada para obter um ponto candidato x_n = x_c + p_c.

O bom desempenho de regiões de confiança está ligado ao controle do tamanho do raio, que determina a região em que o modelo é confiável para se aproximar a função objetivo. Nesta abordagem, vamos determinar o tamanho deste raio de acordo com o decréscimo que a solução de (2.1) causou no problema original que aproxima. Para isto, vamos avaliar a razão entre o decréscimo no modelo e na função objetivo, chamado quociente de redução relativa, dado por

No Algoritmo 1, página 71, apresentamos o algoritmo básico de regiões de confiança. Os parâmetros k, η₁ e η₂ são determinados experimentalmente, embora haja algumas escolhas cláss literatura [2, 11]. Sob esta ótica, escolhemos . Para resolver o subproblema (linha 2), implementamos LMA e MGC. Uma breve descrição destes dois métodos é apresentada a seguir.

2.1. Método de Levenberg-Marquardt

A principal característica deste método é atacar a estrutura particular do problema (1.1). Originalmente, esta estratégia foi proposta por Levenberg [6] e Marquardt [8], e foi considerada precursora da técnica de regiões de confiança (ver também [7]). More [9] propôs uma série de particularidades para implementação, que foram exploradas neste trabalho.

Dado um ponto x_c € R^m, LMA trabalha com a aproximação

em torno de x_c, sendo F_{c ∈ R^m} o resíduo em x_c e J_{c ∈ R}^mxn a jacobiana da função residual avaliada em x_c. Deste modo, o subproblema da região de confiança (2.1) é

um problema de quadrados mínimos lineares com passo restrito, cujas condições de otimalidade são dadas por

O sistema linear (2.3a) está bem definido somente se for definida positiva. Como é ao menos semidefinida positiva, então será definida positiva sempre que λ* > 0. Em particular, podemos ter dificuldades numéricas, por mal condicionamento ou por singularidade, quando λ* = 0. Neste caso, de acordo com (2.3b), temos que a solução de (2.1) pode ser interna, ou seja, ||p*|| < Δ_c. O resultado a seguir estabelece condições para que a solução de (2.1) seja interior à região de confiança.

Proposição 2.1.O problema (2.1) não possui solução global na fronteira se, e somente se, é definida positiva e o passo de Newton é interior, ou seja,.

Demonstração. Ver [4, Corolário 1].

Por outro lado, caso o problema possua solução na fronteira, é sempre possível encontrar _λ que satisfaça || p^*|| = Δc. Assim sendo, consideramos

em que † é a pseudoinversa de Moore-Penrose, e queremos encontrar λ* > 0 tal que

A equação (2.4) é chamada de equação secular. Por conseguinte, nosso problema se reduz a encontrar um zero de uma equação em uma variável. Para tal, usaremos o Método de Newton, com algumas salvaguardas que serão apresentadas na Seção 3. Notemos que _λ* atua como um regularizador, ou seja, ele ajusta o sistema (2.3a) para que sempre exista uma solução _p*. É justamente essa a conexão de LMA com regiões de confiança. Originalmente [6, 8], o parâmetro _λ* heurística. Com regiões de confiança, buscamos o multiplicador de Lagrange da restrição quadrática λ* > 0 de tal maneira que ||p*|| = Δ, salvo o caso particular em que a solução seja interna e então λ^* = 0

2.2. Método de Gradientes Conjugados

Neste método, vamos atacar o problema (1.1) como um problema geral de minimização irrestrita, e resolvê-lo usando a estratégia de região de confiança, considerando, a cada iteração, o subproblema de região de confiança (2.1), tomando

o modelo quadrático de f em torno de um ponto corrente x_c, com Hc G R^nxn, g_c ∈ Rⁿ e f_c ∈ R, respectivamente, a matriz hessiana, o vetor gradiente e a função objetivo f, todos avaliados em xc.

Para um estudo mais detalhado sobre a estratégia de gradientes conjugados, veja [5, Seção 3]. Para a resolução do subproblema, a idéia é produzir uma seqüência {d_i} de direções conjugadas e encontrar a solução para (2.1) iterativamente, por meio de p_i+1 = p_i + λidi. MGC seria uma estratégia η atural se H_c fosse uma matriz definida positiva, mas, na prática, nada sabemos sobre ela. É justamente neste ponto que as regiões de confiança tornam-se essenciais. As condições necessárias e suficientes de otimalidade de (2.1) são dadas por

Toint [13] e Steihaug [12] propuseram algumas modificações no MGC original de tal maneira que se ||g_c + H_cpi||for menor que dada tolerância, então temos uma solução interna, pois pⁱ satisfaz aproximadamente a condição (2.5a) com λ* = 0. Por outro lado, se ||p_i + 1|| > Δ_c, então basta encontrar τ de forma que ||p_i + τ_di|| ₌ Δ_c;. Por fim, se H_c não é definida positiva, então existe uma direção de curvatura negativa ou nula a partir de p_i. Assim sendo, se encontrarmos uma tal direção d_i, basta definir um passo que atinja a borda da região de confiança. Em outras palavras, se , computamos τ tal que ||p_i + τ_di|| ₌ Δ_c.

Deste modo, MGC atua corno GC clás sico quando H_c é definida positiva. Caso contrário, se for encontrada uma direção de curvatura negativa, determina-se um passo que atinja a borda da região de confiança. De qualquer forma, o passo tem seu tamanho restrito e essa também é uma condição de parada do método. No Algoritmo 2 descrevemos MGC.

3. Implementação Computacional

Nesta seção, apresentamos os detalhes de implementação do LMA que, por explorar as peculiaridades do problema (1.1), possui mais particularidades de implementação que o Algoritmo 2, que é um método para resolução do subproblema de regiões de confiança no caso de minimização irrestrita geral.

3.1. Álgebra matricial empregada na solução do subproblema

A equação (2.3a) pode ser vista como um sistema de equações normais para o problema linear de quadrados mínimos:

Tal equivalência permite resolver o problema sem avaliar o produto , reduzindo pela metade a quantidade de operações. Deste modo, se calcularmos a fatoração QR

teremos

A matriz que resulta de (3.2) possui uma estrutura particular: R_c é uma matriz triangular superior, enquanto é uma matriz escalar. Diante disto, é possível obter uma matriz triangular superior equivalente eliminando-se os η termos não-nulos da matriz com uma seqüência de n(n + 1)/2 rotações de Givens aplicadas em cada linha i de , simultaneamente a cada linha i + j de R_c, para j = 0, ..., (n - i), começando da linha n. Expressando tais rotações em uma matriz por meio de todos os sucessivos n(n + 1)/2 produtos, vem que

Disto, obtemos a matriz triangular superior Rλ, de modo que . Por conseguinte, o sistema (2.3a) torna-se menos custoso de ser resolvido para diferentes valores de λ. De fato, a fatoração QR de J_c é calculada uma única vez, no início da iteração interna, enquanto, para cada valor de λ, obtém-se Rλ por (3.3), que consiste em um procedimento mais barato que o cálculo da fatoração QR propriamente dita.

3.2. Raiz da equação secular

O problema (2.4) consiste em encontrar λ* > 0 tal que ||p(λ*)|| = Δ_c. Então, se definimos , estamos interessados em encontrar uma raiz para δ(λ). Todavia, uma análise mais apurada [4, Seção 3.2] nos mostra que a função δ pode possuir polos nos autovalores da matriz e pode não possuir zeros finitos. Em vista disto, uma estratégia mais estável é resolver a equação secular

que, por ser a recíproca de δ, não possui polos finitos e suas raízes acontecem nos valores negativos dos autovalores da matriz . Deste modo, para encontrar λ* > 0 que satisfaça (2.3a), avaliamos primeiro o caso particular em que λ* = 0 (solução interna). Se (2.3a) não possui solução neste caso, buscamos uma raiz positiva para a equação secular usando o método de Newton.

O método de Newton aplicado ao problema (3.4) gera uma seqüência {λ_j} tal que

Em contrapartida, seria interessante reduzir o custo do cálculo de λ_j+1 e ainda obter uma expressão numericamente mais estável. Isso é possível, conforme estabelece o resultado a seguir.

Teorema 3.1.Considere o problema (3.4) de encontrar a raiz da equação secular. Seja R_λa matriz triangular calculada em (3.3), tal que e e seja Δ_c o raio da região de confiança corrente. A seqüência {λ_j} gerada pelo método de Newton definida em (3.5) é equivalente a:

Demonstração. Ver [4, Lema 1].

Por fim, para que o método de Newton não convirja para uma raiz indesejada de φ(λ) (pois estamos interessados em raízes positivas), usamos uma salvaguarda de modo que, a cada iteração, impomos λ_j € [l_j, u_j], send o l₀ = 0 a estimativa para o menor autovalor de e u₀ obtido pelo quociente de Rayleigh

sendo µ_1;µ_n o menor e o maior autovalores de , respectivamente. Como zero é um limitante inferior para µ₁ e estamos procurando ||p|| = Δ_c, vem que

O intervalo da salvaguarda é atualizado a cada iteração da seguinte maneira:

se , então , caso contrário .

Para maiores detalhes, ver [4, Seção 4.2.1].

3.3. Quociente de redução relativa

Nesta subseção, nosso objetivo é eliminar a ocorrência de underflows e overflows no cálculo de (2.2). Usando LMA, podemos calcular o quociente como

em que F_n = F(x_c + p). De (2.3a), temos que

Algumas manipulações com (3.8) em (3.7) implicam que

Portanto, podemos reescrever o quociente (3.7) dividindo todos os termos por ||F_c||² e obtendo a seguinte relação:

Deste modo, o denominador de (3.9) será sempre não-negativo e não irá gerar overflows. Por outro lado, o numerador de (3.9) jamais gerará overflows a não ser que ||F_n || >> ||F_c||. Neste caso, temos que o modelo não foi reduzido com relação à função objetivo, portanto não é necessário calcular ρ(p).

3.4. Atualização do raio da região de confiança

Utilizamos a aproximação inicial , clássica da literatura [2]. O controle do raio da região consiste em três instâncias: manter, aumentar ou reduzir. Os casos em que cada uma é adotada são expostos no Algoritmo 1. Para aumentar o tamanho do raio, fazemos . Para reduzir o tamanho do raio, fazemos . A fim de definir um valor para v*, seguimos Moré [9]: definimos e encontramos o minimizador v* da parábola interpolante

Como m_q é suave e convexa, então temos que m'_q(v*) = 0, ou seja,

Pré-multiplicando (2.3a) avaliado em p* por (p*)^T, temos

donde segue que

Dividindo todos os terms de (3.10) por ||F_c||². multiplicando o numerador e o denominador por -1/2 e substituindo (3.11), que denotamos por α, obtemos

Pelos mesmos motivos expostos na seção anterior, (3.11) é uma expressão numericamente estável, e (3.12) é confiável e barata para encontrar v* e definir Δ_n.

4. Experimentos Numéricos

Ambos os métodos foram implementados no CAS Maxima [1, 15]. O Maxima é um sistema algébrico computacional livre de código aberto. Possui uma linguagem própria, interpretada, além de executar programas escritos em Lisp. Escolhemos o Maxima por ser um sistema gratuito e livre, além de não possuir ferramentas de otimização implementadas, como é o caso do Matlab^®. Os testes foram rodados num computador Intel Dual-Core T3400 2.16 GHz e sistema operacional Windows 7 Professional 64 bits. Os problemas-teste foram extraídos de [10] e estão numerados de acordo com a referência original, bem como as dimensões e os pontos iniciais adotados (ver detalhes em [4, 5]).

Como critérios de parada, adotamos a satisfação relativa das condições de primeira ordem, isto é, x_k é solução em LMA se e em MGC se . Além disso, em LMA impusemos parada se o resíduo fosse pequeno, ou seja, . Por fim, usamos critérios que implicam parada se a variação entre os iterandos fosse da ordem de 10^-12 e se a quantidade de iterações ultrapassasse 200 para LMA e 300 para MGC, todavia nos testes efetuados nenhum problema parou por tais critérios.

Os resultados dos experimentos são expostos na Tabela 1, enquanto os perfis de desempenho [3] na Figura 1. Neles, consideramos VF a função objetivo avaliada na solução, AH a quantidade de avaliações da matriz hessiana (para o MGC), AJ a quantidade de avaliações da matriz Jacobiana (para o LMA), ITE a quantidade de iterações externas, ITI a quantidade de iterações internas (isto é, para resolução dos subproblemas de regiões de confiança) e EI o esforço interno por problema, dado pela razão entre a quantidade de iterações internas e avaliações matriciais. A quantidade de avaliações da função objetivo é igual à quantidade de ITE acrescida de um. Os problemas marcados com * são aqueles que pararam pelo tamanho do resíduo em LMA; os demais pararam por satisfazer condições de primeira ordem.

Tabela 1: Quadro

comparativo entre os resultados das rotinas LMA e MGC.

Figura 1:
Perfis de desempenho dos algoritmos LMA e MGC.

As avaliações matriciais representam quantos dos pontos gerados foram aceitos. Deste modo, o número de pontos rejeitados, para cada problema, é dado pela diferença entre a quantidade de avaliações de funções e matriciais. Neste contexto, o esforço interno por problema representa quantas iterações internas foram despendidas, em média, para cada ponto aceito.

Acerca dos resultados apresentados, notamos que LMA é mais eficiente em avaliações de função, avaliações matriciais e quantidade de iterações internas, enquanto MGC, apenas em quantidade de iterações internas despendidas. Do ponto de vista de robustez, LMA se sobressai na quantidade de iterações internas, enquanto MGC em termos de avaliações, tanto de função quanto matriciais, e de iterações externas. Isto se deve ao fato de que LMA, por aproveitar melhor a estrutura do problema, trabalha mais internamente, consumindo menos iterações internas porém mais avaliações que MGC. Considerando-se, problema a problema, o número de iterações internas despendido para gerar a seqüência de pontos aceitos, o esforço interno de LMA é inferior ao de MGC. Com relação ao número de pontos rejeitados, no entanto, o resultado se inverte: MGC rejeita menos pontos que LMA.

5. Conclusão

O algoritmo LM A é uma proposta atraente para a resolução de (1.1), por ser um método mais barato, se comparado ao método de Newton. Em LMA, as informações de segunda ordem da função residual são substituídas pelo termo regularizador.

Em contrapartida, o método perde desempenho em problemas extremamente não lineadres ou se o resíduo for muito grande na solução_x*.

Em comparação com MGC, LMA aproveita a estrutura do problema (2.3a) e é mais eficiente em alguns quesitos. Todavia, exige muitas fatorações de matrizes, e trabalha mais internamente, o que pode ser ruim em alguns casos. Em vista disso, MGC é um método competitivo, mostrando-se mais robusto que LMA em parte dos quesitos analisados. Gradientes conjugados é bem conhecido por ser uma boa estratégia para resolução de sistemas lineares com matrizes definidas positivas, todavia aliado à estratégia de regiões de confiança pode ser uma alternativa interessante também para problemas não convexos.

Recebido em 28 Setembro 2012

Aceito em 14 Marco 2013

• [1] L.N. Andrade, Maxima: um completo programa de computação algêbrica, Revista do Professor de Matemática, n. 77, 2012.
• [2] A.R. Conn, N.I.L. Gould, Ph.L. Toint, "Trust-region Methods". SIAM, Philadelphia, 2000.
• [3] E.D. Dolan, J.J. More, Benchmarking optimization software with performance profiles, Mathematical Programming, v. 91, pp. 201-213, 2002.
• [4] J.L.C. Gardenghi, S.A. Santos, "Sistemas não-lineares via região de confiança: o algoritmo de Levenberg-Marquardt". Relatório de Pesquisa, 2011. Disponível em http://www.ime.unicamp.br/sites/default/files/rel_pesq/rp03-11.pdf Acesso em 05 mar. 2013.
• [5] J.L.C. Gardenghi, S.A. Santos, "Minimização irrestrita usando gradientes conjugados e regiões de confiança". Relatório de Pesquisa, 2012. Disponível em http://www.ime.unicamp.br/sites/default/files/rel_pesq/rp04-12.pdf Acesso em 05 mar. 2013.
• [6] K. Levenberg, A method for the solution of certain non-linear problems in least squares, The Quarterly of Applied Mathematics 2, pp. 164-168, 1944.
• [7] K. Madsen, An algorithm for the minimax solution of over determined systems of nonlinear equations. Journal of the Institute of Mathematics and its Applications, v. 16(3), pp. 321-328, 1975.
• [8] D.W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, SIAM Journal on Applied Mathematics 11, pp. 431-441, 1963.
• [9] J.J. More, The Levenberg-Marquardt algorithm: implementation and theory. Lecture Notes in Mathematics 630: Numerical Analysis, Springer-Verlag, New York, pp. 105-116, 1978.
• [10] J.J. More, B.S. Garbow, K.E. Hillstrom, Testing unconstrained optimization software, ACM Transactions on Mathematical Software, v. 7, pp. 17-41, 1981.
• [11] J. Nocedal, S.J. Wright, "Numerical Optimization". Springer, New York, 1999.
• [12] T. Steihaug, The conjugate gradient method and trust region in large scale optimization. SIAM Journal on Numerical Analysis, v. 20(3), pp. 626-637, 1983.
• [13] Ph.L. Toint, Towards an efficient sparsity exploiting Newton method for minimization, em "Sparse Matrix and Their Uses" (I. Duff, ed.), Academic Press, pp. 57-88, 1981.
• [14] Y. Yuan, Recent advances in numerical methods for nonlinear equations and nonlinear least-squares. Numerical Algebra, Control and Optimization, v. 1, n. 1, pp. 15-34, 2011.
• [15] http://maxima.sourceforge.net
  » link

1

Este trabalho contou com os apoios do CNPq (304032/2010-7, 106731/2011-4 e 1 51749/20116), FAPESP (2012/05725-0) e PRONEX-Otimização; Trabalho apresentado no XXXIV CNMAC.

Datas de Publicação

• Publicação nesta coleção
  28 Maio 2013
• Data do Fascículo
  Abr 2013

Histórico

• Recebido
  28 Set 2012
• Aceito
  14 Mar 2013