4.1: Mínimos Cuadrados
( \newcommand{\kernel}{\mathrm{null}\,}\)
Introducción
Aprendimos en el capítulo anterior que no esAx=b necesario poseer una solución cuando el número de filas deA excede su rango, es decir,r<m. Como esta situación surge con bastante frecuencia en la práctica, típicamente bajo la apariencia de 'más ecuaciones que incógnitas', establecemos una justificación para lo absurdoAx=b.
Las ecuaciones normales
El objetivo es elegirx tal queAx esté lo más cerca posible deb. Midiendo la cercanía en términos de la suma de los cuadrados de los componentes llegamos al problema de 'mínimos cuadrados' de minimizar
res
(||Ax−b||)2=(Ax−b)T(Ax−b)
sobre todox∈R. El camino hacia la solución es iluminado por el Teorema Fundamental. Más precisamente, escribimos
∀bR,bN,bR∈R(A)∧bN∈N(AT):(b=bR+bN). Al señalar que (i)∀bR,x∈Rn:((Ax−bR)∈R(A)) y (ii)R(A)⊥N(AT) llegamos al Teorema de Pitágoras.
norm2(Ax−b)=(||Ax−bR+bN||)2
=(||Ax−bR||)2+(||bN||)2
Ahora queda claro del Teorema de Pitágoras que el mejorx es el que satisface
Ax=bR
ComobR∈R(A) esta ecuación en efecto posee una solución. Sin embargo, todavía tenemos que especificar cómo se calculabR dadob. Si bien una expresión explícita para la proyecciónbR ortogonal deb sobreR(A), en términos deA yb está a nuestro alcance, estrictamente hablando, no la requeriremos. Para ver esto, notemos que six satisface la ecuación anterior entonces la proyección ortogonal deb sobreR(A), en términos deA yb está a nuestro alcance, estrictamente hablando, no la requeriremos. Para ver esto, notemos que six satisface la ecuación anterior entonces
Ax−b=Ax−bR+bN
=−bN
Como nobN se calcula más fácilmente debR lo que puede afirmar que sólo vamos en círculos. Sin embargo, la información 'práctica' en la ecuación anterior es que(Ax−b)∈AT, es decirAT(Ax−b)=0, es decir,
ATAx=ATb
ComoATb∈R(AT) independientemente deb este sistema, a menudo referido como las ecuaciones normales, de hecho tiene una solución. Esta solución es única siempre y cuando las columnas deATA sean linealmente independientes, es decir, siempre y cuandoN(ATA)=0. Recordando el Capítulo 2, Ejercicio 2, observamos que esto equivale aN(A)={0}
El conjunto dex∈bR para el que el inadaptado(||Ax−b||)2 es más pequeño está compuesto por aquellosx para los que SiempreATAx=ATb hay al menos uno de esosx. Hay exactamente uno de esosx siN(A)={0}.
Como ejemplo concreto, supongamos con referencia a la Figura 1 queA=(110100) yA=(111)

Como nob≠R(A) hayx tal queAx=b. En efecto(||Ax−b||)2=(x1+x2+−1)2+(x2−1)2+1≥1,, con el mínimo obtenido de manera única enx=(01), de acuerdo con la solución única de la ecuación anterior, paraATA=(1112) yATb=(12). Ahora reconocemos, a posteriori, quebR=Ax=(110) es la proyección ortogonal de b sobre el espacio de columna deA.
Aplicación de mínimos cuadrados al problema de la prueba biaxial
Formularemos la identificación de las 20 rigideces de fibra en esta figura anterior, como problema de mínimos cuadrados. Visualizamos la carga, los 9 nodos y midiendo los 18 desplazamientos asociados,x. A partir del conocimientox yf deseamos inferir los componentes deK=diag(k) dondek se encuentra el vector de rigideces de fibra desconocidas. El primer paso es reconocer que
ATKAx=f
puede escribirse como
∀B,B=ATdiag(Ax):(Bk=f)
Aunque conceptualmente simple esto no es de gran utilidad en la práctica, puesB es de 18 por 20 y de ahí que la ecuación anterior posee muchas soluciones. La salida es computark como resultado de más de un experimento. Veremos que, para nuestra pequeña muestra, bastarán 2 experimentos. Para ser precisos, suponemos quex1 es el desplazamiento producido por la cargaf1 mientras quex2 es el desplazamiento producido por la cargaf2. Luego recogimos las piezas asociadas en
B=(ATdiag(Ax1)ATdiag(Ax2))
y
f=(f1f2).
EstoB es 36-por-20 y así el sistemaBk=f está sobredeterminado y por lo tanto maduro para mínimos cuadrados.
Se procede entonces a ensamblarB yf. Suponemosf1 yf2 para corresponder al estiramiento horizontal y vertical
f1=(−100010−100010−100010)T
f2=(0101010101010−10−10−1)T
respectivamente. A los efectos de nuestro ejemplo suponemos que cada unokj=1 exceptok8=5. EnsamblamosATKA como en el Capítulo 2 y resolvemos
ATKAxj=fj
con la ayuda de la pseudoinversa. Para impartir cierta `realidad' a este problema, contaminamos cada unoxj con un 10 por ciento de ruido antes de construirB
BTBk=BTf
observamos que Matlab resuelve este sistema cuando se presenta con K=b\ f
cuando BB es rectangular. Hemos trazado los resultados de este procedimiento en el enlace. La fibra rígida se identifica fácilmente.

Proyecciones
Desde un punto de vista algebraico La ecuación es una elegante reformulación del problema de mínimos cuadrados. Aunque fácil de recordar, desgraciadamente oscurece el contenido geométrico, sugerido por la palabra 'proyección', de Ecuación. Como las proyecciones surgen frecuentemente en muchas aplicaciones hacemos una pausa aquí para desarrollarlas con más cuidado. Con respecto a las ecuaciones normales observamos que siN(A)={0} entonces
x=(ATA)−1ATb
y así la proyección ortogonal de bb sobreR(A) es:
bR=Ax
=A(ATA)−1ATb
Definiendo
P=A(ATA)−1AT
toma la formabR=Pb. De acuerdo con nuestra noción de lo que debería ser una 'proyección', esperamos queP mapee vectores que no estén dentroR(A) mientras dejan vectores yaR(A) ilesos.R(A) De manera más sucinta, esperamos quePbR=bR i.ePbR=PbR.,. Como este último debe sostenerse para todosb∈Rm esperamos que
P2=P
Encontramos que de hecho
P2=A(ATA)−1ATA(ATA)−1AT
=A(ATA)−1AT
=P
También observamos que elP es simétrico. Dignificamos estas propiedades a través de
Una matrizP que satisfaceP2=P se llama proyección. Una proyección simétrica se denomina proyección ortogonal.
Nos hemos esforzado por motivar el uso de la palabra 'proyección'. Quizás se esté preguntando sin embargo qué tiene que ver la simetría con la ortogonalidad. Esto lo explicamos en términos de la tautología
b=Pb−Ib
Ahora bien, siP es una proyección entonces también lo esI−P. Además, siP es simétrico entonces el punto producto deb.
\ [\ begin {align*} (Pb) ^T (I-P) b &= b^ {T} P^ {T} (I-P) b\\ [4pt] &= b^ {T} (P-P^ {2}) b\\ [4pt] &= b^ {T} 0 b\\ [4pt] &= 0\ end {align*}
es decir,Pb es ortogonal a(I−P)b. Como ejemplos de proyecciones no ortogonales ofrecemos
P=(100−1200−14−121)
yI−P. Por último, señalemos que la fórmula centralP=A(ATA)−1AT, es incluso un poco más general de lo que se anuncia. Se ha facturado como la proyección ortogonal sobre el espacio de columna deA. Sin embargo, a menudo surge la necesidad de la proyección ortogonal sobre algún subespacio arbitrario M. La clave para usar el PP antiguo es simplemente darse cuenta de que cada subespacio es el espacio de columna de alguna matriz. Más precisamente, si
{x1,⋯,xm}
es una base para MM entonces claramente si estosxj se colocan en las columnas de una matriz llamadaA entoncesR(A)=M. Por ejemplo, siM es la línea a través de(11)T entonces
P=(11)12(11)
P=12(1111)
es proyección ortogonal sobreM.
Ejercicios
Gilbert Strang se estiró sobre una rejilla a longitudesl=6,7,8 pies bajo fuerzas aplicadas def=1,2,4 toneladas. Asumiendo la ley de Hookel−L=cf, encuentra su cumplimientoc,, y altura originalL,, por mínimos cuadrados.
Con respecto al ejemplo del § 3 señalar que, debido a la generación aleatoria del ruido que mancha los desplazamientos, se obtiene una 'respuesta' diferente cada vez que se invoca el código.
- Escriba un bucle que invoque el código un número de veces estadísticamente significativo y envíe gráficas de barras de la rigidez promedio de la fibra y su desviación estándar para cada fibra, junto con el archivo M asociado.
- Experimentar con diversos niveles de ruido con el objetivo de determinar el nivel por encima del cual se hace difícil discernir la fibra rígida. Explique cuidadosamente sus hallazgos.
Encuentra la matriz queR3 proyecta en la línea abarcada por(101)T.
Encuentra la matriz queR3 proyecta en la línea abarcada por(101)T y(11−1)T.
SiP es laRm proyección de un subespacio k—dimensionalM, ¿cuál es el rangoP y cuál esR(P)?