Algoritmo de Levenberg–Marquardt

Visão geral

O algoritmo de Levenberg–Marquardt (Levenberg–Marquardt (LM) algorithm) é um método clássico para resolver problemas de mínimos quadrados não lineares (nonlinear least-squares). Ele é amplamente usado em computação científica (scientific computing) e em áreas próximas à IA (AI), como visão computacional (computer vision) (por exemplo, ajuste de feixe (bundle adjustment)), porque combina:

• a convergência local rápida do Método de Gauss–Newton (Gauss–Newton Method) quando o modelo é bem comportado, e
• a robustez da Descida do Gradiente (Gradient Descent) quando o iterando atual está longe da solução.

Conceitualmente, o LM pode ser entendido como um método de Gauss–Newton amortecido (damped Gauss–Newton) ou, de forma equivalente, um método de região de confiança (trust-region method) adaptado a mínimos quadrados.

Formulação do problema: mínimos quadrados não lineares

O LM visa objetivos da forma:

[ \min_{\theta \in \mathbb{R}^p} ; F(\theta) = \frac{1}{2}\sum_{i=1}^m r_i(\theta)^2 ;=; \frac{1}{2}|r(\theta)|_2^2 ]

• (\theta) são os parâmetros a otimizar (por exemplo, parâmetros de curva, poses de câmera).
• (r(\theta)\in\mathbb{R}^m) é o vetor de resíduos.
• (J(\theta) = \frac{\partial r}{\partial \theta}\in\mathbb{R}^{m\times p}) é a Jacobiana (Jacobian) dos resíduos.

O gradiente e a Hessiana (Hessian) (aproximada) têm uma estrutura especial:

[ \nabla F(\theta) = J(\theta)^\top r(\theta) ]

[ \nabla^2 F(\theta) = J^\top J + \sum_{i=1}^m r_i(\theta)\nabla^2 r_i(\theta) ]

O método de Gauss–Newton descarta o segundo termo (muitas vezes pequeno próximo de uma boa solução), resultando na aproximação:

[ \nabla^2 F(\theta) \approx J^\top J ]

O LM modifica o Gauss–Newton para ser mais estável e ter melhor convergência global.

De Gauss–Newton para LM: a equação de atualização

Passo de Gauss–Newton

O Gauss–Newton resolve o sistema linear:

[ (J^\top J),\Delta = -J^\top r ]

e atualiza:

[ \theta \leftarrow \theta + \Delta ]

Isso pode ser muito rápido perto do ótimo, mas pode falhar quando (J^\top J) é mal condicionado, singular, ou quando a linearização é ruim.

Passo de Levenberg–Marquardt (Gauss–Newton amortecido)

O LM introduz um parâmetro de amortecimento (damping) (\lambda \ge 0) e resolve:

[ (J^\top J + \lambda I),\Delta = -J^\top r ]

Então:

[ \theta \leftarrow \theta + \Delta ]

Isso se parece com uma equação normal regularizada por ridge (ridge-regularized). O efeito de (\lambda):

• (\lambda) pequeno: o passo se aproxima de Gauss–Newton.
• (\lambda) grande: o passo se aproxima da descida do gradiente (escalada).

Para ver o comportamento de descida do gradiente, observe que quando (\lambda) domina:

[ (J^\top J + \lambda I)^{-1} \approx \frac{1}{\lambda}I \Rightarrow \Delta \approx -\frac{1}{\lambda}J^\top r = -\frac{1}{\lambda}\nabla F(\theta) ]

Assim, o LM interpola entre comportamento de segunda ordem e de primeira ordem.

Variante de escalonamento diagonal de Marquardt

Uma variante prática comum substitui (\lambda I) por (\lambda D), onde (D) é tipicamente (\operatorname{diag}(J^\top J)):

[ (J^\top J + \lambda,\operatorname{diag}(J^\top J)),\Delta = -J^\top r ]

Isso torna o amortecimento sensível à escala (parâmetros com diferentes unidades/magnitudes se comportam melhor).

Interpretação por região de confiança (por que o LM é robusto)

O LM pode ser derivado de um subproblema de região de confiança aplicado ao modelo quadrático local.

Linearize os resíduos em torno de (\theta):

[ r(\theta+\Delta) \approx r + J\Delta ]

Em seguida, minimize o modelo local de mínimos quadrados:

[ \min_\Delta \frac{1}{2}|r + J\Delta|^2 \quad \text{s.t.} \quad |\Delta| \le \delta ]

A restrição de região de confiança impede dar um passo em que a aproximação linear não seja confiável. Nessa visão:

• (\delta) (raio da região de confiança) corresponde indiretamente a (\lambda).
• (\lambda) maior impõe passos menores (mais conservador).
• (\lambda) menor permite passos maiores, mais parecidos com Gauss–Newton.

Essa perspectiva motiva o amortecimento adaptativo com base em quão bem o modelo prevê a melhoria real.

Escolhendo e atualizando o parâmetro de amortecimento

O sucesso prático do LM frequentemente depende de como você atualiza (\lambda). A estratégia típica usa a razão de ganho (gain ratio):

• Redução prevista pelo modelo quadrático: [ \text{pred} = \frac{1}{2}\left(|r|^2 - |r + J\Delta|^2\right) ]
• Redução real: [ \text{act} = \frac{1}{2}\left(|r|^2 - |r(\theta+\Delta)|^2\right) ]
• Razão de ganho: [ \rho = \frac{\text{act}}{\text{pred}} ]

Interpretação:

• Se (\rho \approx 1), o modelo linearizado previu bem a melhoria → confie mais nele → diminua (\lambda).
• Se (\rho \le 0), o passo não ajudou → aumente (\lambda) e tente um passo mais conservador.

Uma heurística comum de atualização

Uma regra amplamente usada (no estilo dos métodos de LM / região de confiança de Moré):

• Se (\rho > 0): aceite o passo (\theta \leftarrow \theta+\Delta)
  • diminua (\lambda) (por exemplo, (\lambda \leftarrow \lambda \cdot \max(\tfrac{1}{3}, 1-(2\rho-1)^3)))
• Se (\rho \le 0): rejeite o passo
  • aumente (\lambda) (por exemplo, (\lambda \leftarrow 2\lambda))

Muitas implementações usam regras multiplicativas mais simples:

• em sucesso: (\lambda \leftarrow \lambda / \nu)
• em falha: (\lambda \leftarrow \lambda \cdot \nu) com (\nu \in [2, 10]).

Inicialização de \(\lambda\)

Uma escolha comum é:

[ \lambda_0 = \tau \cdot \max(\operatorname{diag}(J^\top J)) ]

com (\tau) pequeno (por exemplo, (10^{-3})). Isso torna (\lambda) aproximadamente consistente com a escala de curvatura.

Algoritmo (pseudocódigo prático)

given θ0, λ0 > 0
repeat until convergence:
    r  ← r(θ)
    J  ← J(θ)
    solve (JᵀJ + λI) Δ = -Jᵀr

    if Δ is tiny: stop

    compute actual reduction:
        act = 0.5*(||r||² - ||r(θ+Δ)||²)

    compute predicted reduction using linear model:
        pred = 0.5*(||r||² - ||r + JΔ||²)

    ρ = act / pred

    if ρ > 0:
        θ ← θ + Δ          (accept)
        decrease λ
    else:
        increase λ         (reject and retry)

Os critérios de parada normalmente combinam:

• gradiente pequeno (|J^\top r|),
• passo pequeno (|\Delta|),
• pequena redução no objetivo,
• máximo de iterações.

Intuição de convergência e garantias (na prática)

O LM é popular porque se comporta bem em diferentes regimes:

• Longe do ótimo: (\lambda) grande faz com que ele se comporte como descida do gradiente, que é mais lenta, mas estável.
• Perto de uma boa solução: (\lambda \to 0) e o LM vira Gauss–Newton, frequentemente fornecendo convergência local rápida (muitas vezes quase quadrática) sob condições regulares padrão (Jacobiana de posto completo, resíduos pequenos etc.).

Ressalvas importantes:

• O LM não tem garantia de encontrar um mínimo global para problemas não convexos (muitos problemas de mínimos quadrados são não convexos).
• Como outros métodos locais, ele é sensível à inicialização, embora tipicamente menos frágil do que o Gauss–Newton puro.

Para um contexto mais amplo sobre por que esses problemas podem ser difíceis, veja Otimização Convexa (Convex Optimization) (para o caso “agradável”) e Otimização Estocástica (Stochastic Optimization) (para a prática de AM em larga escala).

Estabilidade numérica e detalhes de implementação

Evite formar equações normais quando possível

Resolver ((J^\top J + \lambda I)\Delta = -J^\top r) formando explicitamente (J^\top J) pode amplificar problemas de condicionamento (as equações normais elevam ao quadrado o número de condição (condition number)).

Abordagens mais estáveis:

• fatoração QR (QR factorization) de (J) (ou de um sistema aumentado)
• decomposição de Cholesky (Cholesky) em (J^\top J + \lambda I) (frequentemente OK porque o amortecimento o torna SPD)
• decomposição em valores singulares (SVD) (mais robusta, mais cara)

Um truque estável comum é resolver o sistema de mínimos quadrados aumentado:

[ \min_\Delta \left|\begin{bmatrix}J \ \sqrt{\lambda}I\end{bmatrix}\Delta + \begin{bmatrix}r \ 0\end{bmatrix}\right|_2^2 ]

Isso evita formar explicitamente (J^\top J) e se mapeia diretamente para solvers baseados em QR.

Escalonamento importa

O LM pode ter desempenho ruim se os parâmetros estiverem em escalas muito diferentes (por exemplo, metros vs radianos). Correções comuns:

• usar o escalonamento diagonal de Marquardt (\lambda \operatorname{diag}(J^\top J)),
• reescalar parâmetros para magnitudes similares,
• normalizar entradas, ou usar parametrizações bem escolhidas (especialmente em geometria).

Esparsidade para problemas grandes

Em aplicações como ajuste de feixe, (J) é grande, mas esparso (sparse). Implementações eficientes de LM exploram esparsidade usando:

• Cholesky esparso,
• truques de complemento de Schur (Schur complement) (por exemplo, eliminar primeiro pontos 3D, resolver para câmeras),
• solvers iterativos com bons pré-condicionadores (preconditioners).

Exemplo prático 1: ajuste de curvas com LM

Suponha que você queira ajustar um decaimento exponencial:

[ y \approx a e^{-bx} + c ]

Parâmetros: (\theta = [a, b, c]). Resíduos para dados ((x_i, y_i)):

[ r_i(\theta) = a e^{-b x_i} + c - y_i ]

Entradas da Jacobiana:

• (\frac{\partial r_i}{\partial a} = e^{-b x_i})
• (\frac{\partial r_i}{\partial b} = -a x_i e^{-b x_i})
• (\frac{\partial r_i}{\partial c} = 1)

Uma iteração mínima de LM (caso denso) fica assim:

import numpy as np

def lm_step(x, y, theta, lam):
    a, b, c = theta
    exp = np.exp(-b * x)
    r = a * exp + c - y                 # residuals (m,)

    J = np.column_stack([
        exp,                            # dr/da
        -a * x * exp,                   # dr/db
        np.ones_like(x)                 # dr/dc
    ])                                  # (m,3)

    A = J.T @ J + lam * np.eye(3)
    g = J.T @ r
    delta = np.linalg.solve(A, -g)
    return theta + delta, r, J, delta

Em uma implementação completa, você:

• avaliaria se o novo (\theta) reduz (|r|^2),
• adaptaria lam com base em (\rho),
• pararia com base em (|g|), (|\Delta|) etc.

Isso ilustra bem o ponto forte do LM: contagem de parâmetros pequena a média, resíduos diferenciáveis e desejo de convergência rápida.

Exemplo prático 2: ajuste de feixe (visão computacional)

O ajuste de feixe refina parâmetros de câmera e posições de pontos 3D para minimizar o erro de reprojeção:

• parâmetros (\theta): intrínsecos/extrínsecos de câmera + pontos 3D
• resíduos: erros de reprojeção em pixel 2D para pontos observados

O objetivo é uma soma de resíduos ao quadrado, muitas vezes com perdas robustas (ver abaixo). O LM é um solver padrão porque:

• a Jacobiana é esparsa e estruturada,
• passos ao estilo Gauss–Newton convergem rapidamente perto da solução,
• o amortecimento melhora a robustez quando a inicialização não é perfeita.

A maioria dos sistemas modernos de BA implementa LM esparso com:

• eliminação via complemento de Schur,
• parametrização cuidadosa de rotações (por exemplo, atualizações via álgebra de Lie (Lie algebra updates)),
• atualizações de amortecimento de região de confiança.

Perdas robustas e “LM no mundo real”

Dados reais frequentemente contêm outliers, então praticantes muitas vezes substituem mínimos quadrados puros por objetivos robustos (robust):

[ \min_\theta \sum_i \rho\left(r_i(\theta)\right) ]

Perdas robustas comuns: Huber, Cauchy, Tukey. Uma abordagem padrão é mínimos quadrados reponderados iterativamente (iteratively reweighted least squares, IRLS), em que cada iteração resolve um problema de mínimos quadrados ponderados:

[ \min_\theta \frac{1}{2}\sum_i w_i r_i(\theta)^2 ]

Então o LM é aplicado aos resíduos/Jacobiana ponderados. Isso é comum em visão e robótica para reduzir sensibilidade a correspondências ruins.

Relação com o treinamento de redes neurais

O LM já foi usado para treinar redes neurais pequenas e médias (especialmente na literatura mais antiga) porque pode convergir em menos iterações do que métodos de primeira ordem. Porém, ele tipicamente é impraticável para redes modernas grandes porque:

• requer Jacobianas (ou produtos Jacobiana-vetor) e resolve sistemas lineares no espaço de parâmetros,
• memória e tempo escalam mal com a quantidade de parâmetros.

Para aprendizado em larga escala, variantes de Descida do Gradiente (SGD, Adam) dominam, enquanto o LM permanece importante em problemas estruturados de mínimos quadrados (visão, geometria, identificação de sistemas (system identification)).

Quando usar LM (e quando não usar)

LM é uma ótima escolha quando

• seu objetivo é naturalmente uma soma de resíduos ao quadrado,
• a dimensão dos parâmetros é moderada (ou a Jacobiana é esparsa e explorável),
• você consegue calcular Jacobianas com precisão (diferenciação analítica ou automática),
• você quer convergência rápida perto de uma boa solução.

LM é menos indicado quando

• o problema é extremamente grande e denso (milhões de parâmetros),
• resíduos não são suaves ou são caros de avaliar,
• você só tem gradientes estocásticos (stochastic gradients) ruidosos (típico de treinamento de deep learning),
• restrições são centrais (o LM é sem restrições por padrão; veja Otimização com Restrições (Constrained Optimization) para ferramentas dedicadas).

Armadilhas comuns e boas práticas

• Escalonamento ruim: reescale parâmetros ou use amortecimento diagonal.
• Jacobianas por diferenças finitas: podem ser ruidosas; prefira diferenciação analítica/automática.
• Equações normais: evite formar (J^\top J) quando a estabilidade for importante; use QR/SVD ou sistemas aumentados.
• Parar cedo demais: verifique tanto a norma do gradiente quanto o tamanho do passo; não apenas a queda da perda.
• Não convexidade (nonconvexity): use boa inicialização (por exemplo, a partir de um método mais simples) e considere multi-start quando viável.

Resumo

O algoritmo de Levenberg–Marquardt é um “pau para toda obra” para mínimos quadrados não lineares:

• Ele atualiza parâmetros resolvendo ((J^\top J + \lambda I)\Delta = -J^\top r).
• O amortecimento (\lambda) permite ao LM interpolar entre Gauss–Newton e descida do gradiente.
• Uma visão de região de confiança explica sua robustez: (\lambda) grande significa passos pequenos e seguros; (\lambda) pequeno significa passos rápidos ao estilo Gauss–Newton.
• O desempenho prático depende fortemente de regras de atualização do amortecimento, escalonamento e álgebra linear numericamente estável.
• Ele é fundamental em aplicações como ajuste de curvas e ajuste de feixe, onde objetivos estruturados de mínimos quadrados são comuns.

Para um contexto mais profundo, o LM se posiciona naturalmente entre o Método de Gauss–Newton e a Descida do Gradiente, combinando os pontos fortes de ambos.

---