Algoritmo de Levenberg-Marquardt - Levenberg–Marquardt algorithm

Em matemática e computação, o algoritmo de Levenberg-Marquardt ( LMA ou apenas LM ), também conhecido como método de mínimos quadrados amortecidos ( DLS ), é usado para resolver problemas de quadrados mínimos não lineares . Esses problemas de minimização surgem especialmente no ajuste de curvas de mínimos quadrados .

O LMA é usado em muitos aplicativos de software para resolver problemas genéricos de ajuste de curvas. No entanto, como acontece com muitos algoritmos de ajuste, o LMA encontra apenas um mínimo local , que não é necessariamente o mínimo global . O LMA interpola entre o algoritmo de Gauss-Newton (GNA) e o método de descida do gradiente . O LMA é mais robusto do que o GNA, o que significa que em muitos casos ele encontra uma solução mesmo que comece muito longe do mínimo final. Para funções bem comportadas e parâmetros iniciais razoáveis, o LMA tende a ser mais lento do que o GNA. O LMA também pode ser visto como Gauss-Newton usando uma abordagem de região de confiança .

O algoritmo foi publicado pela primeira vez em 1944 por Kenneth Levenberg , enquanto trabalhava no Frankford Army Arsenal . Foi redescoberto em 1963 por Donald Marquardt , que trabalhou como estatístico na DuPont , e independentemente por Girard, Wynne e Morrison.

O problema

A principal aplicação do algoritmo de Levenberg-Marquardt é no problema de ajuste de curva de mínimos quadrados: dado um conjunto de pares empíricos de variáveis ​​independentes e dependentes, encontre os parâmetros da curva do modelo de modo que a soma dos quadrados dos desvios seja minimizada : ${\ displaystyle m}$${\ displaystyle \ left (x_ {i}, y_ {i} \ right)}$${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle f \ left (x, {\ boldsymbol {\ beta}} \ right)}$${\ displaystyle S \ left ({\ boldsymbol {\ beta}} \ right)}$

${\ displaystyle {\ hat {\ boldsymbol {\ beta}}} \ in \ operatorname {argmin} \ limits _ {\ boldsymbol {\ beta}} S \ left ({\ boldsymbol {\ beta}} \ right) \ equiv \ operatorname {argmin} \ limits _ {\ boldsymbol {\ beta}} \ sum _ {i = 1} ^ {m} \ left [y_ {i} -f \ left (x_ {i}, {\ boldsymbol {\ beta}} \ right) \ right] ^ {2},}$ que é considerado não vazio.

A solução

Como outros algoritmos de minimização numérica, o algoritmo de Levenberg – Marquardt é um procedimento iterativo . Para iniciar uma minimização, o usuário deve fornecer uma estimativa inicial para o vetor de parâmetro . Em casos com apenas um mínimo, uma estimativa padrão não informada como funcionará bem; em casos com múltiplos mínimos , o algoritmo converge para o mínimo global apenas se a estimativa inicial já estiver um pouco próxima da solução final. ${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle {\ boldsymbol {\ beta}} ^ {\ text {T}} = {\ begin {pmatrix} 1, \ 1, \ \ dots, \ 1 \ end {pmatrix}}}$

Em cada etapa de iteração, o vetor de parâmetros é substituído por uma nova estimativa . Para determinar , a função é aproximada por sua linearização : ${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle {\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}}}$${\ displaystyle {\ boldsymbol {\ delta}}}$${\ displaystyle f \ left (x_ {i}, {\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right)}$

${\ displaystyle f \ left (x_ {i}, {\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right) \ approx f \ left (x_ {i}, {\ boldsymbol {\ beta} } \ right) + \ mathbf {J} _ {i} {\ boldsymbol {\ delta}},}$

Onde

${\ displaystyle \ mathbf {J} _ {i} = {\ frac {\ partial f \ left (x_ {i}, {\ boldsymbol {\ beta}} \ right)} {\ partial {\ boldsymbol {\ beta} }}}}$

é o gradiente (vetor linha neste caso) de em relação a . ${\ displaystyle f}$${\ displaystyle {\ boldsymbol {\ beta}}}$

A soma dos desvios quadrados tem seu mínimo em um gradiente zero em relação a . A aproximação de primeira ordem acima dá ${\ displaystyle S \ left ({\ boldsymbol {\ beta}} \ right)}$${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle f \ left (x_ {i}, {\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right)}$

${\ displaystyle S \ left ({\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right) \ approx \ sum _ {i = 1} ^ {m} \ left [y_ {i} -f \ left (x_ {i}, {\ boldsymbol {\ beta}} \ right) - \ mathbf {J} _ {i} {\ boldsymbol {\ delta}} \ right] ^ {2},}$

ou em notação vetorial,

${\ displaystyle {\ begin {alinhados} S \ left ({\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right) & \ approx \ left \ | \ mathbf {y} - \ mathbf {f } \ left ({\ boldsymbol {\ beta}} \ right) - \ mathbf {J} {\ boldsymbol {\ delta}} \ right \ | ^ {2} \\ & = \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) - \ mathbf {J} {\ boldsymbol {\ delta}} \ right] ^ {\ mathrm {T}} \ left [\ mathbf { y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) - \ mathbf {J} {\ boldsymbol {\ delta}} \ right] \\ & = \ left [\ mathbf {y } - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] ^ {\ mathrm {T}} \ left [\ mathbf {y} - \ mathbf {f} \ left ({ \ boldsymbol {\ beta}} \ right) \ right] - \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] ^ {\ mathrm { T}} \ mathbf {J} {\ boldsymbol {\ delta}} - \ left (\ mathbf {J} {\ boldsymbol {\ delta}} \ right) ^ {\ mathrm {T}} \ left [\ mathbf { y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] + {\ boldsymbol {\ delta}} ^ {\ mathrm {T}} \ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} {\ boldsymbol {\ delta}} \\ & = \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] ^ {\ mathrm {T}} \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] -2 \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right] ^ {\ mathrm { T}} \ mathbf {J} {\ boldsymbol {\ delta}} + {\ boldsymbol {\ delta}} ^ {\ mathrm {T}} \ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J } {\ boldsymbol {\ delta}}. \ end {alinhado}}}$

Tirar a derivada de em relação a e definir o resultado para zero dá ${\ displaystyle S \ left ({\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}} \ right)}$${\ displaystyle {\ boldsymbol {\ delta}}}$

${\ displaystyle \ left (\ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} \ right) {\ boldsymbol {\ delta}} = \ mathbf {J} ^ {\ mathrm {T}} \ esquerda [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right],}$

onde é a matriz Jacobiana , cuja -ésima linha é igual , e onde e são vetores com -ésima componente e respectivamente. A expressão acima obtida para vem do método de Gauss-Newton. A matriz Jacobiana conforme definida acima não é (em geral) uma matriz quadrada, mas uma matriz retangular de tamanho , onde é o número de parâmetros (tamanho do vetor ). A multiplicação da matriz produz a matriz quadrada necessária e o produto matriz-vetor no lado direito produz um vetor de tamanho . O resultado é um conjunto de equações lineares que podem ser resolvidas . ${\ displaystyle \ mathbf {J}}$${\ displaystyle i}$${\ displaystyle \ mathbf {J} _ {i}}$${\ displaystyle \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right)}$${\ displaystyle \ mathbf {y}}$${\ displaystyle i}$${\ displaystyle f \ left (x_ {i}, {\ boldsymbol {\ beta}} \ right)}$${\ displaystyle y_ {i}}$${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle m \ times n}$${\ displaystyle n}$${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle \ left (\ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} \ right)}$${\ displaystyle n \ times n}$${\ displaystyle n}$${\ displaystyle n}$${\ displaystyle {\ boldsymbol {\ delta}}}$

A contribuição de Levenberg é substituir esta equação por uma "versão amortecida":

${\ displaystyle \ left (\ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} + \ lambda \ mathbf {I} \ right) {\ boldsymbol {\ delta}} = \ mathbf {J} ^ {\ mathrm {T}} \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right],}$

onde é a matriz identidade, dando como incremento ao vetor de parâmetros estimado . ${\ displaystyle \ mathbf {I}}$${\ displaystyle {\ boldsymbol {\ delta}}}$${\ displaystyle {\ boldsymbol {\ beta}}}$

O fator de amortecimento (não negativo) é ajustado a cada iteração. Se a redução de for rápida, um valor menor pode ser usado, aproximando o algoritmo do algoritmo de Gauss-Newton , enquanto se uma iteração der redução insuficiente no residual, pode ser aumentado, dando um passo mais perto da direção gradiente descendente. Observe que o gradiente de em relação a é igual . Portanto, para valores grandes de , o passo será dado aproximadamente na direção oposta ao gradiente. Se o comprimento da etapa calculada ou a redução da soma dos quadrados do último vetor de parâmetro cair abaixo dos limites predefinidos, a iteração para e o último vetor de parâmetro é considerado a solução. ${\ displaystyle \ lambda}$${\ displaystyle S}$${\ displaystyle \ lambda}$${\ displaystyle S}$${\ displaystyle {\ boldsymbol {\ beta}}}$${\ displaystyle -2 \ left (\ mathbf {J} ^ {\ mathrm {T}} \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ direita] \ direita) ^ {\ mathrm {T}}}$${\ displaystyle \ lambda}$${\ displaystyle {\ boldsymbol {\ delta}}}$${\ displaystyle {\ boldsymbol {\ beta}} + {\ boldsymbol {\ delta}}}$${\ displaystyle {\ boldsymbol {\ beta}}}$

O algoritmo de Levenberg tem a desvantagem de que, se o valor do fator de amortecimento for grande, a inversão não será usada. Fletcher forneceu a ideia de que podemos dimensionar cada componente do gradiente de acordo com a curvatura, de modo que haja um movimento maior ao longo das direções em que o gradiente é menor. Isso evita convergência lenta na direção de pequeno gradiente. Portanto, Fletcher em seu artigo de 1971 Uma sub-rotina de Marquardt modificada para mínimos quadrados não lineares substituiu a matriz de identidade pela matriz diagonal que consiste nos elementos diagonais de , tornando assim a escala da solução invariante: ${\ displaystyle \ lambda}$${\ displaystyle \ mathbf {J} ^ {\ text {T}} \ mathbf {J} + \ lambda \ mathbf {I}}$${\ displaystyle \ mathbf {I}}$${\ displaystyle \ mathbf {J} ^ {\ text {T}} \ mathbf {J}}$

${\ displaystyle \ left [\ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} + \ lambda \ operatorname {diag} \ left (\ mathbf {J} ^ {\ mathrm {T}} \ mathbf {J} \ right) \ right] {\ boldsymbol {\ delta}} = \ mathbf {J} ^ {\ mathrm {T}} \ left [\ mathbf {y} - \ mathbf {f} \ left ({\ boldsymbol {\ beta}} \ right) \ right].}$

Um fator de amortecimento semelhante aparece na regularização de Tikhonov , que é usada para resolver problemas lineares mal colocados , bem como na regressão de crista , uma técnica de estimativa em estatística .

Escolha do parâmetro de amortecimento

Vários argumentos mais ou menos heurísticos foram apresentados para a melhor escolha para o parâmetro de amortecimento . Existem argumentos teóricos que mostram porque algumas dessas escolhas garantem a convergência local do algoritmo; entretanto, essas escolhas podem fazer com que a convergência global do algoritmo sofra das propriedades indesejáveis ​​da descida mais íngreme , em particular, convergência muito lenta próxima ao ótimo. ${\ displaystyle \ lambda}$

Os valores absolutos de qualquer escolha dependem de quão bem dimensionado é o problema inicial. Marquardt recomendou começar com um valor e um fator . Inicialmente definindo e computando a soma residual dos quadrados após uma etapa do ponto inicial com o fator de amortecimento de e, em seguida, com . Se ambos forem piores do que o ponto inicial, então o amortecimento é aumentado pela multiplicação sucessiva de até que um ponto melhor seja encontrado com um novo fator de amortecimento de para alguns . ${\ displaystyle \ lambda _ {0}}$${\ displaystyle \ nu> 1}$${\ displaystyle \ lambda = \ lambda _ {0}}$${\ displaystyle S \ left ({\ boldsymbol {\ beta}} \ right)}$${\ displaystyle \ lambda = \ lambda _ {0}}$${\ displaystyle \ lambda _ {0} / \ nu}$${\ displaystyle \ nu}$${\ displaystyle \ lambda _ {0} \ nu ^ {k}}$${\ displaystyle k}$

Se o uso do fator de amortecimento resultar em uma redução do resíduo ao quadrado, então este é considerado o novo valor de (e a nova localização ótima é considerada como aquela obtida com este fator de amortecimento) e o processo continua; se o uso resultou em um resíduo pior, mas o uso resultou em um resíduo melhor, então é deixado inalterado e o novo ótimo é tomado como o valor obtido com o fator de amortecimento. ${\ displaystyle \ lambda / \ nu}$${\ displaystyle \ lambda}$${\ displaystyle \ lambda / \ nu}$${\ displaystyle \ lambda}$${\ displaystyle \ lambda}$${\ displaystyle \ lambda}$

Uma estratégia eficaz para o controle do parâmetro de amortecimento, chamada gratificação retardada , consiste em aumentar o parâmetro em uma pequena quantidade para cada degrau de subida e diminuindo em uma grande quantidade para cada degrau de descida. A ideia por trás dessa estratégia é evitar a descida muito rápida no início da otimização, restringindo, portanto, as etapas disponíveis em iterações futuras e, portanto, desacelerando a convergência. Um aumento por um fator de 2 e uma diminuição por um fator de 3 tem se mostrado eficaz na maioria dos casos, enquanto para problemas grandes valores mais extremos podem funcionar melhor, com um aumento de um fator de 1,5 e uma diminuição de um fator de 5.

Aceleração geodésica

Ao interpretar a etapa de Levenberg-Marquardt como a velocidade ao longo de um caminho geodésico no espaço de parâmetros, é possível melhorar o método adicionando um termo de segunda ordem que explica a aceleração ao longo do geodésico ${\ displaystyle {\ boldsymbol {v}} _ {k}}$${\ displaystyle {\ boldsymbol {a}} _ {k}}$

${\ displaystyle {\ boldsymbol {v}} _ {k} + {\ frac {1} {2}} {\ boldsymbol {a}} _ {k}}$

onde está a solução de ${\ displaystyle {\ boldsymbol {a}} _ {k}}$

${\ displaystyle {\ boldsymbol {J}} _ {k} {\ boldsymbol {a}} _ {k} = - f_ {vv}.}$

Uma vez que este termo de aceleração geodésica depende apenas da derivada direcional ao longo da direção da velocidade , ele não requer o cálculo de toda a matriz derivada de segunda ordem, exigindo apenas uma pequena sobrecarga em termos de custo de computação. Uma vez que a derivada de segunda ordem pode ser uma expressão bastante complexa, pode ser conveniente substituí-la por uma aproximação de diferença finita${\ displaystyle f_ {vv} = \ sum _ {\ mu \ nu} v _ {\ mu} v _ {\ nu} \ parcial _ {\ mu} \ parcial _ {\ nu} f ({\ boldsymbol {x}} )}$${\ displaystyle {\ boldsymbol {v}}}$

${\ displaystyle {\ begin {alinhados} f_ {vv} ^ {i} & \ approx {\ frac {f_ {i} ({\ boldsymbol {x}} + h {\ boldsymbol {\ delta}}) - 2f_ { i} ({\ boldsymbol {x}}) + f_ {i} ({\ boldsymbol {x}} - h {\ boldsymbol {\ delta}})} {h ^ {2}}} \\ & = {\ frac {2} {h}} \ left ({\ frac {f_ {i} ({\ boldsymbol {x}} + h {\ boldsymbol {\ delta}}) - f_ {i} ({\ boldsymbol {x} })} {h}} - {\ boldsymbol {J}} _ {i} {\ boldsymbol {\ delta}} \ right) \ end {alinhado}}}$

onde e já foram calculados pelo algoritmo, exigindo, portanto, apenas uma avaliação de função adicional para calcular . A escolha do passo de diferença finita pode afetar a estabilidade do algoritmo, e um valor em torno de 0,1 é geralmente razoável em geral. ${\ displaystyle f ({\ boldsymbol {x}})}$${\ displaystyle {\ boldsymbol {J}}}$${\ displaystyle f ({\ boldsymbol {x}} + h {\ boldsymbol {\ delta}})}$${\ displaystyle h}$

Uma vez que a aceleração pode apontar na direção oposta à velocidade, para evitar que ela estale o método caso o amortecimento seja muito pequeno, um critério adicional na aceleração é adicionado a fim de aceitar uma etapa, exigindo que

${\ displaystyle {\ frac {2 \ left \ | {\ boldsymbol {a}} _ {k} \ right \ |} {\ left \ | {\ boldsymbol {v}} _ {k} \ right \ |}} \ leq \ alpha}$

onde geralmente é fixado em um valor menor que 1, com valores menores para problemas mais difíceis. ${\ displaystyle \ alpha}$

A adição de um termo de aceleração geodésica pode permitir um aumento significativo na velocidade de convergência e é especialmente útil quando o algoritmo está se movendo através de desfiladeiros estreitos na paisagem da função objetivo, onde os passos permitidos são menores e a maior precisão devido à segunda ordem termo dá melhorias significativas.

Exemplo

Mal ajustado

Melhor ajuste

Melhor ajuste

Neste exemplo, tentamos ajustar a função usando o algoritmo Levenberg – Marquardt implementado no GNU Octave como a função leasqr . Os gráficos mostram progressivamente melhor encaixe para os parâmetros , utilizados na curva inicial. Somente quando os parâmetros no último gráfico são escolhidos mais próximos do original, as curvas se ajustam exatamente. Esta equação é um exemplo de condições iniciais muito sensíveis para o algoritmo de Levenberg-Marquardt. Uma razão para esta sensibilidade é a existência de múltiplos mínimos - a função tem mínimos no valor do parâmetro e . ${\ displaystyle y = a \ cos \ left (bX \ right) + b \ sin \ left (aX \ right)}$${\ displaystyle a = 100}$${\ displaystyle b = 102}$${\ displaystyle \ cos \ left (\ beta x \ right)}$${\ displaystyle {\ hat {\ beta}}}$${\ displaystyle {\ hat {\ beta}} + 2n \ pi}$

Veja também

• Região de confiança
• Método Nelder-Mead
• Variantes do algoritmo de Levenberg-Marquardt também foram usadas para resolver sistemas de equações não lineares.

Referências

Leitura adicional

links externos

• A descrição detalhada do algoritmo pode ser encontrada em Receitas Numéricas em C, Capítulo 15.5: Modelos não lineares
• CT Kelley, Iterative Methods for Optimization , SIAM Frontiers in Applied Mathematics, no 18, 1999, ISBN  0-89871-433-8 . Cópia online
• História do algoritmo em notícias SIAM
• Um tutorial de Ananth Ranganathan
• K. Madsen, HB Nielsen, O. Tingleff, Methods for Non-Linear Least Squares Problems (tutorial de mínimos quadrados não lineares; código LM: secante jacobiana analítica )
• T. Strutz: Ajuste de dados e incerteza (uma introdução prática aos mínimos quadrados ponderados e além). 2ª edição, Springer Vieweg, 2016, ISBN  978-3-658-11455-8 .
• HP Gavin, o método Levenberg-Marquardt para problemas de ajuste de curva de mínimos quadrados não lineares ( implementação MATLAB incluída)