Minimos Cuadrados
Minimos Cuadrados
Minimos Cuadrados
Juan Piccini
27 de octubre de 2003
Introducción
DATOS EXPERIMENTALES
82
80
78
76
Concentración
74
72
70
68
66
64
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Tiempo
1
A partir de la misma, vemos que los puntos no están alineados. Esto era de
esperar, puesto que es inevitable cometer errores en la medición, la mezcla de
la sal en el agua puede no ser homogénea, etc. Esta situación es muy común
en todas las ramas de las ciencias e ingenierı́as : efectuamos mediciones de un
cierto fenómeno, obteniendo datos y deseamos obtener un modelo que explique
dichos datos. En el ejemplo anterior, dicho modelo nos permitirá por ejemplo,
predecir el valor de la concentración en tiempos entre medidas. O puede suceder
que disponemos de un modelo, y queremos someterlo a prueba (ver qué tan bien
ajusta los datos provenientes de nuevos experimentos). Pero los datos siempre
vienen con errores . ¿Qué se puede hacer al respecto?.
Para ello, adoptaremos de aquı́ en más una notación conveniente, que nos per-
mita un planteo lo más general posible.
Supongamos que tenemos datos (t1 , b1 ), (t2 , b2 ), ..., (tm , bm ) correspondientes a
una cierta función f que es quien está ”detrás”de los datos (y a la cual nos
gustarı́a conocer). En el ejemplo anterior, t1 , ..., t15 son los tiempos en los cuales
medimos la concentración, y b1 , ..., b15 los valores de la concentración en dichos
tiempos.
En una situación ideal, tendrı́amos bi = f (ti )∀i = 1...m. Sin embargo, en la
realidad lo que tenemos es bi ≈ f (ti ).
Asumiremos que: observación = modelo + error, o lo que es lo mismo,
bi = f (ti ) + ei ∀ i = 1, .., m, donde las ei son variables aleatorias iid ∼ N (0, σ 2 )
(que es la hipótesis habitual sobre los errores).
2
1. Mı́nimos cuadrados, caso lineal
Más especı́ficamente, asumiremos que el modelo es de la forma
f (ti ) = x1 φ1 (ti ) + x2 φ2 (ti ) + ... + xn φn (ti ). Las funciones φj , j = 1...n se llaman
funciones modelo. En nuestro ejemplo, bi ≈ f (ti ),
donde f (ti ) = x1 + x2 ti ∀i = 1...m, con φ1 (t) = 1 ∀t, φ2 (t) = t ∀t. Los coefi-
cientes xj (en el ejemplo x1 = α y x2 = β) son los parámetros del modelo,
los números que querremos hallar. Como nuestro modelo depende de dichos
parámetros, en realidad tendremos f (x1 , x2 , ti ) = x1 + x2 ti , y en general,
f (x1 , x2 , ..., xn , ti ) = x1 φ1 (ti ) + x2 φ2 (ti ) + ... + xn φn (ti ). Notemos que estamos
diciendo que el fenómeno puede explicarse por una función f que es combinación
lineal de {φ1 , ..., φn }, y que los parámetros del modelo son los coeficientes de
dicha combinación. Por eso diremos que tenemos un modelo lineal (lineal en
x). Por ejemplo, f (t) = x1 + x2 t, f (t) = x1 exp(t) + x2 sin(t) + x3 t2 son
ejemplos de modelos lineales, mientras que f (t) = x1 tx2 o f (t) = x1 exp(x2 t)
no lo son.
Volviendo al tema, si tenemos los datos (t1 , b1 ), ..., (tm , bm ) y el modelo
bi ≈ x1 φ1 (ti ) + ... + xn φn (ti ), entonces tendremos las ecuaciones :
3
puntos, en el sentido de obtener un residuo de norma mı́nima. Observemos que
r(x) guarda las distancias entre bi y f (x1 , x2 , ti ) (ver figura). Deberemos hallar
los parámetros x1 x2 tales que la suma de los cuadrados de las distancias de
cada punto a la recta sea mı́nima, esto es, hallar el vector x = (x1 , x2 ) que min-
Pm
imice j=1 [(b − Ax)j ]2 , o lo que es lo mismo, hallar x que minimice la función
k r(x) k22 . Como estamos minimizando una suma de cuadrados, el método se
llamará de mı́nimos cuadrados. Usando notación vectorial, esto es equivale
a resolver
minkb − Axk22 ≡ min(b − Ax)T (b − Ax)
x x
cuando usamos la norma 2 Euclı́dea.
Interpretación de los residuos como distancias de la realidad al modelo
82
b1−b(t1)
80
78
76
concentracion
74
b8−b(t8)
72
70
68
66
64
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
tiempo
4
b
Subespacio
generado por
las columnas
de A
b-Ax
Ax
Recordemos que para que un vector sea ortogonal a un subespacio, basta con
que sea ortogonal a un generador de dicho subespacio, de modo que el vector
bp = Ax que buscamos tiene que ser tal que b − bp sea ortogonal a todas las
columnas de A (que son las filas de la traspuesta de A , AT ). Esto nos lleva a
pedir AT (b − Ax) = 0 o lo que es lo mismo,
AT Ax = AT b
Estas son las llamadas ecuaciones normales, y puede resolverse por escaler-
ización.
Obs. Notemos que AT A es nxn, con n ¿ m. No depende del número de
datos, sino del número de parámetros del modelo. Sin embargo, estas ecuaciones
tienen una desventaja importante: como el número de condición de AT A es
aproximadamente el cuadrado del número de condición de A, si la matriz A tiene
un número de condición alto, entonces el x hallado puede tener un error relativo
grande (¿por qué ?). Veamos un ejemplo : Supongamos que ²mach = 10−6 y que
cond(A) = 1000. Entonces el error relativo en la solución x cuando aplicamos
escalerización al sistema Ax = b, puede llegar a ser tanto como 1000x10−6 =
10−3 , lo cual significa que unos tres dı́gitos deberı́an ser correctos. Por otro
lado, si resolvemos AT Ax = AT b, el error podrı́a llegar a ser 10002 x10−6 = 1,
de modo que la solución no tendrı́a ningún dı́gito correcto. Pero si cond(A) no
5
es muy alto, y ²mach es muy pequeño, entonces las ecuaciones normales ya dejan
de ser una mala idea. Veamos el ejemplo de la sal en el tanque. Los datos (vector
y matriz) son:
b A
81.8504 1.0000 1.0000
78.6934 1.0000 2.0000
78.8205 1.0000 3.0000
77.4579 1.0000 4.0000
77.6739 1.0000 5.0000
76.2863 1.0000 6.0000
74.3694 1.0000 7.0000
72.0555 1.0000 8.0000
73.4642 1.0000 9.0000
71.3341 1.0000 10.0000
70.8463 1.0000 11.0000
70.3758 1.0000 12.0000
69.7654 1.0000 13.0000
68.2146 1.0000 14.0000
65.5288 1.0000 15.0000
15 120
120 1240
y el vector AT b queda :
x 10^3
1.1067
8.5719
x
81.8402
-1.0072
Esto es, la recta buscada tiene ecuación f (t) = 81,84 − 1 t (de hecho es la
recta que hemos graficado antes). Cuando cond(A) es grande,necesitamos una
6
alternativa a las ecuaciones normales para hallar minkb − Axk22 . Recordemos
x
de los cursos de Algebra Lineal, que una matriz Qmx m tal que QT Q = Id se
llama matriz ortogonal . Entre otras cosas, esto implica que las columnas de
Q forman un conjunto ortonormal en Rm , y que la premultiplicación por Q o
por QT preserva la norma euclı́dea. Esto último puede verse fácilmente :
R1 x = (QT b)1
7
A b
1 1 1 15
1 2 4 40
1 3 9 70
1 4 16 130
1 5 25 180
Q
-0.4472 -0.6325 0.5345 -0.0258 -0.3371
-0.4472 -0.3162 -0.2673 0.2809 0.7414
-0.4472 0.0000 -0.5345 -0.6882 -0.2017
-0.4472 0.3162 -0.2673 0.6367 -0.4725
-0.4472 0.6325 0.5345 -0.2036 0.2698
R
-2.2361 -6.7082 -24.5967
0 3.1623 18.9737
0 0 3.7417
0 0 0
0 0 0
x =
1.0000
7.7143
5.7143
Esto es, f (t) = 5,7143 + 7,7143 t + 1 t2 . Graficando junto con los datos, obten-
emos:
8
distancia de frenado vs. velocidad
200
180
160
140
distancia de frenado
120
100
80
60
40
20
0
0 1 2 3 4 5 6
velocidad
t b
0 2.0147
1.0000 0.6689
2.0000 0.3064
3.0000 0.1808
4.0000 0.0020
5.0000 0.0564
6.0000 0.0677
7.0000 -0.0779
8.0000 -0.0714
9.0000 0.0288
10.0000 -0.0199
9
Cantidad de medicamento en sangre versus tiempo
2.5
1.5
Cantidad de medicamento
0.5
−0.5
0 2 4 6 8 10
Tiempo
10
Medicamento en sangre vs./ tiempo. Curva de ajuste, modelo no lineal
2.5
1.5
Cantidad de medicamento
0.5
−0.5
−1 0 1 2 3 4 5 6 7 8 9 10 11
Tiempo
11
datos con menor error tendrán mayor peso. Otras veces los datos vienen con
la misma cantidad de dı́gitos confiables, esto es, tienen el mismo error relati-
vo. En este caso es usual hacer wi = 1/bi (en tanto sea bi 6= 0). Lo bueno
de esto es que cualquier software que resuelva mı́nimos cuadrados sin pesos,
también puede resolver el problema ponderado. En el caso lineal, en lugar de
resolver x = argminkb − Axk22 , resolveremos x = argminkbw − Aw xk22 , donde
x x
(bw )i = wi bi i = 1, ..., m, y ((Aw ))j,k = wj Ajk . Por supuesto que la solución
va a cambiar en la medida que los pesos cambien. Si todos los wi valen uno,
tenemos de nuevo el caso ya visto anteriormente. ¿Y para el caso no lineal?
En lugar de resolver x = argminkb − f (x, t)k22 , tendremos que hallar el x que
x
miminiza (w1 (b1 − f (x, t1 )))2 + ... + (wm (bm − f (x, tm )))2 . Queda como ejercicio
investigar cómo cambia (si es que lo hace) el método de Gauss-Newton.
12