17.3: Mínimos Cuadrados

Medidas de Cercanía

En el apartado anterior, observamos que en general no es posible encontrar una solución a un sistema sobredeterminado. Nuestro objetivo es así encontrar\(z\) tal que\(B z\) esté “cerca” de\(g\), es decir,\(z\) tal que\[B z \approx g,\] para\(B \in \mathbb{R}^{m \times n}, m>n\). Por conveniencia, introduzcamos la función residual, que se define como\[r(z) \equiv g-B z .\] Tenga en cuenta que\[r_{i}=g_{i}-(B z)_{i}=g_{i}-\sum_{j=1}^{n} B_{i j} z_{j}, \quad i=1, \ldots, m .\] Así,\(r_{i}\) es la “extensión” a la que\(i\) -ésima ecuación no\((B z)_{i}=g_{i}\) se satisface. En particular, si\(r_{i}(z)=0\)\(i=1, \ldots, m\), entonces\(B z=g\) y\(z\) es la solución al sistema lineal. Observamos que el residual es una medida de cercanía descrita por\(m\) valores. Es más conveniente tener un único valor escalar para evaluar el grado en que se satisface la ecuación. Una forma sencilla de lograrlo es tomar una norma del vector residual. Diferentes normas dan como resultado diferentes medidas de cercanía, que a su vez producen diferentes soluciones de mejor ajuste.

Comenzar material avanzado

Consideremos los dos primeros ejemplos, ninguno de los cuales vamos a perseguir en este capítulo.

Ejemplo 17.3.1\(\ell_{1}\) minimización

El primer método se basa en medir el residuo en la norma 1. El escalar que representa la extensión del desajuste es\[J_{1}(z) \equiv\|r(z)\|_{1}=\sum_{i=1}^{m}\left|r_{i}(z)\right|=\sum_{i=1}^{m}\left|(g-B z)_{i}\right| .\] El mejor\(z\), denotado por\(z^{*}\), es el\(z\) que minimiza la extensión del desajuste medido en\(J_{1}(z)\), es decir,\[z^{*}=\arg \min _{z \in \mathbb{R}^{m}} J_{1}(z) .\] El\(\arg \min _{z \in \mathbb{R}^{n}} J_{1}(z)\) devuelve el argumento \(z\)que minimiza la función\(J_{1}(z)\). En otras palabras,\(z^{*}\) satisface\[J_{1}\left(z^{*}\right) \leq J_{1}(z), \quad \forall z \in \mathbb{R}^{m}\] Este problema de minimización puede formularse como un problema de programación lineal. El minimizador no es necesariamente único y el procedimiento de solución no es tan sencillo como el resultante de la norma 2. Por lo tanto, aquí no vamos a perseguir esta opción.

Ejemplo 17.3.2\(\ell_{\infty}\) minimización

El segundo método se basa en medir el residuo en la\(\infty\) norma. El escalar que representa la extensión del desajuste es\[J_{\infty}(z) \equiv\|r(z)\|_{\infty}=\max _{i=1, \ldots, m}\left|r_{i}(z)\right|=\max _{i=1, \ldots, m}\left|(g-B z)_{i}\right| .\] Lo mejor\(z\) que minimiza\(J_{\infty}(z)\) es\[z^{*}=\arg \min _{z \in \mathbb{R}^{n}} J_{\infty}(z) .\] Este llamado problema min-max también se puede lanzar como un problema de programación lineal. Nuevamente, este procedimiento es bastante complicado, y la solución no es necesariamente única.

Material Avanzado

Formulación de mínimos cuadrados (\(\ell_{2}\)minimización)

Minimizar el residuo medido en (digamos) la 1-norma o\(\infty\) -norma da como resultado un problema de programación lineal que no es tan fácil de resolver. Aquí mostraremos que medir el residuo en la norma 2 da como resultado un problema de minimización particularmente simple. Además, la solución al problema de minimización es única asumiendo que la matriz\(B\) es de rango completo - tiene columnas\(n\) independientes. Supondremos que efectivamente\(B\) tiene columnas independientes. La función escalar que representa el grado de desajuste para la\(\ell_{2}\) minimización es\[J_{2}(z) \equiv\|r(z)\|_{2}^{2}=r^{\mathrm{T}}(z) r(z)=(g-B z)^{\mathrm{T}}(g-B z) .\] Tenga en cuenta que consideramos el cuadrado de la norma 2 por conveniencia, en lugar de la norma 2 en sí. Nuestro objetivo es encontrar\(z^{*}\) tal aquello\[z^{*}=\arg \min _{z \in \mathbb{R}^{n}} J_{2}(z),\] que equivale a encontrar\(z^{*}\) con\[\left\|g-B z^{*}\right\|_{2}^{2}=J_{2}\left(z^{*}\right)<J_{2}(z)=\|g-B z\|_{2}^{2}, \quad \forall z \neq z^{*} .\] (Nota “arg min” se refiere al argumento que minimiza: entonces “min” es el mínimo y “arg min” es el minimizador.) Tenga en cuenta que podemos escribir nuestra función objetiva\(J_{2}(z)\) como\[J_{2}(z)=\|r(z)\|_{2}^{2}=r^{\mathrm{T}}(z) r(z)=\sum_{i=1}^{m}\left(r_{i}(z)\right)^{2} .\] En otras palabras, nuestro objetivo es minimizar la suma del cuadrado de los residuos, es decir, mínimos cuadrados. Así, decimos que\(z^{*}\) es la solución de mínimos cuadrados al sistema sobredeterminado\(B z=g: z^{*}\) es la\(z\) que hace\(J_{2}(z)\) -la suma de los cuadrados de los residuales- lo más pequeña posible.

Tenga en cuenta que si\(B z=g\) tiene una solución, la solución de mínimos cuadrados es la solución al sistema sobredeterminado. Si\(z\) es la solución, entonces\(r=B z-g=0\) y en particular\(J_{2}(z)=0\), que es el valor mínimo que\(J_{2}\) puede tomar. Así, la solución\(z\) es el minimizador de\(J_{2}: z=z^{*}\). Derivamos ahora un procedimiento para resolver el problema de mínimos cuadrados para un caso más general donde\(B z=g\) no tiene solución.

Por conveniencia, bajamos el subíndice 2 de la función objetiva\(J_{2}\), y simplemente lo denotamos por\(J\). Nuevamente, nuestro objetivo es encontrar\(z^{*}\) tal que\[J\left(z^{*}\right)<J(z), \quad \forall z \neq z^{*} .\] Ampliando la expresión para\(J(z)\), tenemos\[\begin{aligned} J(z) &=(g-B z)^{\mathrm{T}}(g-B z)=\left(g^{\mathrm{T}}-(B z)^{\mathrm{T}}\right)(g-B z) \\ &=g^{\mathrm{T}}(g-B z)-(B z)^{\mathrm{T}}(g-B z) \\ &=g^{\mathrm{T}} g-g^{\mathrm{T}} B z-(B z)^{\mathrm{T}} g+(B z)^{\mathrm{T}}(B z) \\ &=g^{\mathrm{T}} g-g^{\mathrm{T}} B z-z^{\mathrm{T}} B^{\mathrm{T}} g+z^{\mathrm{T}} B^{\mathrm{T}} B z \end{aligned}\] donde hemos utilizado la regla de transposición que nos dice eso\((B z)^{\mathrm{T}}=z^{\mathrm{T}} B^{\mathrm{T}}\). Observamos que\(g^{\mathrm{T}} B z\) es un escalar, por lo que no cambia bajo la operación de transposición. Así, se\(g^{\mathrm{T}} B z\) puede expresar como de\[g^{\mathrm{T}} B z=\left(g^{\mathrm{T}} B z\right)^{\mathrm{T}}=z^{\mathrm{T}} B^{\mathrm{T}} g,\] nuevo por la regla de transposición. La función simplifica\(J\) así a\[J(z)=g^{\mathrm{T}} g-2 z^{\mathrm{T}} B^{\mathrm{T}} g+z^{\mathrm{T}} B^{\mathrm{T}} B z .\] Por conveniencia, definamos\(N \equiv B^{\mathrm{T}} B \in \mathbb{R}^{n \times n}\), para que\[J(z)=g^{\mathrm{T}} g-2 z^{\mathrm{T}} B^{\mathrm{T}} g+z^{\mathrm{T}} N z .\] sea sencillo confirmar que cada término en la expresión anterior es efectivamente un escalar.

La solución al problema de minimización viene dada por\[N z^{*}=d,\] dónde\(d=B^{\mathrm{T}} g\). A la ecuación se le llama la ecuación “normal”, que puede escribirse como\[\left(\begin{array}{cccc} N_{11} & N_{12} & \cdots & N_{1 n} \\ N_{21} & N_{22} & \cdots & N_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ N_{n 1} & N_{n 2} & \cdots & N_{n n} \end{array}\right)\left(\begin{array}{c} z_{1}^{*} \\ z_{2}^{*} \\ \vdots \\ z_{n}^{*} \end{array}\right)=\left(\begin{array}{c} d_{1} \\ d_{2} \\ \vdots \\ d_{n} \end{array}\right) .\] La existencia y singularidad de\(z^{*}\) está garantizada asumiendo que las columnas de\(B\) son independientes.

Proporcionamos a continuación la prueba de que\(z^{*}\) es el minimizador único de\(J(z)\) en el caso en el que\(B\) tiene columnas independientes.

Prueba

Primero mostramos que la matriz normal\(N\) es simétrica positiva definida, es decir,\[x^{\mathrm{T}} N x>0, \quad \forall x \in \mathbb{R}^{n}(x \neq 0),\] asumiendo que las columnas de\(B\) son linealmente independientes. La matriz normal\(N=B^{\mathrm{T}} B\) es simétrica porque\[N^{\mathrm{T}}=\left(B^{\mathrm{T}} B\right)^{\mathrm{T}}=B^{\mathrm{T}}\left(B^{\mathrm{T}}\right)^{\mathrm{T}}=B^{\mathrm{T}} B=N .\] Para mostrar\(N\) es positivo definido, primero observamos\[x^{\mathrm{T}} N x=x^{\mathrm{T}} B^{\mathrm{T}} B x=(B x)^{\mathrm{T}}(B x)=\|B x\|^{2} .\] que es decir,\(x^{\mathrm{T}} N x\) es la 2-norma de\(B x\). Recordemos que la norma de un vector es cero si y sólo si el vector es el vector cero. En nuestro caso,\[x^{\mathrm{T}} N x=0 \quad \text { if and only if } \quad B x=0 .\] Porque las columnas de\(B\) son linealmente independientes,\(B x=0\) si y solo si\(x=0\). Así, tenemos\[x^{\mathrm{T}} N x=\|B x\|^{2}>0, \quad x \neq 0 .\] Así,\(N\) es simétrico positivo definido.

Ahora recordemos que la función a minimizar es\[J(z)=g^{\mathrm{T}} g-2 z^{\mathrm{T}} B^{\mathrm{T}} g+z^{\mathrm{T}} N z .\] Si\(z^{*}\) minimiza la función, entonces para cualquiera\(\delta z \neq 0\), debemos tener\[J\left(z^{*}\right)<J\left(z^{*}+\delta z\right) \text {; }\] Vamos a expandir\(J\left(z^{*}+\delta z\right)\):\[\begin{aligned} J\left(z^{*}+\delta z\right) &=g^{\mathrm{T}} g-2\left(z^{*}+\delta z\right)^{\mathrm{T}} B^{\mathrm{T}} g+\left(z^{*}+\delta z\right)^{\mathrm{T}} N\left(z^{*}+\delta z\right), \\ &=\underbrace{g^{\mathrm{T}} g-2 z^{*} B^{\mathrm{T}} g+\left(z^{*}\right)^{\mathrm{T}} N z^{*}}_{J\left(z^{*}\right)}-2 \delta z^{\mathrm{T}} B^{\mathrm{T}} g+\delta z^{\mathrm{T}} N z^{*}+\underbrace{\left(z^{*}\right)^{\mathrm{T}} N \delta z}_{\delta z^{\mathrm{T}} N^{\mathrm{T}} z^{*}=\delta z^{\mathrm{T}} N z^{*}}+\delta z^{\mathrm{T}} N \delta z, \\ &=J\left(z^{*}\right)+2 \delta z^{\mathrm{T}}\left(N z^{*}-B^{\mathrm{T}} g\right)+\delta z^{\mathrm{T}} N \delta z . \end{aligned}\] Tenga en cuenta que\(N^{\mathrm{T}}=N\) porque\(N^{\mathrm{T}}=\left(B^{\mathrm{T}} B\right)^{\mathrm{T}}=B^{\mathrm{T}} B=N\). Si\(z^{*}\) satisface la ecuación normal,\(N z^{*}=B^{\mathrm{T}} g\), entonces\[N z^{*}-B^{\mathrm{T}} g=0\] y así\[J\left(z^{*}+\delta z\right)=J\left(z^{*}\right)+\delta z^{\mathrm{T}} N \delta z .\] El segundo término siempre es positivo porque\(N\) es positivo definido. Así, tenemos\[J\left(z^{*}+\delta z\right)>J\left(z^{*}\right), \quad \forall \delta z \neq 0,\] o, fijando\(\delta z=z-z^{*}\),\[J\left(z^{*}\right)<J(z), \quad \forall z \neq z^{*} .\] así,\(z^{*}\) satisfaciendo la ecuación normal\(N z^{*}=B^{\mathrm{T}} g\) es el minimizador de\(J\), es decir, la solución de mínimos cuadrados al sistema sobredeterminado\(B z=g\).

Ejemplo 17.3.3\ times 1\) mínimos cuadrados y su interpretación geométrica

Considerar un caso sencillo de un sistema sobredeterminado,\[B=\left(\begin{array}{c} 2 \\ 1 \end{array}\right)(z)=\left(\begin{array}{l} 1 \\ 2 \end{array}\right)\] Debido a que el sistema es\(2 \times 1\), hay un único parámetro escalar\(z\),, a elegir. Para obtener la ecuación normal, primero construimos la matriz\(N\) y el vector\(d\) (ambos son simplemente escalares para este problema):\[\begin{aligned} &N=B^{\mathrm{T}} B=\left(\begin{array}{ll} 2 & 1 \end{array}\right)\left(\begin{array}{l} 2 \\ 1 \end{array}\right)=5 \\ &d=B^{\mathrm{T}} g=\left(\begin{array}{ll} 2 & 1 \end{array}\right)\left(\begin{array}{l} 1 \\ 2 \end{array}\right)=4 . \end{aligned}\] Resolviendo la ecuación normal, obtenemos la solución de mínimos cuadrados\[N z^{*}=d \quad \Rightarrow \quad 5 z^{*}=4 \quad \Rightarrow \quad z^{*}=4 / 5 .\] Esta elección de\(z\) rendimientos\[B z^{*}=\left(\begin{array}{c} 2 \\ 1 \end{array}\right) \cdot \frac{4}{5}=\left(\begin{array}{c} 8 / 5 \\ 4 / 5 \end{array}\right),\] que por supuesto es diferente de\(g\).

El proceso se ilustra en la Figura 17.5. El lapso de la columna de\(B\) es la línea parametrizada por\(\left(\begin{array}{cc}2 & 1\end{array}\right)^{\mathrm{T}} z, z \in \mathbb{R}\). Recordemos que la solución\(B z^{*}\) es el punto en la línea que está más cerca\(g\) en el sentido de mínimos cuadrados, i.e.\[\left\|B z^{*}-g\right\|_{2}<\|B z-g\|, \quad \forall z \neq z^{*}\]

Figura 17.5: Ilustración de mínimos cuadrados en\(\mathbb{R}^{2}\).

Recordando que la\(\ell_{2}\) distancia es la distancia euclidiana habitual, esperamos que el punto más cercano sea la proyección ortogonal de\(g\) sobre la línea\(\operatorname{span}(\operatorname{col}(B))\). La cifra confirma que efectivamente así es el caso. Podemos verificar esto algebraicamente,\[B^{\mathrm{T}}\left(B z^{*}-g\right)=\left(\begin{array}{ll} 2 & 1 \end{array}\right)\left(\left(\begin{array}{c} 2 \\ 1 \end{array}\right) \cdot \frac{4}{5}-\left(\begin{array}{l} 1 \\ 2 \end{array}\right)\right)=\left(\begin{array}{ll} 2 & 1 \end{array}\right)\left(\begin{array}{c} 3 / 5 \\ -6 / 5 \end{array}\right)=0 .\] así, el vector residual\(B z^{*}-g\) y el espacio de columna de\(B\) son ortogonales entre sí. Si bien la ilustración geométrica de la ortogonalidad puede ser difícil para mínimos cuadrados de dimensiones superiores, la condición de ortogonalidad se puede verificar sistemáticamente utilizando el método algebraico.

Consideraciones computacionales

Analicemos el costo computacional de resolver el sistema de mínimos cuadrados. El primer paso es la formulación de la matriz normal,\[N=B^{\mathrm{T}} B,\] que requiere una multiplicación matriz-matriz de\(B^{\mathrm{T}} \in \mathbb{R}^{n \times m}\) y\(B \in \mathbb{R}^{m \times n}\). Debido a que\(N\) es simétrico, solo necesitamos calcular la parte triangular superior de\(N\), que corresponde a la realización de productos internos\(n(n+1) / 2 m\) -vector. Así, el costo computacional es\(m n(n+1)\). Formando el lado derecho,\[d=B^{\mathrm{T}} g,\] requiere una multiplicación matriz-vector de\(B^{\mathrm{T}} \in \mathbb{R}^{n \times m}\) y\(g \in \mathbb{R}^{m}\). Esto requiere\(n m\) -vector productos internos, por lo que el costo computacional es\(2 m n\). Este costo es insignificante en comparación con las\(m n(n+1)\) operaciones requeridas para formar la matriz normal. Por último, debemos resolver el sistema lineal\(n\) -dimensional\[N z=d\] Como veremos en la unidad de álgebra lineal, resolver el sistema lineal definido positivo\(n \times n\) simétrico requiere aproximadamente\(\frac{1}{3} n^{3}\) operaciones utilizando la factorización Cholesky (ya que discutir más a fondo en la Unidad V). Así, el recuento total de operaciones es\[C^{\text {normal }} \approx m n(n+1)+\frac{1}{3} n^{3} .\] Para un sistema derivado de la regresión,\(m \gg n\), por lo que podemos simplificar aún más la expresión a la\[C^{\text {normal }} \approx m n(n+1) \approx m n^{2},\] que es bastante modesta para\(n\) no demasiado grande.

Si bien el método basado en la ecuación normal funciona bien para sistemas pequeños, este proceso resulta ser numéricamente “inestable” para problemas más grandes. Visitaremos la noción de estabilidad más adelante; por ahora, podemos pensar en la estabilidad como una capacidad de un algoritmo para controlar la perturbación en la solución bajo una pequeña perturbación en los datos (o entrada). En general, nos gustaría que nuestro algoritmo fuera estable. A continuación se discute el método de elección.

Comenzar material avanzado

R\) Factorización y el Procedimiento Gram-Schmidt

Un procedimiento más estable para resolver el sistema sobredeterminado es el basado en la\(Q R\) factorización. \(Q R\)la factorización es un procedimiento para factorizar, o descomponer, una matriz\(B \in \mathbb{R}^{m \times n}\) en una matriz ortonormal\(Q \in \mathbb{R}^{m \times n}\) y una matriz triangular superior\(R \in \mathbb{R}^{n \times n}\) tal que\(B=Q R\). Una vez que tenemos tal factorización, podemos simplificar enormemente la ecuación normal\(B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}} g\). Sustitución de la factorización en la ecuación normal rinde\[B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}} g \Rightarrow R^{\mathrm{T}} \underbrace{Q^{\mathrm{T}} Q}_{I} R z^{*}=R^{\mathrm{T}} Q^{\mathrm{T}} g \quad \Rightarrow \quad R^{\mathrm{T}} R z^{*}=R^{\mathrm{T}} Q^{\mathrm{T}} g .\] Aquí, utilizamos el hecho de que\(Q^{\mathrm{T}} Q=I\) si\(Q\) es una matriz ortonormal. La matriz triangular superior es invertible siempre y cuando sus entradas diagonales sean todas distintas de cero (que es el caso de las columnas linealmente independientes), por lo que podemos simplificar aún más la expresión para rendir\[R z^{*}=Q^{\mathrm{T}} g .\] Así, una vez que la factorización esté disponible, necesitamos formar la derecha\(B\) lado\(Q^{\mathrm{T}} g\), que requiere\(2 m n\) operaciones, y resolver el sistema lineal triangular\(n \times n\) superior, que requiere\(n^{2}\) operaciones. Ambas operaciones son económicas. La mayor parte del costo está en factorizar la matriz\(B\) en matrices\(Q\) y\(R\).

Existen dos clases de métodos para producir una\(Q R\) factorización: el procedimiento Gram-Schmidt y la transformación de Householder. Aquí, discutiremos brevemente el procedimiento de Gram-Schmidt. La idea detrás del procedimiento de Gram-Schmidt es convertir sucesivamente las columnas de\(B\) en vectores ortonormales para formar la matriz ortonormal\(Q\). Por conveniencia, denotamos la\(j\) -ésima columna de\(B\) by\(b_{j}\), es decir,\[B=\left(\begin{array}{llll} b_{1} & b_{2} & \cdots & b_{n} \end{array}\right),\] donde\(b_{j}\) es un\(m\) -vector. De igual manera, expresamos nuestra matriz ortonormal como\[Q=\left(\begin{array}{llll} q_{1} & q_{2} & \cdots & q_{n} \end{array}\right) .\] Recall\(q_{i}^{\mathrm{T}} q_{j}=\delta_{i j}\) (Kronecker-delta),\(1 \leq i, j \leq n\).

El procedimiento de Gram-Schmidt inicia con un conjunto que consiste en un solo vector,\(b_{1}\). Construimos un conjunto ortonormal que consiste en un solo vector\(q_{1}\) que abarca el mismo espacio que\(\left\{b_{1}\right\}\). Trivialmente, podemos tomar\[q_{1}=\frac{1}{\left\|b_{1}\right\|} b_{1} .\] O, podemos expresar\(b_{1}\) como\[b_{1}=q_{1}\left\|b_{1}\right\|,\] que es el producto de un vector unitario y una amplitud.

Ahora consideramos un conjunto que consta de las dos primeras columnas de\(B,\left\{b_{1}, b_{2}\right\}\). Nuestro objetivo es construir un conjunto ortonormal\(\left\{q_{1}, q_{2}\right\}\) que abarque el mismo espacio que\(\left\{b_{1}, b_{2}\right\}\). En particular, mantendremos inalterado el que\(q_{1}\) hemos construido en el primer paso, y elegiremos de\(q_{2}\) tal manera que\((i)\) sea ortogonal a\(q_{1}\), y (ii)\(\left\{q_{1}, q_{2}\right\}\) abarca el mismo espacio que\(\left\{b_{1}, b_{2}\right\}\). Para ello, podemos comenzar con\(b_{2}\) y primero retirar el componente en la dirección de\(q_{1}\), es decir\[\tilde{q}_{2}=b_{2}-\left(q_{1}^{\mathrm{T}} b_{2}\right) q_{1} .\] Aquí, recordamos el hecho de que el producto interno\(q_{1}^{\mathrm{T}} b_{2}\) es el componente de\(b_{2}\) en la dirección de \(q_{1}\). Podemos confirmar fácilmente que\(\tilde{q}_{2}\) es ortogonal a\(q_{1}\), es decir,\[q_{1}^{\mathrm{T}} \tilde{q}_{2}=q_{1}^{\mathrm{T}}\left(b_{2}-\left(q_{1}^{\mathrm{T}} b_{2}\right) q_{1}\right)=q_{1}^{\mathrm{T}} b_{2}-\left(q_{1}^{\mathrm{T}} b_{2}\right) q_{1}^{\mathrm{T}} q_{1}=q_{1}^{\mathrm{T}} b_{2}-\left(q_{1}^{\mathrm{T}} b_{2}\right) \cdot 1=0 .\] finalmente, normalizamos\(\tilde{q}_{2}\) para producir el vector de longitud unitaria\[q_{2}=\tilde{q}_{2} /\left\|\tilde{q}_{2}\right\| .\] Con algún reordenamiento, vemos que se\(b_{2}\) puede expresar como\[b_{2}=\left(q_{1}^{\mathrm{T}} b_{2}\right) q_{1}+\tilde{q}_{2}=\left(q_{1}^{\mathrm{T}} b_{2}\right) q_{1}+\left\|\tilde{q}_{2}\right\| q_{2} .\] Usando un producto matriz-vector, podemos expresar esto como\[b_{2}=\left(\begin{array}{ll} q_{1} & q_{2} \end{array}\right)\left(\begin{array}{c} q_{1}^{\mathrm{T}} b_{2} \\ \left\|\tilde{q}_{2}\right\| \end{array}\right) .\] Combinando con la expresión for\(b_{1}\), tenemos\[\left(\begin{array}{ll} b_{1} & b_{2} \end{array}\right)=\left(\begin{array}{ll} q_{1} & q_{2} \end{array}\right)\left(\begin{array}{cc} \left\|b_{1}\right\| & q_{1}^{\mathrm{T}} b_{2} \\ & \left\|\tilde{q}_{2}\right\| \end{array}\right) .\] En dos pasos, hemos factorizado las dos primeras columnas de\(B\) en una matriz\(m \times 2\) ortogonal\(\left(q_{1}, q_{2}\right)\) y una \(2 \times 2\)matriz triangular superior. El procedimiento de Gram-Schmidt consiste en repetir los\(n\) tiempos del procedimiento; mostremos un paso más para mayor claridad.

En el tercer paso, consideramos un conjunto que consta de las tres primeras columnas de\(B,\left\{b_{1}, b_{2}, b_{3}\right\}\). Nuestro objetivo es construir un conjunto ortonormal\(\left\{q_{1}, q_{2}, q_{3}\right\}\). Siguiendo la misma receta que el segundo paso, mantenemos\(q_{1}\) y\(q_{2}\) sin cambios, y elegimos de\(q_{3}\) tal manera que\((i)\) sea ortogonal a\(q_{1}\) y\(q_{2}\), y (ii)\(\left\{q_{1}, q_{2}, q_{3}\right\}\) abarca el mismo espacio que\(\left\{b_{1}, b_{2}, b_{3}\right\}\). Esta vez, partimos de\(b_{3}\), y eliminamos los componentes de\(b_{3}\) en la dirección de\(q_{1}\) y\(q_{2}\), es decir\[\tilde{q}_{3}=b_{3}-\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{1}-\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{2} .\] Una vez más, recordamos que\(q_{1}^{\mathrm{T}} b_{3}\) y\(q_{2}^{\mathrm{T}} b_{3}\) son los componentes de \(b_{3}\)en la dirección de\(q_{1}\) y\(q_{2}\), respectivamente. Podemos volver a confirmar que\(\tilde{q}_{3}\) es ortogonal a\(q_{1}\)\[q_{1}^{\mathrm{T}} \tilde{q}_{3}=q_{1}^{\mathrm{T}}\left(b_{3}-\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{1}-\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{2}\right)=q_{1}^{\mathrm{T}} b_{3}-\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{1}^{\mathrm{T}} \not{q_{1}}-\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{1}^{\mathrm{T}} \not{q_{2}}=0\] y a\(q_{2}\)\[q_{2}^{\mathrm{T}} \tilde{q}_{3}=q_{2}^{\mathrm{T}}\left(b_{3}-\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{1}-\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{2}\right)=q_{2}^{\mathrm{T}} b_{3}-\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{2}^{\mathrm{T}} \not{q_{1}}-\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{2}^{\mathrm{T}} \not{q_{2}}=0\] Podemos expresarnos\(b_{3}\) como\[b_{3}=\left(q_{1}^{\mathrm{T}} b_{3}\right) q_{1}+\left(q_{2}^{\mathrm{T}} b_{3}\right) q_{2}+\left\|\tilde{q}_{3}\right\| q_{3} .\] O, juntando las tres primeras columnas\[\left(\begin{array}{lll} b_{1} & b_{2} & b_{3} \end{array}\right)=\left(\begin{array}{lll} q_{1} & q_{2} & q_{3} \end{array}\right)\left(\begin{array}{ccc} \left\|b_{1}\right\| & q_{1}^{\mathrm{T}} b_{2} & q_{1}^{\mathrm{T}} b_{3} \\ & \left\|\tilde{q}_{2}\right\| & q_{2}^{\mathrm{T}} b_{3} \\ & & \left\|\tilde{q}_{3}\right\| \end{array}\right)\] Podemos ver que repitiendo el procedimiento \(n\)los tiempos resultarían en la ortogonalización completa de las columnas de\(B\).

Contemos el número de operaciones del procedimiento Gram-Schmidt. En el paso\(j\) -ésimo, hay\(j-1\) componentes a eliminar, cada uno requiriendo de\(4 m\) operaciones. Así, el conteo total de operaciones es\[C^{\text {Gram-Schmidt }} \approx \sum_{j=1}^{n}(j-1) 4 m \approx 2 m n^{2} .\] Así, para la solución del problema de mínimos cuadrados, el método basado en Gram-Schmidt es aproximadamente el doble de costoso que el método basado en la ecuación normal para\(m \gg n\). Sin embargo, la estabilidad numérica superior a menudo justifica el costo adicional.

Observamos que existe una versión modificada de Gram-Schmidt, llamada procedimiento Gram-Schmidt modificado, que es más estable que el algoritmo presentado anteriormente. El procedimiento Gram-Schmidt modificado requiere el mismo costo computacional. También existe otro algoritmo de\(Q R\) factorización fundamentalmente diferente, llamado la transformación Householder, que es aún más estable que el procedimiento modificado de Gram-Schmidt. El algoritmo Householder requiere aproximadamente el mismo costo que el procedimiento de Gram-Schmidt.

Material Avanzado

Comenzar material avanzado

Interpretación de mínimos cuadrados: Proyección

Hasta el momento, hemos discutido un procedimiento para resolver un sistema sobredeterminado,\[B z=g,\] en el sentido de mínimos cuadrados. Usando la interpretación de columna del producto matriz-vector, estamos buscando la combinación lineal de las columnas de\(B\) que minimice la norma 2 del residual - el desajuste entre una representación\(B z\) y los datos\(g\). La solución de mínimos cuadrados al problema es Es\[B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}} g \quad \Rightarrow \quad z^{*}=\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}} g .\] decir, la aproximación más cercana de los datos\(g\) usando las columnas de\(B\) es\[g^{\mathrm{LS}}=B z^{*}=B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}} g=P^{\mathrm{LS}} g .\] Nuestra mejor representación de\(g, g^{\mathrm{LS}}\), es la proyección de\(g\) por el proyector \(P^{\mathrm{LS}}\). Podemos verificar que el operador\(P^{\mathrm{LS}}=B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}}\) es efectivamente un proyector:\[\begin{aligned} \left(P^{\mathrm{LS}}\right)^{2}=\left(B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}}\right)^{2} &=B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}} B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}}=B \underbrace{\left(\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}} B\right)}_{I}\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}} \\ &=B\left(B^{\mathrm{T}} B\right)^{-1} B^{\mathrm{T}}=P^{\mathrm{LS}} . \end{aligned}\] De hecho,\(P^{\mathrm{LS}}\) es un proyector ortogonal porque\(P^{\mathrm{LS}}\) es simétrico. Esto concuerda con nuestra intuición; la representación más cercana del\(g\) uso de las columnas de\(B\) resultados de\(g\) proyectarse\(\operatorname{col}(B)\) a lo largo de un espacio ortogonal a\(\operatorname{col}(B)\). Esto queda claramente demostrado para\(\mathbb{R}^{2}\) en la Figura\(17.5\) considerada anteriormente.

Usando el proyector ortogonal sobre\(\operatorname{col}(B), P^{\mathrm{LS}}\), podemos pensar en otra interpretación de la solución de mínimos cuadrados. Primero proyectamos los datos\(g\) ortogonalmente al espacio de columna para formar\[g^{\mathrm{LS}}=P^{\mathrm{LS}} g .\] Luego, encontramos los coeficientes para la combinación lineal de las columnas de\(B\) ese resultado\(P^{\mathrm{LS}} g\), es decir,\[B z^{*}=P^{\mathrm{LS}} g .\] Este problema tiene una solución porque \(P^{\mathrm{LS}} g \in \operatorname{col}(B)\).

Esta interpretación es útil especialmente cuando la\(Q R\) factorización de\(B\) está disponible. Si\(B=Q R\), entonces\(\operatorname{col}(B)=\operatorname{col}(Q)\). Entonces, el proyector ortogonal sobre\(\operatorname{col}(B)\) es lo mismo que el proyector ortogonal sobre\(\operatorname{col}(Q)\) y es dado por\[P^{\mathrm{LS}}=Q Q^{\mathrm{T}} .\] Podemos verificar que de hecho\(P^{\mathrm{LS}}\) es un proyector ortogonal comprobando que es (i) idempotente\(\left(P^{\mathrm{LS}} P^{\mathrm{LS}}=P^{\mathrm{LS}}\right)\), y \((i i)\)simétrico\(\left(\left(P^{\mathrm{LS}}\right)^{\mathrm{T}}=P^{\mathrm{LS}}\right)\), es decir,\[\begin{aligned} &P^{\mathrm{LS}} P^{\mathrm{LS}}=\left(Q Q^{\mathrm{T}}\right)\left(Q Q^{\mathrm{T}}\right)=Q \underbrace{Q^{\mathrm{T}} Q}_{I} Q^{\mathrm{T}}=Q Q^{\mathrm{T}}=P^{\mathrm{LS}}, \\ &\left(P^{\mathrm{LS}}\right)^{\mathrm{T}}=\left(Q Q^{\mathrm{T}}\right)^{\mathrm{T}}=\left(Q^{\mathrm{T}}\right)^{\mathrm{T}} Q^{\mathrm{T}}=Q Q^{\mathrm{T}}=P^{\mathrm{LS}} . \end{aligned}\] Usando la\(Q R\) factorización de\(B\), podemos reescribir la solución de mínimos cuadrados como\[B z^{*}=P^{\mathrm{LS}} g \quad \Rightarrow \quad Q R z^{*}=Q Q^{\mathrm{T}} g .\] Aplicando\(Q^{\mathrm{T}}\) en ambos lados y usando el hecho de que \(Q^{\mathrm{T}} Q=I\), obtenemos\[R z^{*}=Q^{\mathrm{T}} g .\] Geométricamente, estamos proyectando ortogonalmente los datos\(g\) sobre\(\operatorname{col}(Q)\) pero representando la solución proyectada en la base\(\left\{q_{i}\right\}_{i=1}^{n}\) del espacio\(n\) -dimensional (en lugar de en la base estándar de\(\mathbb{R}^{m}\)). Entonces, encontramos los coeficientes\(z^{*}\) que arrojan los datos proyectados.

Material Avanzado

Comenzar material avanzado

límites de error para mínimos cuadrados

Quizás la forma más obvia de medir la bondad de nuestra solución es en términos del residuo\(\left\|g-B z^{*}\right\|\) que indica hasta qué punto\(B z^{*}=g\) se satisfacen las ecuaciones, qué tan bien\(B z^{*}\) predice\(g\). Ya que elegimos\(z^{*}\) minimizar\(\left\|g-B z^{*}\right\|\) podemos esperar que\(\left\|g-B z^{*}\right\|\) sea pequeño. Pero es importante reconocer que en la mayoría de los casos\(g\) solo refleja datos de un experimento en particular mientras que nos gustaría luego usar nuestra predicción para\(z^{*}\) en otros, diferentes, experimentos o incluso contextos. Por ejemplo, el coeficiente de fricción que medimos en el laboratorio se utilizará posteriormente “en el campo” como parte de una predicción del sistema más grande para, digamos, el rendimiento del robot. En este sentido, no sólo el residual podría no ser una buena medida del “error en”\(z\), un residuo menor podría ni siquiera implicar una “mejor predicción” para\(z\). En esta sección, analizamos cómo el ruido y los modelos incompletos (sesgo) se pueden relacionar directamente con nuestra predicción para\(z\).

Tenga en cuenta que, para simplificar la notación, utilizamos el subíndice 0 para representar el superíndice “verdadero” en esta sección.

Límite de error con respecto a la perturbación en los datos,\(g\) (modelo constante)

Consideremos un ajuste de parámetros para un modelo constante simple. Primero, supongamos que hay una solución\(z_{0}\) al sistema sobredeterminado\[\underbrace{\left(\begin{array}{c} 1 \\ 1 \\ \vdots \\ 1 \end{array}\right)}_{B} z_{0}=\underbrace{\left(\begin{array}{c} g_{0,1} \\ g_{0,2} \\ \vdots \\ g_{0, m} \end{array}\right)}_{g_{0}} .\] Porque\(z_{0}\) es la solución al sistema,\(g_{0}\) debe ser un múltiplo constante de\(B\). Es decir, las entradas de\(g_{0}\) deben ser todas iguales. Ahora, supongamos que los datos están perturbados de tal manera que\(g \neq g_{0}\). Con los datos perturbados, es poco probable que el sistema sobredeterminado tenga una solución, por lo que consideramos la solución de mínimos cuadrados\(z^{*}\) al problema\[B z=g .\] Nos gustaría saber cuánta perturbación en los datos\(g-g_{0}\) cambia la solución\(z^{*}-z_{0}\).

Para cuantificar el efecto de la perturbación, primero observamos que tanto la solución original como la solución al sistema perturbado satisfacen la ecuación normal, es decir,\[B^{\mathrm{T}} B z_{0}=B^{\mathrm{T}} g_{0} \quad \text { and } \quad B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}} g .\] Tomando la diferencia de las dos expresiones, obtenemos\[B^{\mathrm{T}} B\left(z^{*}-z_{0}\right)=B^{\mathrm{T}}\left(g-g_{0}\right) .\] For\(B\) con el modelo constante, tenemos \(B^{\mathrm{T}} B=m\), simplificando la expresión a\[\begin{aligned} z^{*}-z_{0} &=\frac{1}{m} B^{\mathrm{T}}\left(g-g_{0}\right) \\ &=\frac{1}{m} \sum_{i=1}^{m}\left(g-g_{0}\right)_{i} . \end{aligned}\] Así, si el “ruido” está cerca de la media cero,\(z^{*}\) está cerca de\(Z_{0}\). De manera más general, podemos demostrar que\[\left|z^{*}-z_{0}\right| \leq \frac{1}{\sqrt{m}}\left\|g-g_{0}\right\| .\] Vemos que la desviación en la solución está delimitada por los datos de perturbación. Por lo tanto, nuestra solución de mínimos cuadrados\(z^{*}\) es una buena aproximación siempre que la perturbación\(\left\|g-g_{0}\right\|\) sea pequeña.

Para probar este resultado, aplicamos la desigualdad Cauchy-Schwarz, es decir,\[\left|z^{*}-z_{0}\right|=\frac{1}{m}\left|B^{\mathrm{T}}\left(g-g_{0}\right)\right| \leq \frac{1}{m}\|B\|\left\|g-g_{0}\right\|=\frac{1}{m} \sqrt{m}\left\|g-g_{0}\right\|=\frac{1}{\sqrt{m}}\left\|g-g_{0}\right\| .\] Recordemos que la desigualdad Cauchy-Schwarz da un límite bastante pesimista cuando los dos vectores no están muy bien alineados.

Estudiemos ahora de manera más formal cómo la alineación de los dos vectores\(B\) y\(g-g_{0}\) afecta el error en la solución. Para cuantificar el efecto recordemos que la solución de mínimos cuadrados satisface\[B z^{*}=P^{\mathrm{LS}} g,\] dónde\(P^{\mathrm{LS}}\) está el proyector ortogonal sobre el espacio de columna de\(B, \operatorname{col}(B)\). Si\(g-g_{0}\) es exactamente la media cero, es decir,\[\frac{1}{m} \sum_{i=1}^{m}\left(g_{0, i}-g_{i}\right)=0,\] entonces\(g-g_{0}\) es ortogonal a\(\operatorname{col}(B)\). Debido a que cualquier perturbación ortogonal a\(\operatorname{col}(B)\) se encuentra en la dirección a lo largo de la cual se realiza la proyección, no afecta\(P^{\mathrm{LS}} g\) (y por lo tanto\(B z^{*}\)), y en particular\(z^{*}\). Es decir, la solución de mínimos cuadrados,\(z^{*}\), to\[B z=g=g_{0}+\left(g-g_{0}\right)\] es\(z_{0}\) si\(g-g_{0}\) tiene media cero. También podemos demostrar que la perturbación de la media cero no tiene influencia en la solución algebraicamente usando la ecuación normal, es decir\[B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}}\left(g_{0}+\left(g-g_{0}\right)\right)=B^{\mathrm{T}} g_{0}+\underline{B}^{\mathrm{T}}\left(g-g_{0}\right)=B^{\mathrm{T}} g_{0} .\] Los datos perturbados\(g\) no ingresan al cálculo de\(z^{*}\) si\(g-g_{0}\) tiene media cero. Por lo tanto, cualquier error en la solución\(z-z_{0}\) debe deberse a la perturbación no media cero en los datos. En consecuencia, el límite basado en la desigualdad Cauchy-Schwarz es más bien pesimista cuando la perturbación se acerca a la media cero.

Límite de error con respecto a la perturbación en los datos,\(g\) (general)

Generalicemos ahora el análisis de perturbación a un sistema general sobredeterminado,\[B z_{0}=g_{0},\] donde\(B \in \mathbb{R}^{m \times n}\) con\(m>n\). Suponemos que\(g_{0}\) se elige de tal manera que existe la solución al sistema lineal. Ahora digamos que el error de medición ha corrompido\(g_{0}\) a\(g=g_{0}+\epsilon\). En particular, suponemos que el sistema lineal\[B z=g\] no tiene solución. Así, en cambio, encontramos la solución\(z^{*}\) de mínimos cuadrados del sistema.

Para establecer los límites de error, primero introduciremos el concepto de valores singulares máximos y mínimos, que nos ayudan a caracterizar el comportamiento de\(B\). Los valores máximos y mínimos singulares de\(B\) son definidos por\[\nu_{\max }(B)=\max _{v \in \mathbb{R}^{n}} \frac{\|B v\|}{\|v\|} \quad \text { and } \quad \nu_{\min }(B)=\min _{v \in \mathbb{R}^{n}} \frac{\|B v\|}{\|v\|} .\] Tenga en cuenta que, debido a que la norma escala linealmente bajo multiplicación escalar, las definiciones equivalentes de los valores singulares son\[\nu_{\max }(B)=\max _{\substack{v \in \mathbb{R}^{n} \\\|v\|=1}}\|B v\| \quad \text { and } \quad \nu_{\min }(B)=\min _{\substack{v \in \mathbb{R}^{n} \\\|v\|=1}}\|B v\| .\] En otras palabras, el valor singular máximo es el estiramiento máximo que \(B\)puede inducir a un vector unitario. De igual manera, el valor singular mínimo es la contracción máxima\(B\) que puede inducir. En particular, recordemos que si las columnas de no\(B\) son linealmente independientes, entonces podemos encontrar una no trivial\(v\) para la cual\(B v=0\). Así, si las columnas de\(B\) son linealmente dependientes,\(\nu_{\min }(B)=0\).

También observamos que los valores singulares están relacionados con los valores propios de\(B^{\mathrm{T}} B\). Recordemos que 2-norma está relacionada con el producto interno por\[\|B v\|^{2}=(B v)^{\mathrm{T}}(B v)=v^{\mathrm{T}} B^{\mathrm{T}} B v,\] lo tanto, a partir del cociente Rayleigh, la raíz cuadrada de los valores propios máximos y mínimos de\(B^{\mathrm{T}} B\) son los valores máximos y mínimos singulares de\(B\).

Cuantifiquemos la sensibilidad del error de solución al lado derecho de dos maneras diferentes. Primero está el condicionamiento absoluto, que se\(\left\|z^{*}-z_{0}\right\|\) relaciona con\(\left\|g-g_{0}\right\|\). El límite viene dado por\[\left\|z^{*}-z_{0}\right\| \leq \frac{1}{\nu_{\min }(B)}\left\|g-g_{0}\right\| .\] Segundo es el condicionamiento relativo, que relaciona la perturbación relativa en la solución\(\| z^{*}-\)\(z_{0}\|/\| z_{0} \|\) y la perturbación relativa en el lado derecho\(\left\|g-g_{0}\right\| /\left\|g_{0}\right\|\). Este encuadernado es dado por\[\frac{\left\|z^{*}-z_{0}\right\|}{\left\|z_{0}\right\|}=\frac{\nu_{\max }(B)}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|} .\] Derivamos estos resultados en breve.

Si la perturbación\(\left\|g-g_{0}\right\|\) es pequeña, esperamos que el error\(\left\|z^{*}-z_{0}\right\|\) sea pequeño siempre y cuando\(B\) esté bien condicionado en el sentido de que no\(\nu_{\max }(B) / \nu_{\min }(B)\) sea demasiado grande. Obsérvese que si\(B\) tiene columnas linealmente dependientes, entonces\(\nu_{\min }=0\) y\(\nu_{\max } / \nu_{\min }\) es infinito; así,\(\nu_{\max } / \nu_{\min }\) es una medida de la independencia de las columnas de\(B\) y de ahí la medida en que podemos determinar de forma independiente los diferentes elementos de\(z\). De manera más general,\(\nu_{\max } / \nu_{\min }\) es una medida de la sensibilidad o estabilidad de nuestras soluciones de mínimos cuadrados a las perturbaciones (por ejemplo, in\(g\)). Como ya hemos visto en este capítulo, y volveremos a ver en el Capítulo 19 dentro del contexto de regresión, podemos hasta cierto punto “controlar”\(B\) a través de la elección de variables, dependencias funcionales y puntos de medición (lo que da lugar al importante campo del “diseño del experimento (s)”); así podemos esforzarnos por controlar a\(\nu_{\max } / \nu_{\min }\) través de buenas elecciones “independientes” y así asegurar una buena predicción de\(z\).

Ejemplo 17.3.4 Medición de Ruido en Ajuste Polinómico

Demostremos el efecto de la perturbación en\(g\) - o el error de medición - en el contexto

a) gran perturbación

b) pequeñas perturbaciones

Figura 17.6: Efecto de la perturbación de datos sobre la solución.

de ajuste polinomial que consideramos anteriormente. Como antes, asumimos que la salida depende de la entrada cuadráticamente según\[y(x)=-\frac{1}{2}+\frac{2}{3} x-\frac{1}{8} c x^{2}\] con\(c=1\). Construimos datos limpios\(g_{0} \in \mathbb{R}^{m}, m=7\), evaluando\(y\) en\[x_{i}=(i-1) / 2, \quad i=1, \ldots, m\] y configurando\[g_{0, i}=y\left(x_{i}\right), \quad i=1, \ldots, m\] Porque sigue\(g_{0}\) con precisión la función cuadrática,\(z_{0}=(-1 / 2,2 / 3,-1 / 8)\) satisface el sistema sobredeterminado\(B z_{0}=g_{0}\). Recordemos que\(B\) es la matriz\(m \times n\) de Vandermonde con los puntos de evaluación\(\left\{x_{i}\right\}\).

Luego construimos datos\(g\) perturbados agregando ruido aleatorio a\(g_{0}\), es decir,\[g_{i}=g_{0, i}+\epsilon_{i}, \quad i=1, \ldots, m\] Entonces, resolvemos para la solución\(z^{*}\) de mínimos cuadrados de\(B z^{*}=g\).

El resultado de resolver el problema de ajuste polinomial para dos niveles de perturbación diferentes se muestra en la Figura\(17.6\). Para el caso de gran perturbación, la perturbación en los datos y el error en la solución -ambos medidos en 2-norma- son\[\left\|g-g_{0}\right\|=0.223 \text { and }\left\|z-z_{0}\right\|=0.072\] En contraste, el caso de pequeña perturbación produce\[\left\|g-g_{0}\right\|=0.022 \text { and }\left\|z-z_{0}\right\|=0.007\] Los resultados confirman que una menor perturbación en los datos da como resultado un error menor en la solución. También podemos verificar los límites de error. El valor singular mínimo de la matriz de Vandermonde es\[\nu_{\min }(B)=0.627 .\] Aplicación del error (absoluto) ligado al gran caso de perturbación rendimientos\[0.072=\left\|z-z_{0}\right\| \leq \frac{1}{\nu_{\min }(B)}\left\|g-g_{0}\right\|=0.356 .\] El límite de error se satisface claramente. El error límite para el caso de pequeña perturbación se satisface de manera similar.

Ahora probamos los límites de error.

Prueba. Para establecer el límite de error absoluto, primero observamos que la solución al problema limpio,\(z_{0}\), y la solución al problema perturbado,\(z^{*}\), satisfacen la ecuación normal, es decir,\[B^{\mathrm{T}} B z_{0}=B^{\mathrm{T}} g_{0} \quad \text { and } \quad B^{\mathrm{T}} B z^{*}=B^{\mathrm{T}} g .\] Tomando la diferencia de las dos ecuaciones\[B^{\mathrm{T}} B\left(z^{*}-z_{0}\right)=B^{\mathrm{T}}\left(g-g_{0}\right) .\] Ahora, multiplicamos ambos lados por\(\left(z^{*}-z_{0}\right)^{\mathrm{T}}\) para obtener\[\begin{aligned} \text { (LHS }) &=\left(z^{*}-z_{0}\right)^{\mathrm{T}} B^{\mathrm{T}} B\left(z^{*}-z_{0}\right)=\left(B\left(z^{*}-z_{0}\right)\right)^{\mathrm{T}}\left(B\left(z^{*}-z_{0}\right)\right)=\left\|B\left(z^{*}-z_{0}\right)\right\|^{2} \\ \text { (RHS }) &=\left(z^{*}-z_{0}\right)^{\mathrm{T}} B^{\mathrm{T}}\left(g-g_{0}\right)=\left(B\left(z^{*}-z_{0}\right)\right)^{\mathrm{T}}\left(g-g_{0}\right) \leq\left\|B\left(z^{*}-z_{0}\right)\right\|\left\|g-g_{0}\right\|, \end{aligned}\] donde hemos invocado la desigualdad Cauchy-Schwarz en el lado derecho. Así, tenemos\[\left\|B\left(z^{*}-z_{0}\right)\right\|^{2} \leq\left\|B\left(z^{*}-z_{0}\right)\right\|\left\|g-g_{0}\right\| \quad \Rightarrow \quad\left\|B\left(z^{*}-z_{0}\right)\right\| \leq\left\|g-g_{0}\right\| .\] Podemos encuadernar el lado izquierdo desde abajo usando la definición del valor mínimo singular\[\nu_{\min }(B)\left\|z^{*}-z_{0}\right\| \leq\left\|B\left(z^{*}-z_{0}\right)\right\| .\] Así, tenemos\[\nu_{\min }\left\|z^{*}-z_{0}\right\| \leq\left\|B\left(z^{*}-z_{0}\right)\right\| \leq\left\|g-g_{0}\right\| \quad \Rightarrow \quad\left\|z^{*}-z_{0}\right\| \leq \frac{1}{\nu_{\min }(B)}\left\|g-g_{0}\right\|,\] cual es el límite de error absoluto deseado.

Para obtener el límite de error relativo, primero dividimos el error absoluto enlazado por\(\left\|z_{0}\right\|\) para obtener\[\frac{\left\|z^{*}-z_{0}\right\|}{\left\|z_{0}\right\|} \leq \frac{1}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|z_{0}\right\|}=\frac{1}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|} \frac{\left\|g_{0}\right\|}{\left\|z_{0}\right\|} .\] Para encuadernar el cociente\(\left\|g_{0}\right\| /\left\|z_{0}\right\|\), tomamos la norma de ambos lados de\(B z_{0}=g_{0}\) e invocamos la definición del valor singular máximo, i.e.\[\left\|g_{0}\right\|=\left\|B z_{0}\right\| \leq \nu_{\max }\left\|z_{0}\right\| \Rightarrow \frac{\left\|g_{0}\right\|}{\left\|z_{0}\right\|} \leq \nu_{\max } .\] Sustituyendo la expresión por el límite anterior\[\frac{\left\|z^{*}-z_{0}\right\|}{\left\|z_{0}\right\|} \leq \frac{1}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|} \frac{\left\|g_{0}\right\|}{\left\|z_{0}\right\|} \leq \frac{\nu_{\max }(B)}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|},\] que es el límite de error relativo deseado. Prueba (usando descomposición de valores singulares). Comenzamos con la descomposición del valor singular de la matriz\(B\),\[B=U \Sigma V^{\mathrm{T}}\] donde\(U\) es una matriz\(m \times m\) unitaria,\(V\) es una matriz\(n \times n\) unitaria, y\(\Sigma\) es una\(m \times n\) matriz diagonal. En particular,\(\Sigma\) consiste en valores singulares de\(B\) y es de la forma\[\Sigma=\left(\begin{array}{cccc} \nu_{1} & & & \\ & \nu_{2} & & \\ & & \ddots & \\ & & & \nu_{n} \\ & & & \end{array}\right)=\left(\begin{array}{c} \widehat{\Sigma} \\ 0 \end{array}\right)\] La descomposición del valor singular existe para cualquier matriz. La solución al problema original viene dada por\[B z=g \quad \Rightarrow \quad U \Sigma V^{\mathrm{T}} z=g \quad \Rightarrow \quad \Sigma V^{\mathrm{T}} z=U^{\mathrm{T}} g .\] La solución al problema de mínimos cuadrados es\[\begin{aligned} z^{*} &=\arg \min _{z}\|B z-g\|=\arg \min _{z}\left\|U \Sigma V^{\mathrm{T}} z-g\right\|=\arg \min _{z}\left\|\Sigma V^{\mathrm{T}} z-U^{\mathrm{T}} g\right\| \\ &=V\left(\arg \min _{\tilde{z}}\|\Sigma \tilde{z}-\tilde{g}\|\right), \end{aligned}\] donde la tercera igualdad se desprende del hecho de que la acción por una matriz unitaria no altera la norma 2, y hemos hecho las sustituciones\(\tilde{z}=V^{\mathrm{T}} z\) y \(\tilde{g}=U^{\mathrm{T}} g\). Observamos que debido a que\(\Sigma\) es diagonal, la norma 2 a minimizar es particularmente simple,\[\Sigma \tilde{z}-\tilde{g}=\Sigma=\left(\begin{array}{ccc} \nu_{1} & & \\ & \ddots & \\ & & \nu_{n} \\ & & \\ & & \end{array}\right)\left(\begin{array}{c} \tilde{z}_{1} \\ \vdots \\ \tilde{z}_{n} \end{array}\right)-\left(\begin{array}{c} \tilde{g}_{1} \\ \vdots \\ \tilde{g}_{n} \\ \tilde{g}_{n+1} \\ \vdots \\ \tilde{g}_{m} \end{array}\right) .\] Obsérvese que elegir\(\tilde{z}_{1}, \ldots, \tilde{z}_{n}\) sólo afecta al primer\(n\) componente del vector residual. Por lo tanto, debemos escoger\(\tilde{z}_{1}, \ldots, \tilde{z}_{n}\) tal que\[\left(\begin{array}{ccc} \nu_{1} & & \\ & \ddots & \\ & & \nu_{n} \end{array}\right)\left(\begin{array}{c} \tilde{z}_{1} \\ \vdots \\ \tilde{z}_{n} \end{array}\right)=\left(\begin{array}{c} \tilde{g}_{1} \\ \vdots \\ \tilde{g}_{n} \end{array}\right) \Rightarrow \quad \tilde{z}_{i}=\frac{\tilde{g}_{i}}{\nu_{i}}, \quad i=1, \ldots, n .\] Al introducir una matriz de\(n \times m\) restricción que extrae las primeras\(n\) entradas de\(\tilde{g}\), podemos escribir concisamente lo anterior como\[\widehat{\Sigma} \tilde{z}=R \tilde{g} \quad \Rightarrow \quad \tilde{z}=\widehat{\Sigma}^{-1} R \tilde{g},\] y la solución al problema de mínimos cuadrados as\[z^{*}=V \tilde{z}^{*}=V \widehat{\Sigma}^{-1} R \tilde{g}=V \widehat{\Sigma}^{-1} R U^{\mathrm{T}} g .\] El límite de número de condición absoluta se obtiene por\[\begin{aligned} \left\|z^{*}-z_{0}\right\| &=\left\|V \widehat{\Sigma}^{-1} R U^{\mathrm{T}}\left(g-g_{0}\right)\right\|=\frac{\left\|V \widehat{\Sigma}^{-1} R U^{\mathrm{T}}\left(g-g_{0}\right)\right\|}{\left\|g-g_{0}\right\|}\left\|g-g_{0}\right\| \\ & \leq\left(\sup _{\delta g} \frac{\left\|V \widehat{\Sigma}^{-1} R U^{\mathrm{T}} \delta g\right\|}{\|\delta g\|}\right)\left\|g-g_{0}\right\| . \end{aligned}\] El término entre paréntesis se limita al señalar que las transformaciones ortogonales preservan la norma 2 y que el operador de restricción no aumenta la norma 2,\[\sup _{\delta g}\left(\frac{\left\|V \widehat{\Sigma}^{-1} R U^{\mathrm{T}} \delta g\right\|}{\|\delta g\|}\right)=\sup _{\delta \tilde{g}}\left(\frac{\left\|V \widehat{\Sigma}^{-1} R \delta \tilde{g}\right\|}{\|U \delta \tilde{g}\|}\right)=\sup _{\delta \tilde{g}}\left(\frac{\left\|\widehat{\Sigma}^{-1} R \delta \tilde{g}\right\|}{\|\delta \tilde{g}\|}\right) \leq \frac{1}{\nu_{\min }(B)} .\] es decir, tenemos el límite de error absoluto deseado \[\left\|z^{*}-z_{0}\right\| \leq \frac{1}{\nu_{\min }(B)}\left\|g-g_{0}\right\| .\]Ahora consideremos el límite de error relativo. Primero observamos que\[\frac{\left\|z^{*}-z_{0}\right\|}{\left\|z_{0}\right\|}=\frac{1}{\nu_{\min }(B)}\left\|g-g_{0}\right\| \frac{1}{\left\|z_{0}\right\|}=\frac{1}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|} \frac{\left\|g_{0}\right\|}{\left\|z_{0}\right\|} .\] El término se\(\left\|g_{0}\right\| /\left\|z_{0}\right\|\) puede delimitar expresando\(z_{0}\) en términos de\(g\) usar la expresión explícita para la solución de mínimos cuadrados, es decir\[\frac{\left\|g_{0}\right\|}{\left\|z_{0}\right\|}=\frac{\left\|B z_{0}\right\|}{\left\|z_{0}\right\|}=\frac{\left\|U \Sigma V^{\mathrm{T}} z_{0}\right\|}{\left\|z_{0}\right\|} \leq \sup _{z} \frac{\left\|U \Sigma V^{\mathrm{T}} z\right\|}{\|z\|}=\sup _{\tilde{z}} \frac{\|U \Sigma \tilde{z}\|}{\|V \tilde{z}\|}=\sup _{\tilde{z}} \frac{\|\Sigma \tilde{z}\|}{\|\tilde{z}\|}=\nu_{\max }(B) .\], tenemos el límite de error relativo\[\frac{\left\|z^{*}-z_{0}\right\|}{\left\|z_{0}\right\|} \leq \frac{\nu_{\max }(B)}{\nu_{\min }(B)} \frac{\left\|g-g_{0}\right\|}{\left\|g_{0}\right\|} .\] Esto concluye el prueba.

Límite de error con respecto a la reducción en el espacio,\(B\)

Consideremos ahora un escenario que ilustre el efecto del sesgo. Nuevamente, comenzamos con un sistema lineal sobredeterminado,\[B_{0} z_{0}=g,\] donde\(B_{0} \in \mathbb{R}^{m \times n}\) con\(m>n\). Suponemos que\(z_{0}\) satisface todas las\(m\) ecuaciones. Recordamos que, en el contexto del ajuste polinómico,\(B_{0}\) es de la forma,\[B_{0}=\left(\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} & \cdots & x_{1}^{n} \\ 1 & x_{2} & x_{2}^{2} & \cdots & x_{2}^{n} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{m} & x_{m}^{2} & \cdots & x_{m}^{n} \end{array}\right),\] donde\(m\) está el número de puntos de datos y\(n\) es el grado de polinomio. Ahora, supongamos que decidimos usar un polinomio de grado\(p\) -ésimo en lugar del polinomio de grado\(n\) -ésimo, donde\(p<n\). En otras palabras, podemos particionar\(B_{0}\) en\[B_{0}=\left(B_{\mathrm{I}} \mid B_{\mathrm{II}}\right)=\left(\begin{array}{cccc|ccc} 1 & x_{1}^{1} & \cdots & x_{1}^{p} & x_{1}^{p+1} & \cdots & x_{1}^{n} \\ 1 & x_{2}^{1} & \cdots & x_{2}^{p} & x_{2}^{p+1} & \cdots & x_{m}^{n} \\ \vdots & \vdots & & \vdots & \vdots & & \vdots \\ 1 & x_{m}^{1} & \cdots & x_{m}^{p} & x_{m}^{p+1} & \cdots & x_{m}^{n} \end{array}\right)\] dónde\(B_{\mathrm{I}} \in \mathbb{R}^{m \times(p+1)}\) y\(B_{\mathrm{II}} \in \mathbb{R}^{m \times(n-p)}\). Entonces podemos resolver el problema de mínimos cuadrados resultante de la primera partición, es decir,\[B_{\mathrm{I}} z^{*}=g .\] Por conveniencia, particionemos también la solución al sistema original\(z_{0}\) en dos partes correspondientes a\(B_{\mathrm{I}}\) y\(B_{\mathrm{II}}\), es decir,\[z_{0}=\left(\begin{array}{c} z_{\mathrm{I}} \\ z_{\mathrm{II}} \end{array}\right),\] donde \(z_{\mathrm{I}} \in \mathbb{R}^{p+1}\)y\(z_{\mathrm{II}} \in \mathbb{R}^{n-p}\). La pregunta es, ¿qué tan cerca están los coeficientes\(z^{*}=\left(z_{1}, \ldots, z_{p-1}\right)\) del sistema reducido en comparación con los coeficientes de la primera partición del sistema original,\(z_{1}\)?

De hecho, podemos encuadernar el error en la solución\(\left\|z^{*}-z_{\mathrm{I}}\right\|\) en términos del “espacio faltante”\(B_{\mathrm{II}}\). En particular, el límite de error absoluto viene dado por\[\left\|z^{*}-z_{\mathrm{I}}\right\| \leq \frac{1}{\nu_{\min }\left(B_{\mathrm{I}}\right)}\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|\] y el límite de error relativo viene dado por\[\frac{\left\|z^{*}-z_{\mathrm{I}}\right\|}{\left\|z_{\mathrm{I}}\right\|} \leq \frac{\nu_{\max }\left(B_{\mathrm{I}}\right)}{\nu_{\min }\left(B_{\mathrm{I}}\right)} \frac{\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|}{\left\|g-B_{\mathrm{II}} z_{\mathrm{II}}\right\|},\] donde\(\nu_{\min }\left(B_{\mathrm{I}}\right)\) y\(\nu_{\max }\left(B_{\mathrm{I}}\right)\) son los valores mínimos y máximos singulares de\(B_{\mathrm{I}}\).

Ejemplo 17.3.5 Efecto de sesgo en el ajuste polinomial

Demostremos el efecto del espacio reducido de solución -o el efecto sesgo- en el contexto del ajuste polinómico. Como antes, la salida depende de la entrada cuadráticamente según\[y(x)=-\frac{1}{2}+\frac{2}{3} x-\frac{1}{8} c x^{2} .\] Recall que\(c\) controla la fuerza de la dependencia cuadrática. Los datos\(g\) se generan evaluando\(y\) en\(x_{i}=(i-1) / 2\) y configurando\(g_{i}=y\left(x_{i}\right)\) para\(i=1, \ldots, m\), con\(m=7\). Particionamos nuestra matriz de Vandermonde para el modelo cuadrático\(B_{0}\) en aquella para el modelo afín\(B_{\mathrm{I}}\) y la única parte cuadrática\(B_{\mathrm{II}}\), i.e.\[B_{0}=\left(\begin{array}{ccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{m} & x_{m}^{2} \end{array}\right)=\left(\begin{array}{cc|c} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ \vdots & \vdots & \vdots \\ 1 & x_{m} & x_{m}^{2} \end{array}\right)=\left(B_{\mathrm{I}} \mid B_{\mathrm{II}}\right) .\]

(a)\(c=1\)

b)\(c=1 / 10\)

Figura 17.7: Efecto de la reducción del espacio sobre la solución.

Como antes, debido a que los datos subyacentes son cuadráticos, podemos hacer coincidir exactamente la función usando el espacio completo\(B_{0}\), es decir,\(B_{0} z_{0}=g\).

Ahora, nos limitamos a afinar funciones, y encontramos la solución de mínimos cuadrados\(z^{*}\) para\(B_{\mathrm{I}} z^{*}=g\). Nos gustaría cuantificar la diferencia en los dos primeros coeficientes del modelo completo\(z_{I}\) y los coeficientes del modelo reducido\(z^{*}\).

La figura\(17.7\) muestra el resultado de ajustar una función afín a la función cuadrática para\(c=1\) y\(c=1 / 10\). Para el\(c=1\) caso, con la fuerte dependencia cuadrática, el efecto de la función cuadrática faltante es\[\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|=1.491\] Esto da como resultado un error de solución relativamente grande de También\[\left\|z^{*}-z_{\mathrm{I}}\right\|=0.406\] observamos que, con el valor singular mínimo de\(\nu_{\min }\left(B_{\mathrm{I}}\right)=1.323\), el límite de error (absoluto) es satisfecho como\[0.406=\left\|z^{*}-z_{\mathrm{I}}\right\| \leq \frac{1}{\nu_{\min }\left(B_{\mathrm{I}}\right)}\left\|B_{\mathrm{II}} z_{I I}\right\|=1.1267\] De hecho, el encuadernado en este caso particular es razonable agudo.

Recordemos que la solución de mínimos cuadrados\(z^{*}\) minimiza el\(\ell_{2}\) residual\[0.286=\left\|B_{\mathrm{I}} z^{*}-g\right\| \leq\left\|B_{\mathrm{I}} z-g\right\|, \quad \forall z \in \mathbb{R}^{2}\] y el residual es en particular menor que el de la solución truncada\[\left\|B_{\mathrm{I}} z_{\mathrm{I}}-g\right\|=1.491\] Sin embargo, el error para la solución de mínimos cuadrados, en términos de predecir los dos primeros coeficientes del polinomio subyacente - es mayor que el de la solución truncada (que por supuesto es cero). Este caso demuestra que minimizar el residuo no necesariamente minimiza el error. Para los\(c=1 / 10\) con una dependencia cuadrática más débil, el efecto de faltar la función cuadrática es\[\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|=0.149\] y el error en la solución es en consecuencia menor ya que\[\left\|z^{*}-z_{\mathrm{I}}\right\|=0.041 .\] Esto concuerda con nuestra intuición. Si los datos subyacentes exhiben una dependencia cuadrática débil, entonces podemos representar bien los datos usando una función afín, es decir,\(\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|\) es pequeña. Entonces, el límite de error (absoluto) sugiere que el pequeño residuo da como resultado un pequeño error.

Ahora probamos el límite de error.

Prueba. Reorganizamos el sistema original como\[B_{0} z_{0}=B_{\mathrm{I}} z_{\mathrm{I}}+B_{\mathrm{II}} z_{\mathrm{II}}=g \quad \Rightarrow \quad B_{\mathrm{I}} z_{\mathrm{I}}=g-B_{\mathrm{II}} z_{\mathrm{II}} .\] Por nuestra suposición, hay una solución\(z_{\mathrm{I}}\) que satisface el sistema\(m \times(p+1)\) sobredeterminado\[B_{\mathrm{I}} z_{\mathrm{I}}=g-B_{\mathrm{II}} z_{\mathrm{II}} .\] El sistema reducido,\[B_{\mathrm{I}} z^{*}=g,\] no tiene una solución en general, por lo que se resuelve en el sentido de mínimos cuadrados. Estos dos casos son idénticos a los casos no perturbados y perturbados del lado derecho considerados el subapartado anterior. En particular, la perturbación en el lado derecho es\[\left\|g-\left(g-B_{\mathrm{II}} z_{\mathrm{II}}\right)\right\|=\left\|B_{\mathrm{II}} z_{\mathrm{II}}\right\|,\] y la perturbación en la solución es\(\left\|z^{*}-z_{\mathrm{I}}\right\|\). La sustitución de las perturbaciones en los límites de error absoluto y relativo establecidos en la subsección anterior arroja los resultados deseados.

Material Avanzado