INTRODUCCIÓN AL AJUSTE DE DATOS

(DATA FITTING)

Germán A. Holguín L.
Profesor Titular Ingeniería Eléctrica
Grupo de Investigación en Gestión de Sistemas Eléctricos, Electrónicos y Automáticos GSEEA
Universidad Tecnológica de Pereira

10 de marzo de 2025

Este documento presenta un conjunto de herramientas matemáticas básicas y necesarias

para resolver el problema de minimización de mínimos cuadrados en el contexto de la
aplicación a la estimación de parámetros y ajuste de modelos matemáticos para datos
experimentales, aún en presencia de ruido y datos atípicos.

Índice
1 SVD (Singular Value Decomposition) 3
1.1 Valor Singular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Caso Real. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 Relación con los valores propios . . . . . . . . . . . . . . . . . . . . . . 4

2 MÍNIMOS CUADRADOS LINEALES 5

2.1 Sistemas Lineales de Ecuaciones . . . . . . . . . . . . . . . . . . . . . . 5
2.1.1 Caso de rango completo . . . . . . . . . . . . . . . . . . . . . . 5
2.1.2 Caso de rango deficiente . . . . . . . . . . . . . . . . . . . . . . 6
2.1.3 Caso de rango desconocido . . . . . . . . . . . . . . . . . . . . . 7
2.2 Pseudo-Inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Solución mediante ecuaciones normales . . . . . . . . . . . . . . . . . . 8
2.4 Ecuaciones Homogéneas . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1
MIE - UTP GSEEA, UTP.

3 RANSAC (Random Sample Consensus) 11

3.1 Definiciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Concepto Básico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Un primer ejemplo conceptual: Modelo Lineal . . . . . . . . . . . . . . 12
3.4 Análisis de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.5 Algoritmo RANSAC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

4 MíNIMOS CUADRADOS NO LINEALES 17

4.1 Gradiente Descendente. GD. . . . . . . . . . . . . . . . . . . . . . . . . 18
4.1.1 Cálculo del gradiente . . . . . . . . . . . . . . . . . . . . . . . . 18
4.1.2 Criterios de parada . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.1.3 Observaciones y conclusiones sobre GD. . . . . . . . . . . . . . 19
4.2 Gauss-Newton. GN. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2.1 Formulación final de GN. . . . . . . . . . . . . . . . . . . . . . . 22
4.2.2 Observaciones y conclusiones sobre GN. . . . . . . . . . . . . . 22
4.3 El método de Levenberg-Marquardt. LM. . . . . . . . . . . . . . . . . . 23
4.3.1 Invariante de ciclo . . . . . . . . . . . . . . . . . . . . . . . . . 23
4.3.2 Test de Calidad ρLM
k+1 . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.3.3 Demostración . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.3.4 Interpretación de ρLM
k+1 . . . . . . . . . . . . . . . . . . . . . . . . 25

Página 2 de 25
MIE - UTP GSEEA, UTP.

1 SVD (Singular Value Decomposition)

La descomposición de valor singular, SVD, es una factorización de matrices reales Rm×n
o complejas Cm×n de la forma
H
Am×n = Um×m Σm×n Vn×n (1)

donde V H es la transpuesta conjugada de V , también conocida como la Matriz Hermi-

tania de V . A es una matriz rectangular real o compleja. U y V son matriz unitarias.
La matriz Σ es una matriz diagonal de valores singulares.
Las matrices unitarias son equivalentes a la versión compleja de las matrices ortogonales
reales. Esto es,
UHU = UUH = I
V HV = V V H = I

1.1 Valor Singular

Los valores singulares de A son números reales no-negativos σ para los cuales existen
vectores de norma unitaria ~u ∈ Cm y ~v ∈ Cn tal que

A~v = σ~u

AH ~u = σ~v

donde ~v se conoce como el vector singular derecho de σ, y ~u se conoce como el vector

singular izquierdo de σ.
La matriz U contiene en sus columnas todos los vectores singulares izquierdos de A y
la matriz V contiene en sus columnas todos los vectores singulares derechos de A.
Por convención, los valores singulares se ubican en la diagonal de la matriz Σ y se
ordenan de mayor a menor. Esto quiere decir que

σ1 ≥ σ2 ≥ σ3 ≥ ... ≥ σn

La columna i de U es el vector singular izquierdo ui que está asociado al valor singular

σi . De la misma forma, el vector columna vi de V es el vector singular derecho asociado
con σi .
Una matriz A de m × n tiene como máximo p = min(m, n) valores singulares distintos.
Un valor singular para el que existan dos vectores izquierdos o derechos, es un valor
singular degradado. Cuando una matriz tiene σi degradados, entonces su factoriza-
ción SVD NO es única.

Página 3 de 25
MIE - UTP GSEEA, UTP.

1.2 Caso A ∈ Rm×n

Cuando la matriz A es real, la descomposición de valor singular se convierte en la forma

T
Am×n = Um×n Dm×n Vn×n (2)

Donde U y V son matrices ortogonales y D es una matriz diagonal.

1.3 Relación con los valores propios

Es importante evitar confundir los valores singulares σ con los valores propios λ.

A = U DV T
T
A A = V DU T U DV T
AT A = V D2 V T
VT = V −1
AT A = V D2 V −1

Por definición, esto es la descomposición de valores propios de AT A. Por tanto, los

valores propios de AT A son las entradas de D2 y las columnas de V son los vectores
propios de AT A.
Los valores singulares de A son la raíz cuadrada de los valores propios de AT A.
Como AT A es simétrica y positiva semi-definida, entonces sus valores propios son reales
y no-negativos. Por tanto, la raíz cuadrada de los valores propios, los valores singulares,
serán valores reales y no-negativos.

Página 4 de 25
MIE - UTP GSEEA, UTP.

2 MÍNIMOS CUADRADOS LINEALES

El método de los mínimos cuadrados lineales se utiliza para encontrar una aproximación
al problema de los sistemas lineales sobre-determinados. De forma clásica, la solución a
este problema implica invertir una matriz de coeficientes, proceso que en general es de
alto costo computacional y altamente vulnerable en sistemas pobremente condicionados.
Sin embargo, en las aproximaciones modernas se imponen los métodos de proyección
ortogonal que no involucran el cálculo explícito de la inversa. Por esta misma razón,
este resumen privilegia el método de la descomposición de valor singular, que además
generaliza adecuadamente a los casos homogéneos y con restricciones.

2.1 Sistemas Lineales de Ecuaciones

Considere un sistema de ecuaciones de la forma

Am×n xn×1 = bm×1 (3)

Donde x ∈ Rn , b ∈ Rm , y la matriz A es una transformación lineal T : Rn → Rm .
Por tanto, dependiendo de los valores específicos de m y n, existen entonces tres posi-
bilidades:

1. m < n. Más incógnitas que ecuaciones.

Sistema sub-determinado. No existe una solución única. Existe un espacio vecto-
rial de soluciones.
2. m = n. Igual número de incógnitas y ecuaciones.
El sistema es determinado (rango completo) si y sólo si det A 6= 0, caso en el cual
el sistema tiene solución única.
3. m > n. Más ecuaciones que incógnitas.
Sistema sobre-determinado. En general, NO TIENE SOLUCIÓN (a menos que
b ∈ span(col(A))).

2.1.1 Caso de rango completo

Cuando m ≥ n, y el rango de A es completo, rank(A) = n, entonces una solución de

mínimos cuadrados está dada por
x∗ = arg mı́nn kAx − bk
x∈R

Sea A = U DV T la SVD de A, por tanto

kAx − bk = kU DV T x − bk
= kU T U DV T x − U T bk
= kDV T x − U T bk

Página 5 de 25
MIE - UTP GSEEA, UTP.

Lo cual se puede re-escribir de forma más simple, si se definen:

y , V Tx (4)

y
b0 , U T b (5)

El problema se convierte entonces en

arg mı́nkDy − b0 k
y

Entonces el nuevo sistema a resolver es de la forma

 0 
  b1
d1  b02 
d  . 
 . 
 2
  
 . 
  y1
 d 3

  y2   0 
...     bn 

  ..  = 

 .  
 
 

 d n 
 yn
 0 
bn+1 
 . 
 .. 
 
0
b0m

Evidentemente, el vector Dy mas cercano a b0 es el vector

 0 T
b1 , b02 , b03 , . . . , b0n , 0, . . . , 0

que se obtiene cuando

b0i
yi = i = 1..n di 6= 0
di

Para que finalmente, de la ecuación (4) se tiene que

x=Vy

El algoritmo 1 resume entonces el método para calcular x utilizando SVD.

2.1.2 Caso de rango deficiente

Cuando m > n y el rango de A no es completo, r = rank (A) < n, existe un espacio

vectorial de soluciones, parametrizado con (n − r) parámetros λ. Por tanto, la solución
es de la forma
x = V y + λr+1 V~r+1 + ... + λn V~n
Donde V~r+1 hasta V~n son las (n − r) últimas columnas de V .
En resumen, el procedimiento para este caso es como se muestra en el algoritmo 2.

Página 6 de 25
MIE - UTP GSEEA, UTP.

Algoritmo 1 Resumen caso de rango completo.

Require: rank(A) = n
1: procedure SolverRangoCompleto(A, b)
2: U, D, V ← svd(A)
3: b0 ← U T b
4: for i ← 1, n do
5: y(i) = b0 (i)/d(i)
6: x←Vy
7: return x

Algoritmo 2 Resumen caso de rango deficiente.

Require: r = rank(A) < n
1: procedure SolverRangoDeficiente(A, b)
2: U, D, V ← svd(A)
3: b0 ← U T b
4: y = ~0
5: for i ← 1, r do
6: y(i) = b0 (i)/d(i)
7: x←Vy
8: return x

2.1.3 Caso de rango desconocido

En la práctica es el rango de la matriz A es casi siempre estimable. Sin embargo, si por

alguna razón no lo es, en todo caso se debe aproximar. Por ejemplo, aquellos valores
singulares que sean muy pequeños, en comparación con aquellos más grandes, pueden
hacerse cero.
σi
< δ ⇒ yi = 0
σ1

2.2 Pseudo-Inversa
Dada una matriz diagonal cuadrada D, se define su Pseudo-Inversa D+ tal que
(
0 di = 0
Di+ =
1/di di 6= 0

Dada una matriz Am×n con m ≥ n donde A = U DV T , se tiene que

A+ = V D+ U T

Página 7 de 25
MIE - UTP GSEEA, UTP.

Y de las ecuaciones (4) y (5) se tiene que

A+ b = V D+ U T b
= V D + b0
=Vy
=x

En conclusión,
x = A+ b (6)

Por tanto, x es la solución de mínimos cuadrados para Ax = b cuando rank(A) = n.

Para los casos de rango deficiente, x = A+ b es la solución que minimiza kxk. Y para
los casos donde m < n, simplemente se extiende A agregando filas de ceros hasta que
A sea cuadrada.

2.3 Solución mediante ecuaciones normales

El problema de mínimos cuadrados lineales puede ser solucionado utilizando una me-
todología que involucra las denominadas ecuaciones normales. El problema en cuestión
es el mismo sistema sobre-determinado que se ha estado estudiando antes:

Ax = b A ∈ Rm×n ; m>n

Como este sistema no tiene solución en general, la tarea consiste en encontrar un valor
x∗ tal que
x∗ = arg mı́nn kAx − bk
x∈R

Nótese que x ∈ Rn , pero Ax ∈ Rm . Ax debe estar en el sub-espacio generador de las

columnas de A. Esto es
Ax ∈ span(Col(A)).
Desde este punto de vista, la tarea se convierte en encontrar el vector Ax, más cercano
a b que pertenezca al span(Col(A)).
Si x es una solución, luego, la diferencia Ax − b debe ser ortogonal con Col(A). Otra
forma de decirlo es que Ax − b es perpendicular a cada una de las columnas de A. Al
ser ortogonal, su producto punto debe ser cero, así:

AT (Ax − b) = 0

Luego,
AT Ax = AT b

AT A x = AT b (7)

Página 8 de 25
MIE - UTP GSEEA, UTP.

La ecuación (7) representa lo que se conoce como un sistema de ecuaciones normales. Es

decir, un sistema lineal de n×n de rango completo. Se puede demostrar que rank(A) = n
si y solo si rank(AT A) = n. Por tanto, AT T es invertible. En conclusión:

−1
x = AT A AT b (8)

Si se comparan las ecuaciones (6) y (8) tenemos la ecuación (9), más conocida como la
pseudo-inversa de Moore-Penrose.

−1 T
A+ = AT A A (9)

El cálculo de la pseudo-inversa de Moore-Penrose involucra invertir una matriz, por lo

tanto, este método de las ecuaciones normales es computacionalmente más eficiente sólo
cuando el valor de n sea relativamente pequeño y la matriz AT A sea bien condicionada.

2.4 Ecuaciones Homogéneas

En los sistemas de ecuaciones homogéneos, el vector b = 0.

Ax = 0

Para m > n y teniendo en cuenta que x = ~0 (solución trivial) no es una solución de

interés, obsérvese que si x es una solución, entonces kx también lo es, ∀k ∈ R.
Se requiere entonces establecer una, o un conjunto de restricciones en x, que permitan
determinar el mejor x. Una restricción razonable es kxk = 1.

arg mı́n kAxk

kxk=1

Utilizando un procedimiento similar al visto en la sección 2.1, y sea A = U DV T ,

entonces

kAxk = kU DV T xk
= kDV T xk

De nuevo, y , V T x, y sabiendo que la norma kV yk es igual a la norma kyk, entonces

el problema queda efectivamente transformado en

arg mı́n kDyk

kyk=1

Donde D es una matriz diagonal con entradas no-negativas en orden descendiente. El

vector y que minimiza el producto kDyk estará dado por y = (0, 0, 0, ..., 1)T . Esto se

Página 9 de 25
MIE - UTP GSEEA, UTP.

Algoritmo 3 Resumen caso homogéneo.

Require: rank(A) = n
1: procedure solverHomogeneo(A)
2: V ← svd(A)
3: x ← V (lastColumn)
4: return x

debe a que el último elemento de la diagonal de D es el valor singular más pequeño de

A. Esta solución cumple de hecho con kyk = 1. Luego,

x=Vy

Pero y = (0, 0, ..., 1)T significa que x es la ultima columna de V .

Esto tiene grandes implicaciones computacionales. Quiere decir que para calcular x en
un sistema homogéneo, sólo es necesario calcular V . No es necesario calcular U . Y la
respuesta es simplemente la última columna de V . Esto se resumen en el algoritmo 3.

Página 10 de 25
MIE - UTP GSEEA, UTP.

3 RANSAC (Random Sample Consensus)

El método lineal de mínimos cuadrados sólo puede tolerar una pequeña cantidad de
ruido en los datos. Es decir, la minimización algebraica

arg mı́n kAxk

kxk=1

es sensible a la varianza del ruido que contamina las muestras. Adicionalmente, si ade-
más existen correspondencias putativas falsas, el método de mínimos cuadrados sim-
plemente arrojará resultados equivocados.
Las correspondencias putativas entonces se pueden dividir en dos grandes categorías, las
típicas (inliers), y las atípicas (outliers). Las correspondencias atípicas podrían además
incluir correspondencias falsas, conocidas también como atípicas extremas.
Los métodos de “estimación robusta” buscan precisamente obtener el modelo que
mejor explica los datos, aún en presencia de dichos atípicos. La técnica más común
implica identificar los datos atípicos para poder eliminarlos de la regresión del modelo.
Existen varios métodos de detección de atípicos (outlier detection methods). Cuando se
procesan los datos de forma recursiva, el método más simple se conoce como z-score.
Este método calcula para cada dato entrante x, cual es su factor de la desviación
estándar σ, dada una media estimada µ,
x−µ
z=
σ
y rechaza aquellos datos que superan un umbral τ establecido. Esto es

x ∈ N (µ, σ) ⇐⇒ z ≤ τ

Por ejemplo, si τ = 3 el dato x supera 3σ, entonces x se considera atípico y se rechaza.

Otra forma de decirlo es, el dato x está “más lejos” de µ que 3σ y por tanto es atípico.
Esta técnica simple, sin embargo, no es escalable a modelos complejos, o a distribuciones
de probabilidad multimodales o no-paramétricas.
En 1981, Fischler y Bolles introdujeron el algoritmo de Consenso por Muestreo Aleatorio
(RANSAC por sus siglas en inglés), que es un método iterativo para detectar datos
atípicos. Inicialmente concebido para el procesamiento de imágenes digitales, RANSAC
demostró rápidamente su generalidad y es hoy por hoy considerado estado del arte para
un rango muy amplio de aplicaciones.

3.1 Definiciones
Símbolo D: Conjunto de datos (correspondencias) disponibles.

Símbolo ηt : Número total de correspondencias en D. ηt = |D|.

Página 11 de 25
MIE - UTP GSEEA, UTP.

Inlier: Dato típico. Es un dato que pertenece al conjunto de datos explicados por el
modelo.

Outlier: Dato atípico. Es un dato que pertenece al conjunto de datos que no pueden
ser explicados por el modelo. Posiblemente datos falsos, o producto de errores de
adquisición, estimación o transmisión.

Símbolo n: En minúscula, es el número mínimo de datos necesario para estimar el

modelo.

Símbolo Ki : Sub-conjunto de datos Ki ⊂ D compuesto por n elementos. Esto es,

|Ki | = n. Este es el conjunto mínimo necesario para regresar un modelo.

Soporte S: Conjunto de datos S ⊂ D, que son inliers. También conocido como “con-
senso”.

Símbolo M : Mínimo número aceptable de inliers.

Símbolo N : En mayúscula, es el número de iteraciones máximo que se repite el pro-

ceso.

Símbolo δ: Tolerancia para aceptar un dato como inlier.

3.2 Concepto Básico

1. Seleccionar aleatoriamente un conjunto mínimo de n datos necesario para cons-
truir un modelo estimado a partir de ellos.

2. Dado el modelo estimado, determinar cual es el “soporte” para ese modelo. Esto
es, determinar |S|.

3. El modelo se acepta, cuando el “soporte” supere el umbral deseado |S| > M , o

cuando se hayan excedido N intentos.

4. Re-estimar el modelo utilizando todos los inliers en el S aceptado.

3.3 Un primer ejemplo conceptual: Modelo Lineal

Suponga que se desea estimar la relación lineal existente entre una variable aleatoria
de entrada y otra variable aleatoria de salida, ambas continuas. En este caso, regresión
lineal consiste en estimar la transformación afín de la forma

y = ax + b

La figura 1 muestra la idea básica de etiquetar los datos atípicos y separarlos de los que
si se ajustan al modelo.

Página 12 de 25
MIE - UTP GSEEA, UTP.

ideal
datos
outliers 150
ajuste

100

50

-10 -5 5 10

-50

Figura 1: Ajuste de un modelo linear en presencia de outliers.

Esta es una transformación afín, ya que en general

~y = Ay + ~b

y por tanto,
A ~b y
    
~y
= T
1 0 1 1

En este caso, el número mínimo de puntos necesario para construir un estimado del
modelo es 2. Por tanto, a partir de un conjunto de datos D dado

D = {(xi , yi ) : xi , yi ∈ R, i = 1, .., n}

RANSAC seguirá el siguiente procedimiento:

1. Se escogen aleatoriamente dos correspondencias (xi , yi ) del conjunto D. Estas dos

correspondencias constituyen Ki .

2. Se utiliza entonces un Ki para estimar los parámetros A y ~b del modelo.

Página 13 de 25
MIE - UTP GSEEA, UTP.

3. Se mide el soporte M del modelo.

4. Se repite este proceso hasta maximizar el número de inliers.

Para calcular el soporte, se utiliza para cada dato i alguna métrica de distancia, siendo
en este caso la norma L1 la más simple,

yi − axi − b = ξ | ξ |< δ

El valor de δ puede ser ajustado utilizando un modelo estadístico del ruido para los
inliers. En su forma más simple, podrá modelarse con una gaussiana de la forma
1 − 12 (4x2 +4y2 )
g(4x, 4y) = e 2σ
2πσ 2

Por ejemplo, si el valor deseado de inliers aceptado es del 99 %, entonces la elección

natural será δ = 3σ. Si se desea ser más estricto, 2σ dejará pasar sólo el 95 %, y 1σ el
68 % respectivamente.

3.4 Análisis de convergencia

Encontrar S a partir de D es un problema NP-Completo. En este tipo de problemas sólo
es posible conocer la respuesta óptima después de probar todas las posibles respuestas.
Si n es el mínimo número de datos para obtener un modelo y ηt es el número total de
datos disponibles, entonces, el número total de sub-conjuntos diferentes de n elementos
que se puede obtener a partir de D está dado por
 
ηt ηt !
= ηt Cn =
n (ηt − n)!n!

Para contextualizar este resultado, tres simples ejemplos:

ηt = 1000 n = 2 ηt Cn = 499500
ηt = 1000 n = 10 ηt Cn = 2.63 × 1023
ηt = 10000 n = 10 ηt Cn = 2.74 × 1033

En la práctica, no es posible comprobar un número tan elevado de combinaciones, por

lo que se espera que en algún algoritmo útil

N  ηt Cn

Pero también, N debe ser lo suficientemente grande tal que en al menos un intento Ki
no contenga outliers.
Para la tarea de determinar un valor adecuado de N para RANSAC, se define:

Página 14 de 25
MIE - UTP GSEEA, UTP.

ε: La probabilidad que un dato tomado al azar sea outlier.

Ejemplo: ε = 0.1 significa que se espera un 10 % de outliers.

ω: La probabilidad que un dato tomado al azar sea inlier ω = 1 − ε.

M ≈ ω ηt = (1 − ε) ηt

p: La probabilidad que al menos uno de los Ki no contenga outliers.

La probabilidad que al escoger aleatoriamente un sub-conjunto Ki de n corresponden-

cias, todas ellas sean inliers, es:
ω n = (1 − ε)n

La probabilidad que al menos una de las (n) correspondencias sea un outlier, es

1 − ω n = 1 − (1 − ε)n

La probabilidad que cada uno de los N intentos tenga al menos un outlier, es

[1 − (1 − ε)n ]N = [1 − ω n ]N

La probabilidad que al menos uno de los N intentos NO contenga outliers, es

p = 1 − [1 − (1 − ε)n ]N

Resolviendo por N , se tiene

1 − p = [1 − (1 − ε)n ]N

ln(1 − p) ln(1 − p)
N= n
=
ln [1 − (1 − ε) ] ln [1 − wn ]
Debido a que ε es en general desconocido, se puede estimar en cada iteración mediante

|Sb |
ε=1−
ηt
donde Sb es el mejor sub-conjunto S encontrado hasta la iteración i. La probabilidad
p se desea que sea muy alta, valores típicos están alrededor de 0.99. Cada vez que se
encuentra un consenso S más grande, N se reduce significativamente.
La figura 2 muestra como decrece N a medida que el soporte crece, para un caso
particular donde n = 2, p = 0.99.

Página 15 de 25
MIE - UTP GSEEA, UTP.

N
500

400

300

200

100

|Sb |
ηtotal
0.1 0.3 0.5 0.7 0.9

Figura 2: N vs Soporte

3.5 Algoritmo RANSAC

El algoritmo 4 resumen de la forma más general posible el proceso seguido por RANSAC
para marcar las correspondencias que pertenecen al modelo como inliers y descartar
los outliers.

Algoritmo 4 Resumen RANSAC.

1: procedure ransac(D, n)
2: i←0
3: N ← +∞
4: Sb ← ∅
5: while i < N and |Sb | < M do
6: Obtener Aleatoriamente Ki
7: Determinar Si
8: if |Si | > |Sb | then
9: Sb ← Si
10: Actualizar N
11: i←i+1
12: Estimar x con todos los inliers.
13: return x

Página 16 de 25
MIE - UTP GSEEA, UTP.

4 MíNIMOS CUADRADOS NO LINEALES

Es muy común que la hipótesis utilizada para relacionar las correspondencias de un
experimento o sistema, sea no lineal. En particular, que los parámetros a encontrar no
sean linealmente separables.
Sea D un conjunto de correspondencias, con |D| = N , de la forma
 
s m
D = (xi , yi ) ∀ i ∈ [0, N ) : x ∈ R , y ∈ R

obtenidos de un sistema multivariado F , como el representado por:

x y
y = F (x, p)
s m

donde la función vectorial F (x, p) ∈ Rm es en general no lineal y el vector p ∈ Rn es

el vector de parámetros que determinan F .
Entonces, se define el error individual, o residuo, de la correspondencia i, como

ri = yi − F (xi , p)

y en total, el vector de errores individuales como

r = [r0 , r1 , r2 , . . . , rN −1 ]T

y entonces, el error cuadrático total, está dado por

C (p) = rT r
N
X −1
= ri2
i=0

= k y − F (x, p) k2

La función C, es llamada función de costo, pues debe ser minimizada para encontrar
aquellos parámetros p que mejor explican las correspondencias (xi , yi ) bajo la hipótesis
F.
Razón por la cual, el problema de optimización asociado se expresa formalmente como

p∗ = arg mı́n C (p)

p

Luego C (p) define una superficie en el espacio (n + 1)-dimensional, C ∈ Rn+1 . Por lo

que se busca una manera de descender iterativamente por esta superficie, comenzando
desde un p0 hasta alcanzar el fondo del valle en p∗ .

Página 17 de 25
MIE - UTP GSEEA, UTP.

4.1 Gradiente Descendente. GD.

Dado el punto pk donde la altura de la superficie es C (pk ), el siguiente paso será

pk+1 = pk − γk ∇C (pk )

Donde ∇C (pk ) es el gradiente de la superficie en pk , y el factor γk representa el tamaño

del paso a dar en la iteración k en la dirección del gradiente negativo.
El signo negativo se requiere porque la dirección del gradiente se da en el sentido del
ascenso más pronunciado, pero se quiere avanzar un paso en el sentido del descenso más
pronunciado (steepest descent).
La implementación más simple utiliza un paso γk constante, que en otros escenarios
se conoce también como la tasa de aprendizaje (learning rate). Sin embargo, γ puede
utilizarse con algún comportamiento predefinido. Típicamente se define con un decai-
miento progresivo. De hecho, en caso que la función de costo crezca entre iteraciones,
esto es, cuando
C(pk+1 ) > C(pk )
esto quiere decir que γk era muy grande, y por tanto ese paso no puede darse, entonces
se divide γk a la mitad y se calcula nuevamente pk+1 .

4.1.1 Cálculo del gradiente

C (p) = r(p)T r(p)

N
X −1
= ri2 (p)
i=0

= k y − F (x, p) k2

∇C (p) = 2JrT (p) r (p)

donde Jr es el Jacobiano del vector residual r,

Jr = J (y − F (x, p))
= −JF

donde
∂f1 ∂f1
 
 ∂p1 · · · ∂pn 
 .. ... .. 
 
 .
JF =  . 

 
 ∂fm ∂fm 
···
∂p1 ∂pn

Página 18 de 25
MIE - UTP GSEEA, UTP.

4.1.2 Criterios de parada

Los criterios de parada para la familia de algoritmos de GD, se pueden clasificar en

absolutos, relativos e híbridos. Naturalmente se prefiere la clase híbrida por defecto, ya
que ofrece mayor flexibilidad y robustez.

Absolutos

a)
k pk+1 − pk k < δ; δ>0

b)
k C (pk+1 ) − C (pk ) k < ; >0

Relativos

a)
k pk+1 − pk k
, < δ; δ>0
k pk k

b)
k C (pk+1 ) − C (pk ) k
< ; >0
k C (pk ) k

Híbridos

a)
k pk+1 − pk k
, < δ; δ>0
máx (1, k pk k)

b)
k C (pk+1 ) − C (pk ) k
< ; >0
máx (1, k C (pk ) k)

4.1.3 Observaciones y conclusiones sobre GD.

A medida que el algoritmo se aproxima al mínimo, ∇C es cada vez más pequeño y por
tanto pk+1 − pk es también cada vez más pequeño, luego, la selección de γk es crítica
para alcanzar p∗ .

Página 19 de 25
MIE - UTP GSEEA, UTP.

4.2 Gauss-Newton. GN.

En 1809, Gauss introdujo una modificación al método de Newton con aplicación espe-
cífica al problema de los mínimos cuadrados no lineales. Gauss estuvo estudiando este
método ampliamente en sus labores de astronomía y seguimiento de trayectorias de
cuerpos celestes.
La idea de Newton fue que en funciones continuas, la forma más rápida de descender
sobre la superficie de búsqueda, NO necesariamente está en la dirección de la pendiente
más pronunciada. Newton propuso que, para llegar al mínimo, se puede saltar directa-
mente hasta el punto que minimiza la aproximación cuadrática de la función objetivo,
en una iteración, y que puede hallarse con la expansión en series de Taylor de segundo
orden en el punto pk .

40
T (2) 6
30
4
20 2
p∗ pk+1
pk
10
0.5 1 1.5 2 2.5 3
−2
−3 −2 −1 1 2 3 4
−10 T (1) −4

(a) Aproximaciones de Taylor. (b) Iteración hacia el mínimo.

f (xk )
f 0 (xk ) =
xk − xk+1
4
f (xk )
xk+1 = xk −
f 0 (xk )
Esta expresión equivale a encontrar
2
la raíz de T (1). Pero se desea encon-
trar donde T (2) es mínima, por tan-
to:
f 0 (xk )
xk+1 = xk − 00
2 2.2 2.4 pk 2.6 f (xk )
siempre y cuando
T (1)
−2 f 00 (xk ) > 0
(c) Raíz más cercana.

Figura 3: Principio de operación de Newton-Raphson.

Página 20 de 25
MIE - UTP GSEEA, UTP.

Joseph Raphson había propuesto el mismo método 45 años antes que la publicación ofi-
cial de Newton. Las publicaciones de Raphson no fueron tan conocidas en su momento.
Su trabajo fue reconocido posteriormente, aún cuando estaba enfocado principalmente a
polinomios y funciones escalares. El método, conocido ahora como de Newton-Raphson,
generaliza al caso Rn .

pk

x y

Para una función de costo C (pk ) con pk ∈ Rn , el método de Newton-Raphson se

formula según:
pk+1 = pk − H −1 ∇C(pk )

donde H es la matriz Hesiana de la función de costo C. Nótese que esto mismo se puede
reescribir como
pk+1 = pk + δpk

donde el paso δpk , en la iteración k, es calculado como la solución al sistema lineal de

la forma
Hδpk = −∇C(pk )

El problema para Gauss con este planteamiento era que se requiere el cálculo explicito
de la matriz Hesiana, lo cual no es trivial, o simplemente no es posible en el caso no
lineal más genérico en Rn . El objetivo es encontrar

p∗ = arg mı́n k y − F (x, p) k2

p

Gauss propuso modificar el cálculo de δpk notando que después de dar el paso δpk ,
F (pk + δpk ) debe estar más cerca al mínimo. Es decir, es una mejor aproximación de
y que F (pk ).
y ≈ F (pk + δpk ) = F (pk ) + r(pk )

Página 21 de 25
MIE - UTP GSEEA, UTP.

Donde el residuo r(pk ), puede escribirse como

r(pk ) = JF (pk ) · δpk

y de donde se puede encontrar δpk al solucionar el sistema lineal de ecuaciones resul-

tante. Nótese que r(pk ) es un vector columna conocido, así:

r(pk ) = y − F (x, pk )

Utilizando la pseudo-inversa de Moore-Penrose, la solución a este sistema está dada por

 −1
δpk = JF (pk )T JF (pk ) JF (pk )T r (pk )

4.2.1 Formulación final de GN.

 −1
pk+1 = pk + JF (pk )T JF (pk ) JF (pk )T r (pk )
Lo que se puede reescribir como

pk+1 = pk + δpk

donde el paso δpk se obtiene, de una mejor manera, resolviendo el sistema lineal de
ecuaciones  
JF (pk )T JF (pk ) δpk = JF (pk )T r (pk )

4.2.2 Observaciones y conclusiones sobre GN.

• La matriz JF (pk ) no puede ser de rango deficiente.

• Si JF (pk )T JF (pk ) es diagonal, entonces GN esta descendiendo por el camino del

gradiente mas pronunciado. Es decir, GN es equivalente a GD.

• La matriz JF (pk )T JF (pk ) debe ser siempre positiva definida.

JF (pk )T JF (pk ) > 0

• Por definición, GN asume que el valor de pk ya es de por si cercano a p∗ .

• Si el punto incial p0 es muy alejado, el método puede divergir.

Página 22 de 25
MIE - UTP GSEEA, UTP.

4.3 El método de Levenberg-Marquardt. LM.

LM trata de combinar las mejores características de GD con las mejores de GN. Esto
es, cuando el punto pk está retirado del objetivo, entonces LM se comporta más como
GD. Pero a medida que pk se acerca al objetivo p∗ , entonces LM se comporta más como
GN.
Para tal fin, LM introduce un término heurístico denominado factor de amortiguamien-
to, al sistema lineal que estima δpk , de la siguiente forma
 
JF (pk )T JF (pk ) + µk I δpk = JF (pk )T r (pk )

De esta forma, si µ = 0, entonces LM es equivalente a GN. Pero,

 cuando µ > 0  es
un valor mucho más grande que los elementos en la diagonal de JF (pk ) JF (pk ) ,
T

entonces LM se comporta más cercano a GD.

El reto es entonces encontrar un buen valor de µk en cada iteración, para lo cual se
introduce el concepto de test de calidad, ρLM
k+1 , y el algoritmo de LM.

Algoritmo 5 Método de Levenberg-Marquardt.

1: procedure LM
2: Determinar un p0
3: Determinar un µ0
4: while no se cumpla la condición de parada do
5: Resolver por δpk
6: Calcular pk+1 = pk + δpk
7: if ρLM
k+1 then . El Test de calidad fue exitoso
8: Calcular µk+1
9: else . El test de calidad ha fallado
10: pk+1 = pk
11: µk+1 = 2µk

Aquí hay algo de razonamiento circular. El valor de δpk depende del valor de µK . Sin
embargo, el mejor valor posible para µk depende del valor de δpk .

4.3.1 Invariante de ciclo

Las invariantes de ciclo (loop invariant) sirven para realizar un análisis de convergencia
de un algoritmo iterativo. En este caso, la invariante es que en cada iteración se produce
un valor más cercano al óptimo deseado.
El ciclo se inicializa con valores para p0 y µ0 . Al comienzo de cada iteración k se tienen
valores para pk y µk , con los que se procede a encontrar δpk , y dar el paso

pk+1 = pk + δpk

Página 23 de 25
MIE - UTP GSEEA, UTP.

Aquí se procede a evaluar la calidad del paso. Si el test de calidad ρLM

k+1 de δpk falla,
entonces se recalcula δpk haciendo:
pk+1 = pk
µk+1 = 2µk

De esta manera, el ciclo siempre se está moviendo en la dirección descendente. El ciclo

finaliza al cumplirse con la condición de parada.

4.3.2 Test de Calidad ρLM

k+1 .

El test de calidad de LM es la relación entre el cambio de la función de costo, y el

cambio producido por un valor particular de µk .
C(pk ) − C(pk+1 )
ρLM
k+1 =
δpTk JF (pk )T r(pk ) + δpTk µK Iδpk

Para el test de calidad, el signo de C(pk ) − C(pk+1 ) es muy importante. Si el signo es

negativo, probablemente el algoritmo pasó sobre el mínimo, entonces se debe regresar
y para el próximo paso usar una dirección más alineada con GD que con GN.

4.3.3 Demostración

Se sabe que
C(p) = rT (p)r(p)
y el cambio en dos iteraciones consecutivas esta dado entonces por
C(pk ) − C(pk+1 ) = rT (pk )r(pk ) − rT (pk+1 )r(pk+1 )

Sin embargo:
r(pk+1 ) = r(pk ) + Jr (pk ) δpk
r(pk+1 ) = r(pk ) − JF (pk ) δpk

Sustituyendo:
C(pk ) − C(pk+1 ) = 2δpTk JF (pk )T r(pk )

−δpTk JF (pk )T JF (pk ) δpk

+δpTk µK Iδpk − δpTk µK Iδpk

Recordando que
 
JF (pk ) JF (pk ) + µk I δpk = JF (pk )T r(pk )
T

Se obtiene:
C(pk ) − C(pk+1 ) = δpTk JF (pk )T r(pk ) + δpTk µK I δpk

Página 24 de 25
MIE - UTP GSEEA, UTP.

4.3.4 Interpretación de ρLM

k+1 .

Esta última ecuación muestra a la izquierda de igualdad el cambio de la función de

costo. A la derecha muestra la predicción del cambio hecha para un valor particular
de µk . Luego, la relación entre ambas, es la relación entre el cambio real obtenido y el
predicho por LM. Por tanto, se evalúan entonces 3 casos así:

k+1 > 0 : C(pk+1 ) < C(pk ). El paso se dio en la dirección correcta.

ρLM

k+1 ≥ 1 : Muy buena dirección. La reducción del costo fue mejor que la predicha.
ρLM

k+1 ≤ 0 : La predicción fue mala. La dirección es incorrecta.

ρLM

Para controlar la predicción, se puede utilizar la expresión

 
1 LM

3
µk+1 = µk · max , 1 − 2 ρk+1 −1
3

De esta forma:
1
ρLM
k+1 > 0 µk+1 = µk
3
ρLM
k+1 ≤ 0 µk+1 = 2µk
Luego, la única pregunta restante es como inicializar µk . Esto se puede hacer con el
valor más grande de la diagonal de JF (pk )T JF (pk ) que corresponde a la dirección en
la cual C(pk ) tiene un descenso más pronunciado.
n  o
T
µ0 = τ · max diag JF (pk ) JF (pk ) 0<τ ≤1

LM puede ser detenido con criterios de parado híbridos para:

• k JF (pk ) k
• k δpk k
• Número máximo de iteraciones

Referencias
[1] Edwin Kah Pin Chong and Stanislaw H. Żak, An introduction to optimization, Fourth edition.,
Wiley series in discrete mathematics and optimization, Wiley, Hoboken, New Jersey, 2013 (eng).
[2] Martin A. Fischler and Robert C. Bolles, Random sample consensus: A paradigm for model fitting
with applications to image analysis and automated cartography, Commun. ACM 24 (June 1981),
no. 6, 381–395.
[3] Richard Hartley and Andrew Zisserman, Multiple view geometry in computer vision, 2nd ed., Cam-
bridge University Press, USA, 2003.
[4] T. Strutz, Data fitting and uncertainty: A practical introduction to weighted least squares and
beyond (Springer, ed.), 2010.

Página 25 de 25