Story Transcript
Cap´ıtulo 8
El Problema de M´ınimos Cuadrados
8.1.
El problema de m´ınimos cuadrados
El problema del ajuste de datos; es decir, descubrir una funci´on matem´atica que pueda explicar de la mejor forma posible el comportamiento de alg´ un mecanismo o grupo de seres u objetos que puede ser medido, y del cual conocemos algunos datos (con sus posibles errores de medici´on), es un problema cl´asico y ha supuesto un reto para la comunidad matem´atica desde su planteamiento por Gauss y Legendre hacia 1800. En lenguaje de ´algebra lineal consiste en encontrar la soluci´on de un sistema lineal Ax b siendo A P Cmn con m ¥ n. En A y b se recogen todos los datos del experimento que se quieren ajustar. Enseguida veremos algunos ejemplos. Sabemos que el sistema tiene soluci´on si y s´olo si b P Im A, condici´on que dif´ıcilmente se cumple si m es mucho mayor que n. Si m ¡ n diremos que el sistema Ax b est´a sobredeterminado; y si no existe ning´ un x P Cn1 tal que Ax b lo que se 159
160
El Problema de M´ınimos Cuadrados
puede intentar es buscar un x de forma que el vector r
Ax b P Cm,
que se llama residuo o vector residual, sea lo “m´as peque˜ no” posible. Lo grande o peque˜ no que sea r lo mediremos mediante una norma. El problema de m´ınimos cuadrados consiste en encontrar x P Cn1 para que el vector residuo tenga la menor norma eucl´ıdea posible. Su formulaci´on precisa ser´ıa la siguiente: Problema de m´ınimos cuadrados: Sea F R o C. Dada una matriz A P Fmn y un vector b P Fm ,(m ¥ n) encontrar x P Cn1 para que }Ax b}2 sea m´ınimo. La elecci´on de la norma eucl´ıdea se puede justificar desde varios puntos de vista: hist´orico, geom´etrico, estad´ıstico, . . . ; pero sobre todo porque es la norma habitual, conduce a los algoritmos m´as sencillos y as´ı como todas las normas son funciones continuas de los vectores, la norma eucl´ıdea es, adem´as, diferenciable. Como, por a˜ nadidura, la funci´on }Ax b}2 alcanza su m´aximo absoluto en un m´aximo local, este m´aximo puede calcularse igualando a cero las derivadas parciales de dicha funci´on. Este proceso conduce a lo que se llaman las ecuaciones normales del problema de m´ınimos cuadrados. Estas ecuaciones nos las encontraremos m´as adelante procedentes de un contexto completamente diferente.
Ejemplo 8.1 El ejemplo m´as t´ıpico de ajuste de datos por m´ınimos cuadrados es el c´alculo del polinomio de interpolaci´ on: dados m puntos del plano pxi , yi q, i 1, . . . , m, de forma que xi xj para i j, se trata de encontrar un polinomio de grado a lo m´as m 1, ppxq a0 a1 x am1 xm1 , que pase por los n puntos.
Se trata, entonces, de encontrar los coeficientes a0 , a1 ,. . . , am1 , para que ppxi q yi ,
i 1, . . . , m.
Esto equivale a resolver el sistema lineal a0
a1 x i
am1 xim1
yi ,
i 1, . . . , m
(8.1)
8.1 El problema de m´ınimos cuadrados
161
cuya matriz de coeficientes es
1 1 x1 x21 xm 1 m 1 2 j 1 1 x2 x2 x2 A xi .. .. .. . . . .. . . . . m1 2 1 xm xm xm Esta es una matriz de Vandermonde cuyo determinante es det A
¹
pxi xj q.
i j
Por consiguiente A es invertible y el sistema (8.1) tiene una u ´nica soluci´on. En la Figura 8.1 se muestra la gr´afica del polinomio de interpolaci´on que pasa por los 11 puntos pi, 0q con i 5, 4, . . . , 1, 0, 1, . . . , 4, 5. Se trata de un polinomio de grado 10. A su derecha se ha escrito el c´odigo de MATLAB que obtiene los coeficientes del polinomio de interpolaci´on: c Azb y que da la gr´afica: plot(t, polyval(p,t),a,b,’r*’,’markersize’,10);. El comando polyval(p,t) devuelve un vector: el valor del polinomio p en las componentes del vector t. N´otese el uso de las sentencias fliplr y flipud para obtener la matriz de los coeficientes del sistema y los coeficientes del polinomio de interpolaci´on en la forma apropiada. La sentencia zeroaxes hace que los ejes cartesianos se dibujen en la forma que se acostumbra a utilizar en la pizarra: dos l´ıneas perpendiculares que se cortan en p0, 0q. 5 4 3 2 1 −6
−4
−2
2 −1 −2 −3
4
6
Figura 8.1: a=-5:5; A=fliplr(vander(a)); b=[0 0 0 1 1 1 0 0 0 0 0]’; c=Azb; p=flipud(c); t=linspace(-5.05,5.05); plot(t, polyval(p,t),a,b,’r*’,...; ’markersize’,10); zeroaxes
Indudablemente el polinomio de interpolaci´on se ajusta perfectamente a los datos, pero a menudo ´estos proceden de experimentos y lo que se pretende es buscar una
162
El Problema de M´ınimos Cuadrados
curva que interprete los datos obtenidos. Como las mediciones se realizan en momentos puntuales, tal gr´afica deber´ıa reflejar, salvo que se tengan motivos para pensar lo contrario, una comportamiento suave entre dos datos consecutivos, y no la fluctuaci´on que muestra el polinomio de interpolaci´on de nuestro ejemplo. Quiz´a un polinomio de menor grado pueda mostrar un mejor comportamiento. Tal polinomio no pasar´a por algunos de los puntos. Pero la medici´on de todo proceso real conlleva errores que podr´ıan explicar tal fen´omeno. Ejemplo 8.2 (Ajuste por m´ınimos cuadrados) Con los mismos datos del ejemplo anterior calcular el polinomio de grado 7 que mejor se ajusta a los datos en el sentido de los m´ınimos cuadrados. Siguiendo los mismos pasos que en el ejemplo anterior lo que buscamos es un polinomio de grado a lo m´as n 1 ¤ m 1 ppxq a0
a1 x
an1 xn1
tal que }ppxq y }2 sea m´ınimo. Aqu´ı, ppxq es el vector cuya i-´esima componente es ppxi q y el vector y es el que tiene por i-´esima componente yi . Ahora bien, si
1 x1 x21 x1n1 1 x x22 x2n1 2 A .. .. .. . . .. . .. . . 2 n1 1 xm xm xm
es la matriz de Vandermonde truncada y c pa0 , a1 , . . . , an1 q es el vector de los coeficientes de p tenemos que ppxq Ac. As´ı que se trata de hallar un vector c en el que se alcance m´ınn }Ax y }2 .
P
x C
En la Figura 8.2 se presenta la gr´afica producida por MATLAB del polinomio de grado 7 que mejor se ajusta a los datos pai , bi q en el sentido de los m´ınimos cuadrados. El problema de m´ınimos cuadrados correspondiente se calcula en el interior de la sentencia ployfit. El resto del c´odigo de MATLAB que se muestra es el mismo que en el apartado anterior. Una observaci´on final antes de abordar la soluci´on del problema de m´ınimos cuadrados y los algoritmos correspondientes. Una variante del problema de encontrar
8.2 La soluci´on del problema de m´ınimos cuadrados
Figura 8.2:
5
a=-5:5 A=fliplr(vander(a)) b=[0 0 0 1 1 1 0 0 0 0 0]’; p=polyfit(a,b,7); t=linspace(-5.3,5.3); plot(t, polyval(p,t),a,b,’r*’,... ’markersize’,10); axis([-6 6 -3 5]); zeroaxes
4 3 2 1 −6
−4
−2
163
2
4
6
−1 −2 −3
el polinomio de grado a lo m´as n 1 que mejor se ajusta a una nube de m puntos, pxi, yiq, distribu´ıdos en el plano es el de encontrar la funci´on φpxq c1 φ1 pxq
c2 φ2 pxq
cn φn pxq,
que mejor se ajusta a una de tales nubes de puntos. Aqu´ı φ1 , φ2 ,. . . , φn son funciones dadas, o de las que uno sospecha que una combinaci´on lineal de ellas puede ajustarse bien a los datos. En este caso el problema se reduce a calcular el vector c donde se alcanza el m´ınimo: m´ın }Ax y }2 ,
P
x C
siendo A P Fmn la matriz cuyo elemento en la posici´on pi, j q es φj pxi q y siendo y el vector cuyas componetes son y1 , . . . , ym .
8.2.
La soluci´ on del problema de m´ınimos cuadrados
En esta secci´on demostramos el resultado que nos da la soluci´on del problema de m´ınimos cuadrados. Geom´etricamente la situaci´on es muy simple (v´ease la Figura 8.3): el vector de Im A cuya distancia a b es m´ınima es la proyecci´on ortogonal de b sobre Im A. Y la distancia m´ınima es la norma del residuo. Esto es lo que nos dice el siguiente Teorema.
164
El Problema de M´ınimos Cuadrados
r=Ax0−b b
Ax
Im A
0
Figura 8.3: La soluci´on del problema de m´ınimos cuadrados Teorema 8.3 Sean A P Fmn , F R o C, y b P Fm1 , m ¥ n. Sea PA la proyecci´on ortogonal sobre Im A. Entonces, Ax0 cumple
}Ax0 b}2 m´ ın }Ax b}2 xPF n
si y s´olo si se cumple cualquiera de las siguientes condiciones que son equivalentes:
PAb. b Ax0 P pIm AqK . A Ax0 A b.
(i) Ax0 (ii) (iii)
Adem´as, la soluci´on x0 es u ´nica si y s´olo si rang A n. En la demostraci´on se usa el teorema de Pit´agoras: Si x, y P Cn son ortogonales entonces }x y }22 }x}22 }y }22 . En efecto, }x y }22 px y q px y q x x y x y y x x x y y }x}22 }y }22 , donde hemos usado que x y y x 0 porque x e y son ortogonales. Demostraci´ on.- Todo se basa en la demostraci´on de la idea geom´etrica expuesta m´as arriba: m´ınn }Ax b}2 }PA b b}2 .
P
x C
Vemaos que esto es as´ı. En efecto
}Ax b}22 }Ax PAb
PA b b}22
}Ax PAb}22 }PAb b}22
8.3 Algoritmos para calcular la soluci´on del problema de m´ınimos cuadrados
165
porque Ax PA b P Im A y PA b b pI PA qb P pImAqK y aplicando el Teorema de Pit´agoras. Por lo tanto, para todo x P Cn
}Ax b}2 ¥ }PAb b}2. Pero como PA b P Im A, existe x1 P Cn tal que Ax1 PA b; i.e., }Ax1 b}2 }PA bb}2 . Por lo tanto m´ın }Ax b}2 }PA b b}2 , tal y como dese´abamos mostrar. xPC n
Ahora, si Ax0 es el vector que hace m´ınima la distancia de Ax a b entonces ın }Ax b}22 }Ax0 b}22 }Ax0 PA b}22 }PA b b}22 . }PAb b}22 xm´ PC n
De aqu´ı deducimos que }Ax0 PA b}2
0; es decir, Ax0 PAb.
S´olo queda demostrar la equivalencia entre las tres condiciones. Por una parte
ñ b Ax0 b PAb pI PAqb P pIm AqK. Y si b Ax0 P pIm AqK entonces PA pb Ax0 q 0 de modo que PA b PA Ax0 . Pero como Ax0 P Im A tenemos que PA Ax0 Ax0 . Esto demuestra la equivalencia entre PA b Ax0
las condiciones (i) y (ii).
Finalmente, como pIm AqK Ker A tenemos que b Ax0 A pb Ax0 q 0; i. e., A Ax0 A b.
P pIm AqK si y s´olo si
Falta demostrar que la soluci´on del problema de m´ınimos cuadrados es u ´nica si y s´olo si rang A n. Ahora bien, rang A n si y s´olo si A A es invertible. Una forma de ver esto es, por ejemplo, la siguiente: rang A n si y s´olo si σn pAq 0. Como los valores singulares de A son las ra´ıces cuadradas positivas de los valores propios de A A, σn 0 si y s´olo si todos los valores propios de A A son distintos de cero; i.e. detpA Aq 0. Pero A A es invertible si y s´olo si el sistema A Ax A b tiene soluci´on u ´nica.
8.3.
Algoritmos para calcular la soluci´ on del problema de m´ınimos cuadrados
El Teorema 8.3 nos da las claves para calcular un vector x0 que solucione el problema de m´ınimos cuadrados. En primer lugar, el sistema A Ax A b recibe el nombre de
166
El Problema de M´ınimos Cuadrados
ecuaciones normales del problema de m´ınimos cuadrados. Es el sistema que aparece al calcular el m´ınimo local de la funci´on f pxq }Ax b}2 . Es decir, el sistema
Bf pxq 0, B xi
i 1, . . . , n
da lugar al sistema A Ax A b. Como la funci´on f es convexa, alcanza su m´ınimo absoluto en un m´ınimo local. Para resolver el sistema A Ax A b num´ericamente, tendremos en cuenta la estructura especial de la matriz A A: es herm´ıtica (sim´etrica en el caso real). Adem´as es definida positiva porque x A Ax }Ax}2 ¡ 0 para x 0. Ahora bien toda matriz herm´ıtica definida positiva admite una factorizaci´on de Cholesky (variante de la factorizaci´on LU cuando la matriz es sim´etrica o herm´ıtica). En MATLAB el comando chol(A) devuelve la factorizaci´on de Cholesky de A si ´esta es herm´ıtica definida positiva. Recordemos que una factorizaci´on de Cholesky de una matriz herm´ıtica definida positiva, A, es una factorizaci´on de la forma A LL siendo L una matriz triangular superior. Esta factorizaci´on es especialmente apropiada para resolver sistemas lineales mediante sustituci´on hacia adelante y hacia atr´as. Teniendo en cuenta todo esto podemos dar un primer algoritmo para la resoluci´on del problema de m´ınimos cuadrados.
Algoritmo m´ınimos cuadrados via ecuaciones normales Dada A P Fmn y dado b P Fm 1. F´ormense las matrices A A y A b. 2. Calc´ ulese la factorizaci´on de Cholesky de A A LL , (chol(A))
Ab (z o sustituci´on hacia adelante) Resu´elvase L x y.(z o sustituci´on hacia atr´as)
3. Resu´elvase Ly 4.
8.3 Algoritmos para calcular la soluci´on del problema de m´ınimos cuadrados
167
Cuando A es herm´ıtica definida positiva pero mal condicionada, el algoritmo de Cholesky puede dar resultados muy inexactos. Por lo tanto, el algoritmo anterior no es adecuado para resolver el problema de m´ınimos cuadrados. Una posibilidad ser´ıa usar la factorizaci´on QR de A A para conseguir la factorizaci´on de Cholesky de A A en vez del algoritmo de Cholesky. En efecto, si A QR es la factorizaci´on QR reducida de A, entonces A A R Q QR R R LL , donde hemos utilizado que Q Q In , por tener Q columnas ortonormales, y hemos puesto L R . Esto nos muestra que si A es de rango completo, la factorizaci´on de Cholesky de A A es u ´nica. Esta alternativa, sin embargo, carece de sentido porque la ventaja del algoritmo de Choleski es que su coste es n2 (mientras que el de la factorizaci´on QR ser´ıa de orden c´ ubico) y hay otros algoritmos que, cuando el de Choleski no es aconsejable, usan la factorizaci´on QR y s´olo precisan resolver un sistema triangular. Presentamos a continuaci´on uno de tales algoritmos. En el segundo algoritmo se trata de conseguir la soluci´on, x0 , como soluci´on del sistema Ax0 PA b. Para ello, debemos conseguir PA , que es la proyecci´on ortogonal sobre Im A. Recordemos que PA QQ , donde Q es una matriz cuyas columnas son una base ortonormal de Im A. Recordemos que si A QR es una factorizaci´on QR de A entonces las columnas de Q son una base ortonormal de Im A y, por consiguiente, QQ es la proyecci´on ortogonal sobre Im A.
Algoritmo m´ınimos cuadrados via factorizaci´ on QR Dada A P Fmn y dado b P Fm 1. H´allese la factorizaci´on QR, reducida, de A 2. Calc´ ulese Q b. 3. Resu´elvase Rx Q b (z o sustituci´on hacia atr´as)
Estos dos algoritmos pueden presentar problemas si la matriz A es singular o muy pr´oxima a serlo; i.e., κpAq es muy grande. Si la situaci´on es ´esta, tenemos una alternativa: los valores singulares.
168
El Problema de M´ınimos Cuadrados
8.3.1.
El problema de m´ınimos cuadrados y la inversa MoorePenrose
La inversa generalizada de Moore-Penrose se puede definir como la matriz que soluciona el problema de m´ınimos cuadrados. Es decir, de la misma que la soluci´on del sistema Ax b es x A1 b siempre que A sea invertible, la soluci´on del problema de m´ınimos cuadrados es x0 A: b. Veremos que esto es, en efecto as´ı, y que este hecho nos proporciona un algoritmo alternativo para resolver el problema de m´ınimos cuadrados. Recordemos que x0 es soluci´on del problema se m´ınimos cuadrados si y s´olo si es soluci´on del sistema Ax PA b, siendo PA la proyecci´on ortogonal sobre Im A. Recordemos tambi´en que si r rang A y A U ΣV es una descomposici´on completa de A en valores singulares, entonces las r primeras columnas de U son una base ortonormal de Im A. (Proposici´on 3.10). Sea Ur la submatriz de U formada por sus r primeras columnas. Como son una base ortonormal de Im A, la matriz Ur Ur es la proyecci´on ortogonal sobre Im A. Es decir, PA Ur Ur . As´ı pues, si ponemos Σ
Σr 0 0 0
entonces el sistema Ax PA b es equivalente a
Ur Σr 0 V x Ur Ur b. Y este sistema es equivalente a Σr 0 V x Ur b dado que Ur Ur c Ur b y y V x. Entonces,
Ax PA b ô Σr 0 y Si Σr
Ir . Pongamos
c.
Diagpσ1, . . . , σr q, la soluci´on general de este sistema es y
c1{σ1
c2 {σ2
cr {σr yr
1
yn
T
con yr 1 ,. . . , yn , n´ umeros arbitrarios. Si el rango de A es completo, n, la soluci´on del sistema queda completamente determinada; cosa que ya hab´ıamos demostrado. Pero si rang A n entonces hay infinitas soluciones del problema de m´ınimos cuadrados. Entre ellas, se suele escoger la soluci´on de norma m´ınima. Y ´esta es la que se consigue haciendo yr 1 yn 0. Finalmente, como y V x, tenemos que x V y.
8.3 Algoritmos para calcular la soluci´on del problema de m´ınimos cuadrados
169
Adem´as, }x}2 }V y }2 }y }2 . En definitiva, el vector x0 de norma m´ınima que soluciona el problema de m´ınimos cuadrados es:
c1 {σ1 c2 {σ2 .. .
c1 {σ1 c {σ 2 2 x0 V y0 V cr {σr Vr .. , . 0 . cr {σr .. 0 donde Vr es la submatriz de V formada por sus r primeras columnas. Ahora vamos a “volver hacia atr´as” a fin de recuperar la soluci´on del problema en t´erminos de los datos originales: A y b. Recordemos, antes de nada, que con las notaciones introducidas A Ur Σr Vr es una descomposici´on reducida de A en valores singulares. En consecuencia A:
Vr Σr 1Ur
es la inversa generalizada o Moore-Penrose de A. Ahora bien,
c1 {σ1 c {σ 2 2 1 1 x0 Vr .. Vr Σ r c Vr Σr Ur b, . cr {σr
porque c Ur b. As´ı pues, el vector de norma m´ınima que soluciona el problema de m´ınimos cuadrados es x0
Vr Σr 1Urb A:b
tal y como hab´ıamos anunciado. Este proceso, adem´as, nos proporciona un algoritmo para calcular la soluci´on del problema de m´ınimos cuadrados en el caso m´as general.
170
El Problema de M´ınimos Cuadrados
Algoritmo m´ınimos cuadrados via valores singulares Dada A P Fmn y dado b P Fm 1. Calcular la descomposici´on reducida de A en valores singulares: A U ΣV , U P Fmr , V P Fnr y Σ Diagpσ1 , . . . , σr q. 2. Calcular c U b. 3. Calcular x V Σ1 c. Este es el vector de menor norma que soluciona el problema de m´ınimos cuadrados.
8.4.
El condicionamiento del problema de m´ınimos cuadrados
El condicionamiento del problema de m´ınimos cuadrados es un tema importante porque tiene implicaciones no triviales en el estudio de la estabilidad de los algoritmos para este problema. Como es habitual analizaremos en detalle el condicionamiento del problema y estudiaremos la estabilidad de los algoritmos mediante ejemplos significativos. Recordemos que el problema de m´ınimos cuadrados consiste en lo siguiente (supondremos en lo sucesivo que la matriz del sistema tiene rango completo): Dada A P Fmn de rango completo y m ¥ n, y dado b P Fm , calcular x P Fn para que }Ax b}2 sea m´ınima.
(8.2)
Ya sabemos que, en este caso, la soluci´on del problema es u ´nica y viene dada por x A: b donde A: P Fnn es la pseudoinversa o inversa generalizada de Moore-Penrose de A (ver 3.6). Adem´as, si y Ax es el vector en Im A m´as pr´oximo a b entonces y
P b,
8.4 El condicionamiento del problema de m´ınimos cuadrados
171
siendo P la proyecci´on ortogonal sobre Im A. Nuestro objetivo es estudiar el condicionamiento del Problema (8.2) respecto a perturbacione en A y en b. Es decir, los datos del problema son la matriz A y el vector b. La soluci´on es el vector de coeficientes x o el correspondiente vector y Ax. As´ı pues Datos: A, b. Soluciones: x, y. Tenemos as´ı, en realidad, cuatro posibles cuestiones de condicionamiento: Error relativo en x o y respecto a peque˜ nas perturbaciones en los datos b o A. El objetivo en esta secci´on es dar un resultado fundamental sobre el condicionamiento del problema de m´ınimos cuadrados. Para entender el enunciado y como apoyo a la demostraci´on necesitamos introducir algunos conceptos y un par de resultados auxiliares. En primer lugar debemos recordar que para una matriz A P Fmn de rango completo n el n´ umero de condici´on es κ2 pAq }A}2 }A: }2
b θ
σσ1 . n
r=Ax−b y=Ax=Pb
Im A
Figura 8.4: El problema de m´ınimos cuadrados. En segundo lugar, necesitamos el conceto de ´angulo entre dos vectores de Fm . Lo definimos, como es habitual, a partir de la desigualdad de Cauchy-Schwartz: si x, y P Fm entonces |xy| ¤ }x}2}y}2,
172
El Problema de M´ınimos Cuadrados
de modo que
1 ¤ }x}x }yy} ¤ 1. 2
2
Se define, entonces el a´ngulo, θ, de x, y como θ
arc cos }x}x }yy} 2
Equivalentemente cos θ
}x}x }yy} 2
.
2
. 2
En nuestro caso, para calcular el ´angulo de b e y debemos hacer el producto escalar b y b P b. Teniendo en cuenta que P es una proyecci´on (P 2 P ) ortogonal pP P ) resulta que b y
bP b bP P b bP P b }P b}22 }y}22.
As´ı pues
}}yb}}2 . 2 Para el seno usamos la identidad sen2 θ 1 cos2 θ: }y}22 }b}22 }y}22 . sen2 θ 1 }b}22 }b}22 Ahora bien, b y b y siendo y y b y ortogonales (y P Im A y b y b P b pIm P qb P pIm AqK ). Por consiguiente, usando el Teorema de Pit´agoras, }b}22 }y}22 }b y}22 y }b y }2 sen θ }b} cos θ
2
Todas estos resultados generalizan los habituales en R2 (Figura 8.4) En tercer lugar, para cualquier norma inducida }Ax} ¤ }A} }x}. Necesitaremos una medidad de lo lejos o cerca que est´a }y } }Ax} de su m´aximo valor posible: η
2 }x}2 }A}}y2 }}x}2 }A}}Ax } . 2
2
Estos par´ametros tienen el siguiente rango de variaci´on: 1 ¤ κpAq 8 0 ¤ θ
¤ π {2
1¤η
¤ κpAq
8.4 El condicionamiento del problema de m´ınimos cuadrados
173
Todas estas acotaciones son claras excepto, quiz´a, la u ´ltima que viene de la siguiente observaci´on: }x}2 }A1Ax}2 ¤ }A1}2}Ax}2, de modo que
} A}2 }A1 }2 }Ax}2 κ2pAq. η¤ }Ax}2
El resultado previo que necesitamos es el siguiente Lema 8.4 Para A P Fmn de rango completo n se tiene:
}pAAq1A}2 σ1 .
(8.3)
}pAAq1}2 σ12
(8.4)
n
siendo σ1
¥ ¥ σn los valores singulares de A.
n
Demostraci´ on.- La propiedad (8.4) es una consecuencia inmediata de que los valores singulares de A son las ra´ıces cuadradas positivas de los valores propios de A A, y de que 1{σn2 es el mayor valor singular de pA Aq1 (ver Porposiciones 3.12 y 3.16). Para probar la propiedad (8.3) se usa el Teorema SVD para ver que los valores singulares de pA Aq1 A y los de A: coinciden. Podemos ahora plantear y demostrar el resultado principal Teorema 8.5 Sean b P Fm y A P Fmn con b 0 y rang A n ¤ m. Para el problema de m´ınimos cuadrados (8.2) se cumplen las siguientes propiedades, todo respecto de la norma `2 : (i) El n´ umero de condici´on del problema de calcular y 1 . cos θ
Ax P Fn respecto de b es (8.5)
174
El Problema de M´ınimos Cuadrados
(ii) El n´ umero de condici´on del problema de calcular x P Fn respecto de b es κ2 pAq . η cos θ
(8.6)
(iii) Sea x˜ la soluci´on del problema de m´ınimos cuadrados relativo a minimizar }2 σn entonces }pA δAqx˜ b}2. Si y˜ pA δAqx˜ y }}δA A}2 σ1
}y˜ y}2 ¤ κ2pAq }y}2 cos θ
Op2 q.
(8.7)
(iv) En las mismas condiciones del apartado anterior
}x˜ x}2 ¤ κ pAq 2 }x}2
κ2 pAq2 tan θ η
Op2 q.
(8.8)
(v) Sea x˜ es la soluci´on del problema de m´ınimos cuadrados relativo a minimizar }pA δAqx˜ pb δbq}2. Si sen θ 1 y m´ax
"
}δA}2 , }δb}2 * σn }A}2 }b}2 σ1
¥ ¥ σn los valores singulares de A, entonces
* }x˜ x}2 ¤ "κ pAq 1 κ2 pAq2 tan θ 1 2 }x}2 η cos θ η
siendo σ1
Op2 q.
(8.9)
En el caso especial en que m n, el problema de m´ınimos cuadrados se reduce al de la resoluci´on de un sistema de ecuaciones lineales. En tal caso θ 0 y las cotas de (8.6) y (8.8) se reducen a κ2 pAq{η y κ2 pAq recuperando los resultados vistos en la Lecci´on 4. Demostraci´ on.- Demostraremos las propiedades (i) y (v). La propiedad (ii) se demuestra de manera muy similar a la propirdad (i), la propiedad (iv) es una caso particular de la (v) y la propiedad (iii) requiere una demostraci´on independiente pero parecida a la (iv). (i) La relaci´on entre b e y es y
Pb
8.4 El condicionamiento del problema de m´ınimos cuadrados
175
siendo P la proyeccci´on ortogonal sobre Im A. Vimos en la Secci´on 4.3 del Cap´ıtulo 4 que el n´ umero de condici´on del problema de multiplicar Ax para A dada y x variable: f : Fn Ñ Fm x ; b Ax es κ
}f 1pxq} }A} }x} , }f pxq}{}x} }Ax}
(8.10)
En nuestro caso se trata del n´ umero de condici´on del problema y b para P dado. Entonces }P }2}b}2 κpbq }y }2
P b respecto de
Como P es una proyecci´on ortogonal P QQ para alguna matriz Q P Fnm con r una matriz unitaria. Entonces columnas ortonormales. Sea U Q Q U P U
Q I 0 r QQ Q Q n r 0 0 Q
Por lo tanto, los valores singulares de P son todos iguales a 1 y en particular
}P }2 σ1pP q 1 }y}2 . En conclusi´on: Finalmente, recordemos que cos θ }b}2 }P }2}b}2 1 , κpbq }y}2 cos θ tal y como se deseaba demostrar.
}2 σn 1 δA y f 1 δb. Por hip´otesis, }}δA A}2 σ1 y como }A}2 σ1 , tenemos que }δA}2 σn . Por el Teorema 3.18 resulta que rangpA δAq n y consecuentemente rangpA tE q n para todo t P r0, s. Se (v) En primer lugar definimos E
sigue entonces que el sistema
pA
tE q pA
tE qxptq pA
tE q pb
tf q
(8.11)
176
El Problema de M´ınimos Cuadrados
tiene una u ´nica soluci´on para cada t P r0, s. Invertir una matriz es una funci´on diferenciable en los elementos de la matriz. Por lo tanto x es una funci´on diferenciable de t en r0, s. Por la definici´on de diferencial x1 p0q
xpq xp0q
Op2 q.
Como xpq x˜, xp0q x, b 0 y sen θ
1, tenemos que x 0 y }x˜ x}2 }x1p0q}2 Op2q. (8.12) }x}2 }x}2 Necesitamos una estimaci´on de }x1 p0q}2 . Para ello calculamos x1 p0q en (8.11). Primero derivamos E pA tE qxptq
tE q Exptq
pA
pA
tE q pA tE qx1 ptq E ppb tf q
tE q f,
pA
y sustitu´ımos t 0 recordando que xp0q x: E Ax As´ı
A Ex
x1 p0q pA Aq1 A pf
A Ax1 p0q A f
E b.
Exq pAAq1E pb Axq.
Tomando normas:
}x1p0q}2 ¤ }pAAq1A}2p}f }2 }E }2}x}2q }pAAq1}2}E }2}b Ax}2. Por una parte, }E }2
1 }δA}2 ¤ }A}2. Y de la misma forma }f }2 ¤ }b}2. Tambi´en
}A}2 }A}2, as´ı que
}x1p0q}2 ¤ }pAAq1A}2p}A}2}x}2 }b}2q }pAAq1}2}A}2}b Ax}2 ¤ }2 . ¤ }pAAq1A}2}A}2 }}Ab}}2 }x}2 }pAAq1}2}A}22 }b }AAx } 2
Por otra parte, por el Lema 8.4, κ2 pAq2 . Entonces
2
}pAAq1A}2}A}2 κ2pAq y }pAAq1}2}A}22
}x1p0q}2 ¤ κ pAq }b}2{}Ax}2 2 }x}2 }A}2}x}2{}Ax}2
1
κ2 pAq2
}b Ax}2{}Ax}2 . }A}2}x}2{}Ax}2
8.5 Estabilidad de los algoritmos para el problema de m´ınimos cuadrados
Ahora bien, η
}2}x}2 , cos θ }Ax}2 }A}Ax } }b} 2
y tan θ
2
177
}2 . En consecuencia }b}AxAx } 2
}x1p0q}2 ¤ κ pAq 1 2 }x}2 η cos θ
1
κ2 pAq2
tan θ η
sustituyendo en (8.12)
}x˜ x}2 ¤ "κ pAq 1 2 }x}2 η cos θ
1
κ2 pAq2 tan θ η
*
Op2 q,
tal y como se quer´ıa demostrar.
8.5.
Estabilidad de los algoritmos para el problema de m´ınimos cuadrados
Un estudio experimental de la estabilidad de los algoritmos para la resoluci´on del problema de m´ınimos cuadrados lineal es el objetivo de una de las pr´acticas obligatorias con MATLAB que habr´a de realizarse en este curso. En ella se plantea un problema de m´ınimos cuadrados cuya soluci´on debe ser calculada con los tres algoritmos vistos en la Secci´on 8.3. A su vez, para los algoritmos basados en la factorizaci´on QR, se utilizar´an los tres algoritmos vistos en las Lecciones 6 y 7; es decir, los algoritmos cl´asico y modificado de Gram-Schmidt y el algoritmo de Householder. A trav´es de dicho experimento se comprobar´a que los algoritmos basados en la resoluci´on de las ecuaciones normales y el factorizaci´on de Cholesky as´ı como los basados en la factorizaci´on QR obtenida con los algoritmos de Gram-Schmidt pueden dar resultados con mucho error si la matriz del sistema est´a mal condicionada. Por lo tanto, estos algoritmos son inestables. No obstante, hay resultados (ver [11, Sec 20.4]) que muestran que para matrices bien condicionadas el algoritmo basado en la resoluci´on de las ecuaciones normales es estable hacia atr´as y para valores de m mucho mayores que n, el coste operacional sensiblemente menor que el m´etodo basado en la factorizaci´on QR mediante reflexiones de Householder. En tale situaciones, el m´etodo basado en la resoluci´on de las ecuaciones normales ser´ıa el preferido. Por otra parte, a pesar de que el m´etodo basado en la factorizaci´on QR mediante el algoritmo modificado de Gram-Schmidt es inestable cuando se implementa de la manera indicada en la Secci´on 8.3, hay una forma de hacerlo que lo hace estable
178
El Problema de M´ınimos Cuadrados
hacia atr´as (ver [11, Sec 20.3]). Finalmente, los algoritmos basados en la factorizaci´on QR mediante reflexiones de Householder y valores singulares son estables hacia atr´as para sistemas en los que A es una matriz de rango completo. A continuaci´on se enuncian los teoremas que lo hace expl´ıto. Teorema 8.6 Supongamos que el problema de m´ınimos cuadrados con una matriz A de rango completo se resuelve mediante el algoritmo QR por reflexiones de Householder o mediante la descomposici´on en valores singulares en un ordenador que cumple los axiomas (5.2) y (5.4) de la Lecci´on 5. Este algoritmo es estable hacia atr´as: para cada A P Fmn (m ¥ n y rang A n) existe una perturbaci´on δA tal que }δA} Op q M }A} p producida por el algoritmo satisface y la soluci´on x
}pA
p b}2 δAqx
m´ ın }pA xPF n
δAqx b}2 .
Para terminar, conviene mencionar que si la matriz A no es de rango completo, sabemos que la soluci´on no est´a determinada de forma u ´nica y el u ´nico algoritmo completamente estable es el basado en la descomposici´on de A en valores singulares.