\documentclass[12pt]{article}
\usepackage{graphicx,amsmath,amsfonts,amssymb,epsfig,euscript}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}

\newtheorem{ejer}{Ejercicio}
\newcommand{\bej}{\begin{ejer}\rm}
\newcommand{\fej}{\end{ejer}}

\topmargin-2cm \vsize 29.5cm \hsize 21cm
\setlength{\textwidth}{16.75cm}\setlength{\textheight}{23.5cm}
\setlength{\oddsidemargin}{0.0cm}
\setlength{\evensidemargin}{0.0cm}

\begin{document}
 \centerline{{\small Universidad de Buenos Aires - Facultad de Ciencias Exactas y Naturales - Depto. de Matemática}}
 
 \vskip 0.2cm
 \hrulefill
 \vskip 0.2cm

 \centerline{{\bf\Huge {\sc Optimización}}}
 \vskip 0.2cm
 \centerline{\ttfamily Primer Cuatrimestre 2015}
 \hrulefill

 \bigskip
 \centerline{\bf Trabajo Práctico N$^\circ$ 1: Ajuste de círculos.}
 \bigskip

 
 El objetivo de este trabajo es implementar y comparar algoritmos que ajusten conjuntos de datos con círculos, siguiendo diferentes metodologías. 
 A lo largo de todo el texto, asumiremos que los datos vienen dados por vectores ${\bf x} = (x_1,\dots,x_n)^t$ e ${\bf y} = (y_1,\dots,y_n)^t$, de manera tal que los puntos $(x_i,y_i)$ caigan aproximadamente en un círculo. 
 
\subsection*{Ajuste algebraico}
Un círculo puede definirse implíctamente con la ecuación:
$$\alpha (x^2 + y^2) + \beta x + \gamma y + \delta = 0.$$
Para garantizar unicidad en la escritura, asumiremos $\alpha = 1$. Además, según veremos más adelante, resultará útil reescribir los parámetros, notando:
$$ x^2+y^2 -2ax-2by -c = 0.$$

Si todos los datos cayeran sobre un círculo tendríamos que la ecuación se satisface para todo $(x_i,y_i)$. Asumiendo que esto no sucede (los datos contienen error), buscaremos valores de $a$, $b$ y $c$ que minimicen:
$$\sum_i (x_i^2 + y_i^2 -2ax_i-2by_i-c)^2.$$
O, equivalentemente:
\begin{equation}\label{algebraico}
\min_{a,b,c} \big\|2a{\bf x} + 2b{\bf y} + c - {\bf p}\big\|, \quad\quad\textrm{donde }\; {\bf p} =\big(x_1^2+y_1^2, \dots,x_n^2+y_n^2\big)^t. 
\end{equation}

Este enfoque se conoce como \emph{ajuste algebraico} del círculo, dado que se propone ajustar los parámetros de la ecuación algebraica que lo define. 

\bej Observar que \eqref{algebraico} es un problema de cuadrados mínimos y escribir la matriz del ajuste. Implementar un algoritmo que reciba como datos dos vectores ${\bf x}$ e ${\bf y}$ y resuelva \eqref{algebraico} utilizando la descomposición $QR$ implementada en el Ejercicio 18 b) de la Práctica 1.\fej 




\subsection*{Ajuste geométrico}

El ajuste algebraico, al ceñirse a un modelo de cuadrados mínimos clásico, puede resolverse de manera simple y veloz. Sin embargo, no tiene en cuenta la naturaleza geométrica del problema: no minimiza ninguna distancia. Esto hace que en algunos casos se obtengan resultados muy distintos de los esperados. El \emph{ajuste geométrico} consiste en minimizar las distancias de los datos al círculo. 

 Notamos $C(a,b,r)$ al círculo con centro $(a,b)$ y radio $r$ y 
 $$d_i = \sqrt{(x_i-a)^2+(y_i-b)^2}-r,$$
 la distancia de $(x_i,y_i)$ a $C(a,b,r)$. El problema de determinar los parámetros $a,b,r$ que mejor ajusten los datos ${\bf x}$,${\bf y}$ está dado por:
 \begin{equation}\label{geometrico}
 \min F(a,b,r)
 \end{equation}
donde 
$$ F(a,b,r) = \sum_{i=1}^n d_i^2 = \sum_{i=1}^n\Big(\sqrt{(x_i-a)^2+(y_i-b)^2}-r\Big)^2.$$
Es importante observar que este es un problema de cuadrados mínimos \emph{no lineal}, por lo cual no es posible aplicar de manera directa las técnicas habituales de cuadrados mínimos. 

En una primera aproximación podríamos modificar el funcional de manera que resulte posible aplicar las técnicas usuales de cuadrados mínimos (lineales). Esto puede lograrse, por ejemplo, considerando:
\begin{equation}\label{problema2}
\widetilde{F}(a,b,r) = \sum_{i=1}^n \tilde{d}_i^2 = \sum_{i=1}^n  \Big((x_i-a)^2+(y_i-b)^2-r^2\Big)^2
\end{equation}
Este enfoque, sin embargo, no aporta ningún avance significativo:
\bej Probar que minimizar $\widetilde{F}$ es equivalente a resolver el ajuste algebraico \eqref{algebraico}.
\fej



Para resolver \eqref{geometrico} utilizaremos el método iterativo de Gauss-Newton, que describimos a continuación.
El objetivo es minimizar $$\sum_{i=1}^n d_i(\theta)^2, \quad \theta = (a,b,r).$$
Notamos ${\bf d}(\theta) = (d_1(\theta),\dots,d_n(\theta))^t$.
Idealmente (si los datos coincidieran en un mismo círculo) tendríamos que la solución optima $\theta^*$ satisface ${\bf d}(\theta^*) = {\bf 0}$.

Desarrollando el polinomio de Taylor a orden 1 tenemos que:
$${\bf d}(\theta+{\bf h}) \simeq {\bf d}(\theta)+J(\theta){\bf h}\sim {\bf 0},$$
donde $J$ es la matriz diferencial de ${\bf d}$. El método de Gauss-Newton toma una solución inicial $\theta_0$, resuelve el problema de cuadrados mínimos:
\begin{equation}\label{sistema}
J(\theta_0){\bf h} = -{\bf d}(\theta_0), 
\end{equation}
toma $\theta_1 = \theta_0 + {\bf h}$, e itera el procedimiento, calculando una sucesión $\theta_i$ y deteniéndose cuando se satisface algún criterio de convergencia. 

\bej Estudiar, analíticamente, la iteración que daría el método de Newton-Raphson para el problema \eqref{geometrico}. Comparar con el método de Gauss-Newton descripto anteriormente. ¿Qué ventajas presenta Gauss-Newton? \fej 

\bej Calcular analíticamente $J(\theta)$ para el problema \eqref{geometrico} e implementar el algoritmo de Gauss-Newton, resolviendo \eqref{sistema} a través de la descomposición en valores singulares de $J$, $J=U\Sigma V^t$. (Usar el comando \verb+svd+ de Matlab). \fej

\subsection*{Datos singulares}
Cuando la matriz $\Sigma$ es casi singular el error numérico del método anterior crece, haciendo que la solución resulte poco confiable. Es importante remarcar que este inconveniente \emph{no} proviene del uso de la descomposición en valores singulares, sino que es intrínseco a la matriz $J$ y, más esencialmente, a los datos {\bf x} e {\bf y}. El hecho de que $\Sigma$ sea casi singular es un síntoma de que $J$ está cerca de las matrices de rango menor que 3, lo que hace que el problema \eqref{sistema} no tenga solución única. Para resolver este inconveniente Levenberg y Marquardt propusieron una modificación del algoritmo de Gauss-Newton. La idea es la siguiente: se elige un parámetro $\lambda$ y se resuelve \eqref{sistema} reemplazando:
$$\Sigma = \Sigma + \lambda I.$$
De este modo se \emph{corrige} $\Sigma$ para salvar su (casi) singularidad. Cuanto mayor es el valor de $\lambda$, más se parece el vector ${\bf h}$ resultante al dado por el método del gradiente; cuanto menor es $\lambda$, más parecida es la iteración a la correspondiente al método de Gauss-Newton puro. La elección de $\lambda$ no resulta en absoluto trivial. Marquardt sugirió una mecánica iterativa, que actualiza el valor de $\lambda$ en cada paso, según el comportamiento de $J$: 

Se toma un valor inicial $\lambda_0$ (por ejemplo $\lambda_0 = 1$) y un parámetro de correción $\alpha>1$. Se considera la solución inicial $\theta=(a,b,r)$ y se calculan los valores $\theta'=(a',b',r')$ y $\theta''=(a'',b'',r'')$, correspondientes a realizar una iteración del método con $\lambda = \lambda_0$ y con $\lambda = \frac{\lambda_0}{\alpha}$, respectivamente. Para todas estas soluciones se computa el valor del funcional $F$. Si tanto $\theta'$ como $\theta''$ arrojan valores de $F$ \emph{mayores} que el correspondiente a $\theta$, entonces se agranda $\lambda$ iterando: $\lambda = \alpha\lambda$ hasta hallar un $\lambda$ que mejore el valor del funcional. Si, en cambio, $F(\theta'')<F(\theta)$, se actualiza $\lambda = \frac{\lambda}{\alpha}$, mientras que si $F(\theta')<F(\theta)$, se conserva el valor de $\lambda$ actual. 

\bej Implementar un programa que ejecute el método de Gauss-Newton en la variante de Levenberg-Marquardt. 
\fej 

Para probar todos los algoritmos implementados, en la página de la materia está disponible el archivo \verb+datos.mat+. Al cargar \verb+datos.dat+ en Matlab se generan varias matrices de nombre: \verb+datos_i+ ($i=1,\dots,10$). Cada una de estas matrices de $n_i\times 2$ corresponde a un ejemplo y contiene en su primer columna los valores de ${\bf x}$ y en la segunda los valores de ${\bf y}$.

\subsection*{Yapa}
\bej Probar los algoritmos con los datos propuestos.
\begin{itemize}
 \item ¿Qué se observa? ¿Para qué conjuntos de datos el ajuste algebraico resulta \emph{aceptable}? ¿Qué particularidad tienen los datos cuyo ajuste \emph{no} es bueno?
 \item Observar que en \eqref{algebraico} se está minimizando $\sum_i(r_i^2-r^2)=\sum_i(r_i-r)(r_i+r)$, donde $r_i$ es la distancia del $i-$ésimo dato al centro y $r$ es el radio del círculo. La solución al problema deberá equilibrar los aportes de los dos factores involucrados en esta suma.
 \item Dado un conjunto de puntos, ¿cómo conviene elegir el centro del círculo para minimizar el factor $(r_i+r)$? (No hace falta hacer cuentas: se espera una respuesta cualitativa). 
 \item ¿Qué explicación podría darse para el ajuste algebraico en los casos en que da mal?
\end{itemize}
\fej 

 
 \end{document}

 
 