Minimos Cuadrados

Aproximaci´ on por m´ınimos cuadrados Juan Piccini 27 de octubre de 2003 Este material se basa en el libro ”Numerical Methods and Software”de David Kahaner, Cleve Moler y Stephen Nash, Ed. Prentice Hall International Introducci´ on Consideremos el siguiente ejemplo: colocamos una cierta cantidad de sal en un tanque vac´ıo, y luego comenzamos a llenarlo con agua. Cada cierto tiempo, digamos 5 segundos, medimos la concentraci´on de sal en el agua Es de esperar que la concentraci´on decrezca en forma lineal a medida que el tanque se va llenando (¿por qu´e?). Si llamamos f (t) a la concentraci´on de sal en el tanque medida en el tiempo t, tendremos una relaci´on del tipo f (t) = α+β t. Si graficamos concentraci´on versus tiempo, obtenemos la siguiente gr´afica

DATOS EXPERIMENTALES 82

80

78

Concentración

76

74

72

70

68

66

64

0

1

2

3

4

5

6

7

8 Tiempo

1

9

10

11

12

13

14

15

16

A partir de la misma, vemos que los puntos no est´an alineados. Esto era de esperar, puesto que es inevitable cometer errores en la medici´on, la mezcla de la sal en el agua puede no ser homog´enea, etc. Esta situaci´on es muy com´ un en todas las ramas de las ciencias e ingenier´ıas : efectuamos mediciones de un cierto fen´omeno, obteniendo datos y deseamos obtener un modelo que explique dichos datos. En el ejemplo anterior, dicho modelo nos permitir´a por ejemplo, predecir el valor de la concentraci´on en tiempos entre medidas. O puede suceder que disponemos de un modelo, y queremos someterlo a prueba (ver qu´e tan bien ajusta los datos provenientes de nuevos experimentos). Pero los datos siempre vienen con errores . ¿Qu´e se puede hacer al respecto?. Para ello, adoptaremos de aqu´ı en m´as una notaci´on conveniente, que nos permita un planteo lo m´as general posible. Supongamos que tenemos datos (t1 , b1 ), (t2 , b2 ), ..., (tm , bm ) correspondientes a una cierta funci´on f que es quien est´a ”detr´as”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´on, y b1 , ..., b15 los valores de la concentraci´on en dichos tiempos. En una situaci´on ideal, tendr´ıamos bi = f (ti )∀i = 1...m. Sin embargo, en la realidad lo que tenemos es bi ≈ f (ti ). Asumiremos que: observaci´on = 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´otesis habitual sobre los errores).

2

1.

M´ınimos cuadrados, caso lineal

M´as 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 coeficientes xj (en el ejemplo x1 = α y x2 = β) son los par´ ametros del modelo, los n´ umeros que querremos hallar. Como nuestro modelo depende de dichos par´ametros, 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´omeno puede explicarse por una funci´on f que es combinaci´on lineal de {φ1 , ..., φn }, y que los par´ametros del modelo son los coeficientes de dicha combinaci´on. 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 : x1 φ1 (t1 ) + x2 φ2 (t1 ) + · · · + xn φn (t1 ) ≈ b1 x1 φ1 (t2 ) + x2 φ2 (t2 ) + · · · + xn φn (t2 ) ≈ b2 ··· x1 φ1 (tm ) + x2 φ2 (tm ) + · · · + xn φn (tm ) ≈ bm o, en notaci´on matricial,     φ1 (t1 ) φ2 (t1 ) · · · φn (t1 ) x1 b1      φ1 (t2 ) φ2 (t2 ) · · · φn (t2 )   x2   b2      ···    ··· ··· ···     ···  ≈  ···     ··· ··· ···   ···   ···  ··· φ1 (tm ) φ2 (tm ) · · · φn (tm ) xn bm

       

esto es, Ax ≈ b, donde la matriz A = ((aij )) = φj (ti )∀i = 1...m, j = 1...n representa a la funci´on f (f (x) = Ax); postulada como explicaci´on de los datos, b es el vector de observaciones, x es el vector de par´ ametros. Obs.Notemos que m es el n´ umero de datos, mientras que n es el n´ umero de funciones-modelo. En general tendremos que m À n, por lo que el sistema ser´a sobredeterminado. Esto significa que en general el sistema Ax = b ser´a incompatible (en nuestro ejemplo, no existe una recta que pase por todos los puntos), y por tanto el vector de residuos r(x) = b − Ax tendr´a norma positiva. Busquemos entonces una recta que pase ”lo m´as cerca”posible de todos los

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´ametros 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 minPm imice j=1 [(b − Ax)j ]2 , o lo que es lo mismo, hallar x que minimice la funci´on k r(x) k22 . Como estamos minimizando una suma de cuadrados, el m´etodo se llamar´a de m´ınimos cuadrados. Usando notaci´on 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

concentracion

76

74 b8−b(t8) 72

70

68

66

64

0

1

2

3

4

5

6

7 8 tiempo

9

10

11

12

13

14

15

Antes de seguir, recordemos qu´e significa que un sistema lineal Ax = b sea incompatible. Cuando decimos que existe un vector x0 tal que Ax0 = b, esto significa que podemos expresar al vector b como combinaci´on lineal de las columnas de A, y la soluci´on x0 no es otra cosa que los coeficientes usados en dicha combinaci´on. Si el sistema es incompatible, entonces nuestro vector b no est´a en el subespacio generado por las columnas de A. Lo que buscaremos entonces, es, dentro de dicho subespacio, el vector bp ”m´as cercano ”(de norma m´ınima) a b. Como dicho vector bp es de la forma Ax para alg´ un x ∈ Rn , el bp buscado ser´a aquel tal que b − bp sea ortogonal al subespacio generado por las columnas de A 4

b

Subespacio generado por las columnas de A

b-Ax

Ax o

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 escalerizaci´on. Obs. Notemos que AT A es nxn, con n ¿ m. No depende del n´ umero de datos, sino del n´ umero de par´ametros del modelo. Sin embargo, estas ecuaciones tienen una desventaja importante: como el n´ umero de condici´on de AT A es aproximadamente el cuadrado del n´ umero de condici´on de A, si la matriz A tiene un n´ umero de condici´on alto, entonces el x hallado puede tener un error relativo grande (¿por qu´e ?). Veamos un ejemplo : Supongamos que ²mach = 10−6 y que cond(A) = 1000. Entonces el error relativo en la soluci´on x cuando aplicamos escalerizaci´on 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´on no tendr´ıa ning´ un d´ıgito correcto. Pero si cond(A) no

5

es muy alto, y ²mach es muy peque˜ no, 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 81.8504 78.6934 78.8205 77.4579 77.6739 76.2863 74.3694 72.0555 73.4642 71.3341 70.8463 70.3758 69.7654 68.2146 65.5288

A 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000

Tenemos que cond(A) = 19,3133 y cond(AT A) = 373,0033. Como ²mach ≈ 10−16 , entonces el error relativo no es algo que pueda preocuparnos en este caso. La matriz AT A queda : 15 120

120 1240

y el vector AT b queda : x 10^3 1.1067 8.5719 Resolviendo las ecuaciones normales, hallamos x 81.8402 -1.0072 Esto es, la recta buscada tiene ecuaci´on 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´on por Q o por QT preserva la norma eucl´ıdea. Esto u ´ltimo puede verse f´acilmente : kQxk22 = (Qx)T Qx = xT QT Qx = xT x = kxk22 A partir del teorema de Gram-Schmidt (ortogonalizaci´on de bases), puede probarse que dada Amxn , existen Qmxm ortogonal y Rmxn triangular superior, tales que A = QR. Esta factorizaci´on de la matriz A en t´erminos de las otras dos se conoce como descomposici´ on QR. Entonces, minkb − Axk22 = minkb − QRxk22 x

x

= minkQT b − Qt QRxk22 = minkQT b − Rxk22 x

x

Como R es triangular superior, todas à sus filas!desde la n+1 hasta laÃm son nulas, ! R1 (QT b)1 T de modo que si pensamos a R como , y a Q b = como 0 (QT b)2 entonces, minkQT b − Rxk22 = mink(QT b)1 − R1 xk22 + k(QT b)2 k22 x

x

Como el u ´ltimo sumando no depende de x, no tenemos control sobre ´el, de modo que minimizaremos el primer sumando. Para ello alcanza con resolver el sistema nxn R1 x = (QT b)1 para hallar el x buscado, donde R1 es la ”parte de arriba”de R , la submatriz tama˜ no nxn de R formada por las primeras n filas de R, y (QT b)1 es el vector nx1 formado por las primeras n filas de (QT b)1 . Las soluciones halladas por este m´etodo son m´as seguras que las obtenidas mediante las ecuaciones normales, y adem´as hallar la factorizaci´on no requiere m´as operaciones que la resoluci´on de las ecuaciones normales. Veamos un ejemplo. Supongamos que tenemos un modelo f (x, t) = x1 + x2 t + x3 t2 , con 3 par´ametros, y 5 observaciones, correspondiente a la distancia de frenado de un autom´ovil seg´ un su velocidad. Las matrices A y el vector b son:

7

A 1 1 1 1 1

1 2 3 4 5

1 4 9 16 25

b 15 40 70 130 180

mientras que las matrices Q y R son: Q -0.4472 -0.4472 -0.4472 -0.4472 -0.4472

-0.6325 -0.3162 0.0000 0.3162 0.6325

0.5345 -0.2673 -0.5345 -0.2673 0.5345

-0.0258 0.2809 -0.6882 0.6367 -0.2036

-0.3371 0.7414 -0.2017 -0.4725 0.2698

R -2.2361 0 0 0 0

-6.7082 3.1623 0 0 0

-24.5967 18.9737 3.7417 0 0

( como ejercicio escribir Ax − b). Resolviendo R1 x = (QT b)1 , tendremos x = 1.0000 7.7143 5.7143

Esto es, f (t) = 5,7143 + 7,7143 t + 1 t2 . Graficando junto con los datos, obtenemos:

8

distancia de frenado vs. velocidad 200

180

160

distancia de frenado

140

120

100

80

60

40

20

0

2.

0

1

2

3 velocidad

4

5

6

M´ınimos cuadrados, caso no lineal

Veamos ahora otro ejemplo : Se desea medir la cantidad de una cierta sustancia en sangre (por ejemplo un medicamento). Se sabe que este tipo de fen´omenos siguen la ley de decaimiento exponencial , f (x, t) = x1 exp(x2 t). Los datos son los siguientes : t 0 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000

b 2.0147 0.6689 0.3064 0.1808 0.0020 0.0564 0.0677 -0.0779 -0.0714 0.0288 -0.0199

Ahora el residuo es ri (x) = bi − f (x, t). El problema sigue siendo el mismo: hallar el vector 9

Cantidad de medicamento en sangre versus tiempo 2.5

Cantidad de medicamento

2

1.5

1

0.5

0

−0.5

0

2

4

6

8

10

Tiempo

x = argminkb − f (x, t)k22 = argminkr(x)k22 x

x

Pero ahora r(x) es una funci´on no lineal en x. Veremos el M´ etodo de Gauss-Newton, el cual consiste en linealizar el problema y aplicar lo visto para el caso lineal. Supongamos una aproximaci´on xk a la soluci´on. Desarrollando por Taylor, tendremos que: r(x) ≈ r(xk ) + Jr (xk )(x − xk ) donde (Jr (xk ))i,j = ∂ri (xk )∂xj es la matriz Jacobiana de r. Llamando ρk = k xk+1 − x , tendremos r(xk+1 ) = r(xk ) + Jr (xk )(ρk ) Buscamos ρk que minimice kr(xk+1 )k22 . Notemos que si llamamos b = r(xk+1 ), A = −Jr (xk ), ρk = x, estamos en el caso lineal. Aplicamos entonces m´ınimos cuadrados lineal para hallar ρk , entonces obtenemos xk+1 = xk + ρk . Luego repetimos todo : desarrollamos r(xk+1 ) por Taylor, hallamos ρk+1 que minimiza kr(xk+2 )k22 , resolviendo de nuevo un caso lineal, etc. etc. En cada paso, obtenemos un nuevo vector xN tal que kr(xN )k22 es menor que la correspondiente al anterior vector. ¿Cuando parar? Un criterio posible es detenernos cuando kr(xk+2 )k22 ya ”no decrece m´as”. En nuestro ejemplo, obtenemos f (x, t) = 2 exp(−t), cuya gr´afica vemos junto a los datos originales. 10

Medicamento en sangre vs./ tiempo. Curva de ajuste, modelo no lineal 2.5

2

Cantidad de medicamento

1.5

1

0.5

0

−0.5 −1

3.

0

1

2

3

4

5 Tiempo

6

7

8

9

10

11

M´ınimos cuadrados ponderados

Hasta ahora no hemos cuestionado los datos que uno intenta ajustar. Sin embargo, supongamos por un momento que las m observaciones b1 , b2 , ..., bm no son confiables en el mismo grado. Tal vez fueron tomadas en diferentes d´ıas con distintos apartatos, o por distintas personas. En astronom´ıa, por ejemplo, si los datos son distancias a estrellas o galaxias, naturalmente los datos correspondientes a objetos lejanos tendr´an mayor error que los correspondientes a objetos cercanos. Uno va a confiar m´as en los datos ”bien recogidos”que en los otros, m´as ”sospechosos”. Sin embargo, si hay alguna informaci´on contenida en los datos ”sospechosos”, no estamos dispuestos a ignorarla confiando solamente en los datos ”buenos”(desperdiciar informaci´on es casi tan malo como desperdiciar comida). El arreglo m´as sencillo es adjudicar diferentes pesos wi , i = 1, ..., m a las m observaciones bi , i = 1, ..., m, y elegir el x que minimice la suma ponderada de cuadrados 2 w12 (b1 − f (x, t1 )2 + .. + wm (bm − f (x, tm )2

De este modo, mientras mayor sea el peso adjudicado a un dato, el proceso de minimizaci´on tratar´a con mayor intensidad de hacer peque˜ na la distancia correspondiente. As´ı, si adjudicamos los pesos mayores a los datos m´as confiables, ´estos tendr´an mayor peso, y la curva resultante se ver´a m´as atra´ıda por dichos puntos. En general, si tenemos una estimaci´on ei del error de la observaci´on bi , entonces una opci´on puede ser wi = 1/ei , de este modo los 11

datos con menor error tendr´an mayor peso. Otras veces los datos vienen con la misma cantidad de d´ıgitos confiables, esto es, tienen el mismo error relativo. 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´en 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´on 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´omo cambia (si es que lo hace) el m´etodo de Gauss-Newton.

12