\documentclass{TEMA}

\usepackage{algorithm2e}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[brazil]{babel}
\usepackage[latin1]{inputenc}
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{subfigure}

% Formatação de estilos do package algorithm2e
\RestyleAlgo{ruled}
\SetAlgorithmName{Algoritmo}{\autoref}{Lista de Algoritmos}
\SetAlgoCaptionSeparator{:}
\SetAlgoCaptionLayout{vlined}

\newcommand{\up}[1]{\raisebox{1.5ex}[0pt]{#1}}

\begin{document}

\title {Métodos de Regiões de Confiança para Resolução do Problema de
  Quadrados Mínimos: Implementação e Testes Numéricos\thanks{Este
    trabalho contou com os apoios do CNPq (304032/2010-7,
    106731/2011-4 e 151749/2011-6), FAPESP (2012/05725-0) e
    PRONEX-Otimização; Trabalho apresentado no XXXIV CNMAC.}}

\author
    {J.L.C. GARDENGHI%
     \thanks{john@ime.usp.br}\,,
     Departamento de Ciência de Computação, 
     Instituto de Matemática e Estatística, IME,  
     USP - Universidade de São Paulo, 
     05508-090 São Paulo, SP, Brasil.
     \\ \\
     S.A. SANTOS%
     \thanks{sandra@ime.unicamp.br}\,,
     Departamento de Matemática Aplicada,
     Instituto de Matemática, Estatística e Computação Científica, IMECC, 
     UNICAMP - Universidade Estadual de Campinas, 
     13083-859 Campinas, SP, Brasil.}

\criartitulo

\runningheads {Gardenghi e Santos}{Métodos de Regiões de Confiança para Quadrados Mínimos}

\begin{abstract} 

  {\bf 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 \textit{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.

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

\end{abstract}

\newsec{Introdução}

Neste trabalho, consideramos o problema de quadrados mínimos
\begin{equation}\label{eq:prob}
	\min_{x \in \mathbb{R}^n} \frac{1}{2} \| F(x) \|^2 \equiv f(x),
\end{equation}
sendo~$F: \mathbb{R}^n \rightarrow \mathbb{R}^m$ a \textit{função
  residual}, ao menos duas vezes continuamente di\-fe\-ren\-ciá\-vel,
$F(x)$ o \textit{resíduo do problema} em~$x$ e~$\| \cdot \|$ 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 \gg n$. Para um recente panorama do estado da arte deste
problema, tanto do ponto de vista teórico quanto prático,
veja~\cite{yuan} e suas referências.

Na próxima seção, apresentamos dois métodos para resolução
de~(\ref{eq:prob}): 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.

\newsec{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) \mbox{ s.a. } x \in \mathbb{R}^n$, são estratégias iterativas
baseadas em modelos quadráticos para a função objetivo $f$ em torno de
um dado ponto~\cite{cgt}. Em geral, seja~$x_c \in \mathbb{R}^n$ o
ponto corrente, define-se um modelo quadrático~$q_c$ para~$f$ em torno
de~$x_c$ e o subproblema é formulado como
\begin{equation}\label{eq:subprob}
	\min_{p \in \mathbb{R}^n} q_c(p) \mbox{ s.a. } \| p \| \leq \Delta_c,
\end{equation}
em que $0 < \Delta_c \in \mathbb{R}$ é o raio da região de
confiança. A solução~$p_c$, 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~(\ref{eq:subprob}) 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
\begin{equation}\label{eq:apred}
	\rho(p) \stackrel{\mbox{\scriptsize{def}}}{=} \frac{f(x_c) - f(x_c + p)}{q_c(0) - q_c(p)}.
\end{equation}

No Algoritmo~\ref{alg:rc}, página~\pageref{alg:rc}, apresentamos o
algoritmo básico de regiões de confiança. Os parâmetros $\kappa,
\eta_1$ e $\eta_2$ são determinados experimentalmente, embora haja
algumas escolhas clássicas na literatura~\cite{cgt,nocedal}. Sob esta
ótica, escolhemos $\kappa = 10^{-3}$, $\eta_1=\frac{1}{4}$ e
$\eta_2=\frac{3}{4}$. Para resolver o subproblema
(linha~\ref{algrc:subprob}), implementamos LMA e MGC. Uma breve
descrição destes dois métodos é apresentada a seguir.

\begin{algorithm}[Ht]
  \DontPrintSemicolon
  \SetAlgoLined
  \caption{Regiões de Confiança}
  \label{alg:rc}
  \Entrada{Ponto inicial $x_0 \in \mathbb{R}^n$, raio inicial da
região de confiança $\Delta_0$, constantes $\kappa \in
[0,\frac{1}{4}]$, $0<\eta_1\leq\eta_2<1$.}
  \BlankLine
  \nl \Para{$k = 0, 1, \ldots$}
  	  {
	  	\nl \label{algrc:subprob}Encontre $p_k$ como solução (aproximada) do subproblema~(\ref{eq:subprob})\;
	  	\nl $x_{k+1} \leftarrow x_k + p_k$\;
	  	\nl Avalie o quociente de redução relativa $\rho(p_k)$~(\ref{eq:apred})\;
	  	\nl \eSe{$\rho(p_k) > \kappa$}
	  		{
	  			\nl Aceite o ponto candidato $x_{k+1}$ e avalie critérios de parada\;
	  		}
	  		{
	  			\nl Descarte o ponto candidato e considere $x_{k+1} = x_k$\;
	  		}
	  	\BlankLine
	  	\nl \lSe{$\rho(p_k) \leq \eta_1$}{Obtenha $\Delta_{k+1}$ reduzindo $\Delta_k$}\;
	  	\nl \Senao{
		  		\nl \lSe{$\rho(p_k) \geq \eta_2$}{Obtenha $\Delta_{k+1}$ aumentando $\Delta_k$}\;
		  		\nl \lSenao{$\Delta_{k+1} \leftarrow \Delta_k$}}
	  }
\end{algorithm}

\subsection{Método de Levenberg-Marquardt}

A principal característica deste método é atacar a estrutura
particular do problema~(\ref{eq:prob}). Originalmente, esta estratégia
foi proposta por Levenberg~\cite{levenberg} e
Marquardt~\cite{marquardt}, e foi considerada precursora da técnica de
regiões de confiança (ver também~\cite{madsen}). Moré~\cite{more78}
propôs uma série de particularidades para implementação, que foram
exploradas neste trabalho.

Dado um ponto $x_c \in \mathbb{R}^n$, LMA trabalha com a aproximação
\[
q_c(p) \stackrel{\mbox{\scriptsize{def}}}{=} \frac{1}{2} \| J_c p + F_c\|^2
\]
em torno de $x_c$, sendo $F_c \in \mathbb{R}^m$ o resíduo em $x_c$ e
$J_c \in \mathbb{R}^{m \times n}$ a jacobiana da função residual
avaliada em $x_c$. Deste modo, o subproblema de região de
confiança~(\ref{eq:subprob}) é um problema de quadrados mínimos
lineares com passo restrito, cujas condições de otimalidade são dadas
por
\begin{subequations}\label{eq:lma-opt}
	\begin{eqnarray}
		\label{eq:lma-opt1} & (J_c^TJ_c + \lambda^* I)p^* = -J_c^TF_c & \\
		\label{eq:lma-opt2} & \lambda^* (\| p^* \| - \Delta_c) = 0 & \\
		\label{eq:lma-opt3} & \lambda^* \geq 0, \|p^*\| \leq \Delta_c. &
	\end{eqnarray}
\end{subequations}

O sistema linear~(\ref{eq:lma-opt1}) está bem definido somente se
$(J_c^TJ_c + \lambda^* I)$ for definida positiva. Como $J_c^T J_c$ é
ao menos semidefinida positiva, então $(J_c^TJ_c + \lambda^* I)$ será
definida positiva sempre que $\lambda^* > 0$. Em particular, podemos
ter dificuldades numéricas, por mal condicionamento ou por
singularidade, quando $\lambda^* = 0$. Neste caso, de acordo
com~(\ref{eq:lma-opt2}), temos que a solução de~(\ref{eq:subprob})
pode ser interna, ou seja, $\| p^* \| < \Delta_c$. O resultado a
seguir estabelece condições para que a solução de~(\ref{eq:subprob})
seja interior à região de confiança.

\begin{propTEMAp}\label{prop:solfront}
  O problema~(\ref{eq:subprob}) não possui solução global na fronteira
  se, e somente se, $J_c^T J_c$ é definida positiva e o passo de
  Newton é interior, ou seja, $\| (J_c^T J_c)^{-1} J_c^T F_c \| <
  \Delta_c$.
\end{propTEMAp}

\begin{proof}
Ver~\cite[Corolário 1]{rel1}.
\end{proof}

Por outro lado, caso o problema possua solução na fronteira, é sempre
possível encontrar $\lambda^*$ que satisfaça $\| p^* \| =
\Delta_c$. Assim sendo, consideramos
\[
p(\lambda) \stackrel{\mbox{\scriptsize{def}}}{=} - (J_c^T J_c + \lambda I)^\dagger J_c^T F_c,
\]
em que $\dagger$ é a pseudoinversa de Moore-Penrose, e queremos
encontrar $\lambda^* \geq 0$ tal que
\begin{equation}\label{eq:lma-sec}
  \| p(\lambda^*) \| = \Delta_c.
\end{equation}

A equação~(\ref{eq:lma-sec}) é chamada de \textit{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~\ref{sec:imp} Notemos que $\lambda^*$ atua como um
regularizador, ou seja, ele ajusta o sistema~(\ref{eq:lma-opt1}) para
que sempre exista uma solução $p^*$. É justamente essa a conexão de
LMA com regiões de confiança. Originalmente
\cite{levenberg,marquardt}, o parâmetro $\lambda^*$ era computado de
forma heurística. Com regiões de confiança, buscamos o multiplicador
de Lagrange da restrição quadrática $\lambda^* \geq 0$ de tal maneira
que $\| p^* \| = \Delta_c$, salvo o caso particular em que a solução
seja interna e então $\lambda^* = 0$.

\subsection{Método de Gradientes Conjugados}

Neste método, vamos atacar o problema~(\ref{eq:prob}) 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~(\ref{eq:subprob}), tomando
\[
	q_c(p) \stackrel{\mbox{\scriptsize{def}}}{=} \frac{1}{2} p^T H_c p + g_c^T p + f_c
\]
o modelo quadrático de $f$ em torno de um ponto corrente $x_c$, com
$H_c \in \mathbb{R}^{n \times n}$, $g_c \in \mathbb{R}^n$ e $f_c \in
\mathbb{R}$, respectivamente, a matriz hessiana, o vetor gradiente e a
função objetivo $f$, todos avaliados em $x_c$.

Para um estudo mais detalhado sobre a estratégia de gradientes
conjugados, veja~\cite[Seção 3]{rel2}. Para a resolução do
subproblema, a ideia é produzir uma sequência $\{ d_i \}$ de direções
conjugadas e encontrar a solução para~(\ref{eq:subprob})
iterativamente, por meio de $p_{i+1} = p_i + \alpha_i d_i$. MGC seria
uma estratégia natural 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~(\ref{eq:subprob}) são dadas por
\begin{subequations}\label{eq:cg-opt}
	\begin{eqnarray}
		\label{eq:cg-opt1} & (H_c + \lambda ^* I)p^* = -g_c, & \\
		\label{eq:cg-opt2} & \lambda ^*(\| p^* \| - \Delta_c) = 0, & \\
		\label{eq:cg-opt4} & \lambda^* \geq 0, \| p^* \| \leq \Delta_c, & \\
		\label{eq:cg-opt3} & (H_c + \lambda ^* I) \mbox{ é semidefinida positiva}. &
	\end{eqnarray}
\end{subequations}
Toint~\cite{toint} e Steihaug~\cite{steihaug} propuseram algumas
modificações no MGC original de tal maneira que se $\|g_c+H_cp_{i}\|$
for menor que dada tolerância, então temos uma solução interna, pois
$p_i$ satisfaz aproximadamente a condição~(\ref{eq:cg-opt1}) com
$\lambda^*=0$. Por outro lado, se $\|p_{i+1}\|>\Delta_c$, então basta
encontrar $\tau$ de forma que $\|p_{i}+\tau d_{i}\|=\Delta_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 $d_i^T H_c d_i
\leq 0$, computamos $\tau$ tal que $\|p_{i}+\tau d_{i}\|=\Delta_c$.

Deste modo, MGC atua como GC clássico 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~\ref{alg:mgc}
descrevemos MGC.\\

\begin{algorithm}[H]
  \DontPrintSemicolon
  \SetAlgoLined
  \caption{Gradientes Conjugados com Regiões de Confiança}
  \label{alg:mgc}
  \Entrada{Ponto corrente $x_c \in \mathbb{R}^n$, a matriz
    hessiana $H_c \in \mathbb{R}^{n \times n}$ e o gradiente $g_c \in
    \mathbb{R}^n$ da função objetivo avaliados em $x_c$, $\Delta_c$ o raio
    da região de confiança, $\varepsilon > 0$ uma tolerância para o erro
    de otimalidade.}
  \Saida{$p$ direção para computar o novo ponto candidato $x_n = x_c + p$.}
  \BlankLine
  \nl $p_0 \leftarrow 0, r_0 \leftarrow g_c, d_0 \leftarrow r_0$\;
  \nl \Para{$i = 0, 1, \ldots$}
  {
    \nl \Se{$d_i^T H_c d_i \leq 0$}
    {
      \nl Encontre $\tau$ tal que $\| p_i + \tau d_i \| = \Delta_c$\;
      \nl \Retorna $p_i + \tau d_i$\;
    }
    \nl $\displaystyle \alpha_i \leftarrow (r_i^T r_i)/(d_i^T H_c d_i)$\;
    \nl $p_{i+1} \leftarrow p_i + \alpha_i d_i$\;
    \nl \Se{$\| p_{i+1} \| > \Delta_c $}
    {
      \nl Encontre $\tau$ tal que $\| p_i + \tau d_i \| = \Delta_c$\;
      \nl \Retorna $p_i + \tau d_i$\;
    }
    \nl $r_{i+1} \leftarrow r_i + \alpha_i H_c d_i$\;
    \nl \lSe{$\| r_{i+1} \| \leq \varepsilon$}{\Retorna $p_{i+1}$\;}
    \nl $\beta_i \leftarrow (r_{i+1}^T r_{i+1})/(r_i^T r_i)$\;
    \nl $d_{i+1} \leftarrow r_{i+1} + \beta_i d_i$\;
  }
\end{algorithm}

\newsec{Implementação Computacional}\label{sec:imp}

Nesta seção, apresentamos os detalhes de implementação do LMA que, por
explorar as peculiaridades do problema~(\ref{eq:prob}), possui mais
particularidades de implementação que o Algoritmo~\ref{alg:mgc}, que é
um método para resolução do subproblema de regiões de confiança no
caso de minimização irrestrita geral.

\subsection{Álgebra matricial empregada na solução do subproblema}

A equação~(\ref{eq:lma-opt1}) pode ser vista como um sistema de
equações normais para o problema linear de quadrados mínimos:
\begin{equation}\label{eq:linqm}
\left(
\begin{array}{c}
J_c \\ \sqrt{\lambda} I
\end{array}
\right)p \cong
-\left(
\begin{array}{c}
F_c \\ 0
\end{array}
\right).
\end{equation}
Tal equivalência permite resolver o problema sem avaliar o produto
$J_c^T J_c$, reduzindo pela metade a quantidade de operações. Deste
modo, se calcularmos a fatoração~$QR$
\[
J_c = Q \left(
  \begin{array}{c}
    R_c \\ 0
  \end{array}
\right)
\]
teremos
\begin{equation}\label{eq:qr}
  \left(
    \begin{array}{cc}
      Q^T & 0 \\ 0 & I
    \end{array}
  \right) \left(
    \begin{array}{c}
      J_c \\ \sqrt{\lambda} I
    \end{array}
  \right) = \left(
    \begin{array}{c}
      R_c \\ 0 \\ \sqrt{\lambda} I
    \end{array}
  \right).
\end{equation}

\vspace{-0.1cm}

A matriz que resulta de~(\ref{eq:qr}) possui uma estrutura particular:
$R_c$ é uma matriz triangular superior, enquanto $\sqrt{\lambda} I$ é
uma matriz escalar. Diante disto, é possível obter uma matriz
triangular superior equivalente eliminando-se os $n$ termos não-nulos
da matriz $\sqrt{\lambda} I$ com uma sequência de $n(n+1)/2$ rotações
de Givens aplicadas em cada linha $i$ de $\sqrt{\lambda} I$,
simultaneamente a cada linha $i+j$ de $R_c$, para
$j~=~0,~\ldots,~(n-i)$, começando da linha $n$. Expressando tais
rotações em uma matriz $G_\lambda \in \mathbb{R}^{(m+n) \times (m+n)}$
por meio de todos os sucessivos $n(n+1)/2$ produtos, vem que
\begin{equation}\label{eq:giv}
  G_\lambda^T \left(
    \begin{array}{c}
      R_c \\ 0 \\ \sqrt{\lambda} I
    \end{array}
  \right) = \left(
    \begin{array}{c}
      R_\lambda \\ 0 \\ 0
    \end{array}
  \right).
\end{equation}
Disto, obtemos a matriz triangular superior $R_\lambda$, de modo
que~$R_\lambda^T~R_\lambda~=~J_c^T~J_c~+~\lambda~I$. Por conseguinte,
o sistema~(\ref{eq:lma-opt1}) torna-se menos custoso de ser resolvido
para diferentes valores de $\lambda$. 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 $\lambda$, obtém-se $R_\lambda$
por~(\ref{eq:giv}), que consiste em um procedimento mais barato que o
cálculo da fatoração $QR$ propriamente dita.

\subsection{Raiz da equação secular}

O problema~(\ref{eq:lma-sec}) consiste em encontrar $\lambda^* \geq 0$
tal que $\| p(\lambda^*) \| = \Delta_c$. Então, se definimos
$\delta(\lambda) \stackrel{\mbox{\scriptsize{def}}}{=} p(\lambda) -
\Delta_c$, estamos interessados em encontrar uma raiz para
$\delta(\lambda)$. Todavia, uma análise mais
apurada~\cite[Seção~3.2]{rel1} nos mostra que a função $\delta$ pode
possuir polos nos autovalores da matriz $J_c^T J_c + \lambda I$ e pode
não possuir zeros finitos. Em vista disto, uma estratégia mais estável
é resolver a equação secular
\begin{equation}\label{eq:sec}
  \phi(\lambda) \stackrel{\mbox{\scriptsize{def}}}{=} \frac{1}{\| p(\lambda) \|} - \frac{1}{\Delta_c} = 0,
\end{equation}
que, por ser a recíproca de $\delta$, não possui polos finitos e suas
raízes acontecem nos valores negativos dos autovalores da matriz
$J_c^T J_c + \lambda I$. Deste modo, para encontrar $\lambda^* \geq 0$
que satisfaça~(\ref{eq:lma-opt1}), avaliamos primeiro o caso
particular em que $\lambda^* = 0$ (solução
interna). Se~(\ref{eq:lma-opt1}) 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~(\ref{eq:sec}) gera uma
sequência $\{ \lambda_j \}$ tal que
\begin{equation}\label{eq:lambda}
  \lambda_{j+1} = \lambda_j - \frac{\phi(\lambda_j)}{\phi'(\lambda_j)}.
\end{equation}
Em contrapartida, seria interessante reduzir o custo do cálculo de
$\lambda_{j+1}$ e ainda obter uma expressão numericamente mais
estável. Isso é possível, conforme estabelece o resultado a seguir.

\begin{teoTEMA}\label{teo:sec}
  Considere o problema~(\ref{eq:sec}) de encontrar a raiz da equação
  secular. Seja $R_\lambda$ a matriz triangular calculada
  em~(\ref{eq:giv}), tal que $R_\lambda^TR_\lambda p_j = -J_c^TF_c$ e
  \linebreak $R_\lambda^Tq_j = p_j$ e seja $\Delta_c$ o raio da região
  de confiança corrente. A sequência $\{ \lambda_j \}$ gerada pelo
  método de Newton definida em~(\ref{eq:lambda}) é equivalente a:
\begin{equation}\label{eq:nlambda}
  \lambda_{j+1} = \lambda_{j} + \left( \frac{\| p_j \|}{\| q_j \|} \right)^2 \left( \frac{\| p_j \| - \Delta_c}{\Delta_c} \right).
\end{equation}
\end{teoTEMA}

\begin{proof}
  Ver~\cite[Lema~1]{rel1}.
\end{proof}

Por fim, para que o método de Newton não convirja para uma raiz
indesejada de $\phi(\lambda)$ (pois estamos interessados em raízes
positivas), usamos uma salvaguarda de modo que, a cada iteração,
impomos $\lambda_j \in [ \ell_j, u_j ]$, sendo $\ell_0 = 0$ a
estimativa para o menor autovalor de $J_c^T J_c + \lambda I$ e $u_0$
obtido pelo quociente de Rayleigh
\[
(\mu_1 + \lambda)^2 \leq \frac{p^T(J_c^TJ_c + \lambda I)^T(J_c^TJ_c +
  \lambda I)p}{p^Tp} \leq (\mu_n + \lambda)^2,
\]
sendo $\mu_1, \mu_n$ o menor e o maior autovalores de $J_c^T J_c$,
respectivamente. Como zero é um limitante inferior para $\mu_1$, e
estamos procurando $\| p \| = \Delta_c$, vem que
\[
\lambda \leq \frac{\| J_c^TF_c \|}{\Delta_c} = u_0.
\]
O intervalo da salvaguarda é atualizado a cada iteração da seguinte
maneira:
\[
\mbox{se $\delta(\lambda_j) \geq 0$, então } \ell_{j+1} = \lambda_j +
\frac{\| q_j \|^2}{\| p_j \|}, \mbox{ caso contrário } u_{j+1} =
\lambda_j.
\]
Para maiores detalhes, ver~\cite[Seção~4.2.1]{rel1}.

\subsection{Quociente de redução relativa}

Nesta subseção, nosso objetivo é eliminar a ocorrência de
\textit{underflows} e \textit{overflows} no cálculo
de~(\ref{eq:apred}). Usando LMA, podemos calcular o quociente como
\begin{equation}\label{eq:lma-apred}
  \displaystyle
  \rho(p) = \frac{\| F_c \|^2 - \| F_n \|^2}{\| F_c \|^2 - \| F_c + J_cp \|^2},
\end{equation}
em que $F_n = F(x_c + p)$. De~(\ref{eq:lma-opt1}), temos que
\begin{subequations}\label{eq:rels}
  \begin{eqnarray}
    \label{eq:rel1} J_c^TJ_cp & = & -J_c^TF_c - \lambda p \\
    \label{eq:rel2} -F_c^TJ_c & = & \lambda p^T + p^TJ_c^TJ_c.
  \end{eqnarray}
\end{subequations}
Algumas manipulações com~(\ref{eq:rels}) em~(\ref{eq:lma-apred})
implicam que
\[
\| F_c \|^2 - \| F_c + J_cp \|^2 = \| J_cp \|^2 + 2\lambda\| p \|^2.
\]
Portanto, podemos reescrever o quociente~(\ref{eq:lma-apred})
dividindo todos os termos por $\| F_c \|^2$ e obtendo a seguinte
relação:
\begin{equation}\label{eq:lma-qred2}
  \rho(p) = \frac{1 - \left( \displaystyle \frac{\| F_n \|}{\| F_c \|} \right)^2}{\left( \displaystyle \frac{\| J_cp \|}{\| F_c \|} \right)^2 + 2\lambda \left( \displaystyle \frac{\| p \|}{\| F_c \|} \right)^2}.
\end{equation}

Deste modo, o denominador de~(\ref{eq:lma-qred2}) será sempre
não-negativo e não irá gerar \textit{overflows}. Por outro lado, o
numerador de~(\ref{eq:lma-qred2}) jamais gerará \textit{overflows} a
não ser que $\| F_n \| \gg \| F_c \|$. Neste caso, temos que o modelo
não foi reduzido com relação à função objetivo, portanto não é
necessário calcular $\rho(p)$.

\subsection{Atualização do raio da região de confiança}

Utilizamos a aproximação inicial $\Delta_0 = \| J_0^T F_0 \|/10$,
clássica da literatura~\cite{cgt}. 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~\ref{alg:rc}. Para
aumentar o tamanho do raio, fazemos $\Delta_n = 2 \| p^* \|$. Para
reduzir o tamanho do raio, fazemos $\Delta_n=\vartheta^* \Delta_c$. A
fim de definir um valor para $\vartheta^*$, seguimos
Moré~\cite{more78}: definimos \linebreak $\gamma(\vartheta)
\stackrel{\mbox{\scriptsize{def}}}{=} 1/2 \| F(x_c + \vartheta p^*)
\|^2$ e encontramos o minimizador $\vartheta^*$ da parábola
interpolante
\[
m_q(\vartheta) \stackrel{\mbox{\scriptsize{def}}}{=} [\gamma(1) -
\gamma(0) - \gamma'(0)]\vartheta^2 + \gamma'(0)\vartheta + \gamma(0).
\]
Como $m_q$ é suave e convexa, então temos que $m_q'(\vartheta^*) = 0$,
ou seja,
\begin{equation}\label{eq:vartheta}
  \vartheta^* = \displaystyle -\frac{\gamma'(0)}{2[\gamma(1) - \gamma(0) - \gamma'(0)]}
             = \displaystyle -\frac{(p^*)^TJ_c^TF_c}{\| F_n \|^2 - \| F_c \|^2 - 2(p^*)^TJ_c^TF_c}.
\end{equation}
Pré-multiplicando~(\ref{eq:lma-opt1}) avaliado em $p^*$ por $(p^*)^T$,
temos
\[
-(p^*)^T(J_c^TJ_c + \lambda^*I)p^* = (p^*)^TJ_c^TF_c,
\]
donde segue que
\begin{equation}\label{eq:reltheta}
  \frac{(p^*)^TJ_c^TF_c}{\|F_c\|^2} = -\left[\left(\frac{\|J_cp^*\|}{\|F_c\|}\right)^2 + \lambda^*\left(\frac{\|p^*\|}{\|F_c\|}\right)^2\right].
\end{equation}
Dividindo todos os termos de~(\ref{eq:vartheta}) por $\| F_c \|^2$,
multiplicando o numerador e o denominador por $-1/2$ e
substituindo~(\ref{eq:reltheta}), que denotamos por $\alpha$, obtemos
\begin{equation}\label{eq:nvartheta}
  \vartheta^* = \frac{\displaystyle \frac{1}{2}\alpha}{\displaystyle \alpha + \frac{1}{2}\left[1 - \left(\frac{\| F_n \|}{\| F_c \|}\right)^2\right]}.
\end{equation}
Pelos mesmos motivos expostos na seção anterior, (\ref{eq:reltheta}) é
uma expressão numericamente estável, e~(\ref{eq:nvartheta}) é
confiável e barata para encontrar $\vartheta^*$ e definir $\Delta_n$.

\newsec{Experimentos Numéricos}

Ambos os métodos foram implementados no \textit{CAS
  Maxima}~\cite{rpmmax,maxima}. O \textit{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 \textit{Maxima} por ser um sistema gratuito e
livre, além de não possuir ferramentas de otimização implementadas,
como é o caso do Matlab\textsuperscript{\textregistered}. 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~\cite{mgh} e estão numerados de acordo com a referência
original, bem como as dimensões e os pontos iniciais adotados (ver
detalhes em~\cite{rel1, rel2}).

Como critérios de parada, adotamos a satisfação relativa das condições
de primeira ordem, isto é, $x_k$ é solução em LMA se $\left| x_k^T
  J_k^T F_k/\max\{ \| F_k \|, 1 \} \right| \leq 10^{-6}$ e em MGC se
$\| g_k \| / \| g_0 \| \leq 10^{-6}$. Além disso, em LMA impusemos
parada se o resíduo fosse pequeno, ou seja, $\| F_k \| \leq
10^{-6}$. 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~\ref{tab:res},
enquanto os perfis de desempenho~\cite{pp} na
Figura~\ref{fig:pp}. Neles, consideramos \texttt{VF} a função objetivo
avaliada na solução, \texttt{AH} a quantidade de avaliações da matriz
hessiana (para o MGC), \texttt{AJ} a quantidade de avaliações da
matriz Jacobiana (para o LMA), \texttt{ITE} a quantidade de iterações
externas, \texttt{ITI} a quantidade de iterações internas (isto é,
para resolução dos subproblemas de regiões de confiança) e \texttt{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 \texttt{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.

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.

\begin{table}[!htb]
  \centering
  {\footnotesize
    \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
      \hline
      &          \multicolumn{5}{c|}{MGC}    &         \multicolumn{5}{c|}{LMA}     \\ \cline{2-11}
      \up{Função} &       VF       & AH & ITE & ITI & EI &       VF       & AJ & ITE & ITI & EI \\
      \hline \hline
           1$^*$  & $0.34654$ E-08 & 22 & 22  &  43 & 2.0&     $0.0$    & 14 &  18 &  51 & 3.6  \\ \hline
           2      & $0.69988$ E 01 & 8  &  7  &  14 & 1.8& $0.69988$ E 01 & 28 &  40 &  95 & 3.4  \\ \hline
           6      & $0.11151$ E 02 & 9  &  8  &  16 & 1.8& $0.11151$ E 02 & 9  &  19 &  49 & 5.4  \\ \hline
           8      & $0.90635$ E-01 & 14 & 14  &  42 & 3.0& $0.90635$ E-01 & 6  &   5 &   5 & 0.8  \\ \hline
           10     & $0.93779$ E 01 & 199& 216 & 826 & 4.2& $0.93779$ E 01 & 124& 143 & 299 & 2.4  \\ \hline
           12$^*$ & $0.92759$ E-04 & 13 & 13  &  34 & 2.6& $0.62754$ E-10 & 12 &  12 &  34 & 2.8  \\ \hline
           13     & $0.75374$ E-03 & 13 & 12  &  48 & 3.7& $0.19361$ E-03 & 9  &   8 &   8 & 0.9  \\ \hline
           15     & $0.17535$ E-01 & 11 & 11  &  35 & 3.2& $0.17536$ E-01 & 11 &  12 &  43 & 3.9  \\ \hline
           16     & $0.29295$ E 03 & 8  & 7   &  33 & 4.1& $0.29295$ E 03 & 24 &  40 & 127 & 5.3  \\ \hline
           17     & $0.73924$ E-02 & 34 & 40  & 201 & 5.9& $0.73924$ E-02 & 16 &  18 &  35 & 2.2  \\ \hline
           19     & $0.20034$ E 00 & 19 & 22  & 179 & 9.4& $0.20034$ E 00 & 14 &  14 &  38 & 2.7  \\ \hline
           20     & $0.22180$ E-03 & 13 & 12  & 149 &11.5& $0.38213$ E-04 & 5  &   4 &   6 & 1.2  \\ \hline
           27$^*$ & $0.89312$ E-05 & 7  & 6   &  13 & 1.9& $0.19146$ E-09 & 11 &  14 &  23 & 2.1  \\ \hline
           32     & $0.67082$ E 01 & 4  & 3   &   3 & 0.8& $0.67082$ E 01 & 5  &   4 &   7 & 1.4  \\ \hline
           33     & $0.34826$ E 01 & 2  & 2   &   1 & 0.5& $0.34826$ E 01 & 3  &   2 &   6 & 2.0  \\ \hline
           34     & $0.36917$ E 01 & 2  & 2   &   1 & 0.5& $0.36917$ E 01 & 2  &   1 &   2 & 1.0  \\ \hline \hline
         Média    &        --      &23.6& 24.8&102.4& 3.5&       --       &18.3& 22.1& 51.8& 2.6 \\ \hline
    \end{tabular}
  }
  \caption{Quadro comparativo entre os resultados das rotinas LMA e MGC.}
  \label{tab:res}
\end{table}

\begin{figure}[tb]
  \centering
  \subfigure[Avaliações de Função.]{\label{fig:af}\includegraphics[width=6.2cm,keepaspectratio]{perf-af.eps}} \
  \subfigure[Avaliações Matriciais.]{\label{fig:ag}\includegraphics[width=6.2cm,keepaspectratio]{perf-ag.eps}}\\
  \subfigure[Iterações Externas.]{\label{fig:ITE}\includegraphics[width=6.2cm,keepaspectratio]{perf-itext.eps}} \
  \subfigure[Iterações Internas.]{\label{fig:itint}\includegraphics[width=6.2cm,keepaspectratio]{perf-itint.eps}}%
  \caption{Perfis de desempenho dos algoritmos LMA e MGC.}
  \label{fig:pp}
\end{figure}

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 sequê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.

\newsec{Conclusão}

O algoritmo LMA é uma proposta atraente para a resolução
de~(\ref{eq:prob}), 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 $\lambda I$. Em
contrapartida, o método perde desempenho em problemas extremamente não
lineares ou se o resíduo for muito grande na solução $x^*$.

Em comparação com MGC, LMA aproveita a estrutura do
problema~(\ref{eq:lma-opt1}) 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.

\begin{abstract} 

  {\bf 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.

  {\bf Keywords}. Least squares, trust region, computational
  implementation.

\end{abstract}

\begin{thebibliography}{99}

  \vspace{-1mm}

\bibitem{rpmmax} L.N. Andrade, Maxima: um completo programa de
  computação algébrica, {\em Revista do Professor de Matemática},
  n.~77, 2012.

\bibitem{cgt} A.R. Conn, N.I.L. Gould, Ph.L. Toint,
  ``Trust-region Methods''. SIAM, Philadelphia, 2000.

\bibitem{pp} E.D. Dolan, J.J. Moré, Benchmarking optimization software
  with performance profiles, {\em Mathematical Programming}, v.~91,
  pp.~201-213, 2002.
		
\bibitem{rel1} 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.

\bibitem{rel2} 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.
	
\bibitem{levenberg} K. Levenberg, A method for the solution of certain
  non-linear problems in least squares, {\em The Quarterly of Applied
    Mathematics 2}, pp. 164-168, 1944.

\bibitem{madsen} K. Madsen, An algorithm for the minimax solution of
  overdetermined systems of nonlinear equations. {\em Journal of the
    Institute of Mathematics and its Applications}. v.~16(3),
  pp.~321-328, 1975.

\bibitem{marquardt} D.W. Marquardt, An algorithm for least-squares
  estimation of nonlinear parameters, {\em SIAM Journal on Applied
    Mathematics 11}, pp. 431-441, 1963.
	
\bibitem{more78} J.J. Moré, The Levenberg-Marquardt algorithm:
  implementation and theory. {\em Lecture Notes in Mathematics 630:
    Numerical Analysis}, Springer-Verlag, New York, pp. 105-116, 1978.
	
\bibitem{mgh} J.J. Moré, B.S. Garbow, K.E. Hillstrom, Testing
  unconstrained optimization software, {\em ACM Transactions on
    Mathematical Software}, v.~7, pp. 17-41, 1981.
	
\bibitem{nocedal} J. Nocedal, S.J. Wright, ``Numerical
  Optimization''. Springer, New York, 1999.
	
\bibitem{steihaug} T. Steihaug, The conjugate gradient method and
  trust region in large scale optimization. {\em SIAM Journal on
    Numerical Analysis}, v.~20(3), pp.~626-637,~1983.
	
\bibitem{toint} 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.
	
\bibitem{yuan} Y. Yuan, Recent advances in numerical methods for
  nonlinear equations and nonlinear least-squares. {\em Numerical
    Algebra, Control and Optimization}, v.~1, n.~1, pp.~15-34, 2011.

\bibitem{maxima} http://maxima.sourceforge.net.
	
\end{thebibliography}

\end{document}