¿Cómo resuelve SVD el problema de mínimos cuadrados?

Question

¿Cómo puedo demostrar que la solución de mínimos cuadrados para $$\text{minimizar} \quad \|Ax-b\|_2$$ es $A^{+} b$, donde $A^{+}$ es la pseudoinversa de $A$?

Answer 1

La pseudoinversa de Moore-Penrose es una consecuencia natural de aplicar la descomposición de valores singulares al problema de mínimos cuadrados. La SVD resuelve el problema de mínimos cuadrados en dos componentes: (1) una parte del espacio de rango que se puede minimizar y (2) un término de espacio nulo que no se puede eliminar, es decir, un error residual. La primera parte naturalmente creará la solución seudoinversa.

Definir SVD

Comenzar con una matriz no nula $\mathbf{A}\in\mathbb{C}^{m\times n}_{\rho}$, donde el rango de la matriz $1\le\rho y$\rho. La descomposición de valores singulares, garantizada de existir, es $$\mathbf{A} = \mathbf{U} \, \Sigma \, \mathbf{V}^{*} = \left[ \begin{array}{cc} \color{blue}{\mathbf{U}_{\mathcal{R}}} & \color{red}{\mathbf{U}_{\mathcal{N}}} \end{array} \right] % \left[ \begin{array}{c} \mathbf{S} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{array} \right] % \left[ \begin{array}{c} \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \\ \color{red}{\mathbf{V}_{\mathcal{N}}}^{*} \end{array} \right].$$ La matriz codominio $\mathbf{U}\in\mathbb{C}^{m\times m}$ y la matriz de dominio $\mathbf{V}\in\mathbb{C}^{n\times n}$ son unitarias: $$\mathbf{U}^{*}\mathbf{U} = \mathbf{U}\mathbf{U}^{*} = \mathbf{I}_{m}, \quad \mathbf{V}^{*}\mathbf{V} = \mathbf{V}\mathbf{V}^{*} = \mathbf{I}_{n}.$$ Los vectores columna de las matrices del dominio proporcionan bases ortonormales para los cuatro subespacios fundamentales: $$\begin{array}{ll} % matriz & subespacio \\\hline % \color{blue}{\mathbf{U}_{\mathcal{R}}}\in\mathbb{C}^{m\times\rho} & \color{blue}{\mathcal{R}\left(\mathbf{A}\right)} = \text{span}\left\{\color{blue}{u_{1}},\dots,\color{blue}{u_{\rho}}\right\}\\ % \color{blue}{\mathbf{V}_{\mathcal{R}}}\in\mathbb{C}^{n\times\rho} & \color{blue}{\mathcal{R}\left(\mathbf{A}^{*}\right)} = \text{span}\left\{\color{blue}{v_{1}},\dots,\color{blue}{v_{\rho}}\right\}\\ % \color{red}{\mathbf{U}_{\mathcal{N}}}\in\mathbb{C}^{m\times m-\rho} & \color{red}{\mathcal{N}\left(\mathbf{A^{*}}\right)} = \text{span}\left\{\color{red}{u_{\rho+1}},\dots,\color{red}{u_{m}}\right\}\\ % \color{red}{\mathbf{V}_{\mathcal{N}}}\in\mathbb{C}^{n\times n-\rho} & \color{red}{\mathcal{N}\left(\mathbf{A}\right)} = \text{span}\left\{\color{red}{v_{\rho+1}},\dots,\color{red}{v_{n}}\right\}\\ % \end{array}$$ Hay $\rho$ valores singulares que están ordenados y son reales: $$\sigma_{1} \ge \sigma_{2} \ge \dots \ge \sigma_{\rho}>0,$$ y son la raíz cuadrada de los autovalores no nulos de las matrices de productos $\mathbf{A}^{*}\mathbf{A}$ y $\mathbf{A}\mathbf{A}^{*}$. Estos valores singulares forman la matriz diagonal de valores singulares $$\mathbf{S} = \text{diagonal} (\sigma_{1},\sigma_{2},\dots,\sigma_{\rho}) = \left[ \begin{array}{ccc} \sigma_{1} \\ & \ddots \\ && \sigma_{\rho} \end{array} \right] \in\mathbb{R}^{\rho\times\rho}.$$ La matriz $\mathbf{S}$ está incrustada en la matriz sabot $\Sigma\in\mathbb{R}^{m\times n}$ cuya forma garantiza la conformabilidad.

Definir solución de mínimos cuadrados

Considerar que el sistema lineal $$\mathbf{A} x = b$$ no tiene una solución exacta, por lo que generalizamos la pregunta y pedimos la mejor solución vectorial $x$. Utilizando la norma $2-$, pedimos la solución de mínimos cuadrados que minimiza $r^{2}(x) = \lVert \mathbf{A} x - b \rVert_{2}^{2}$, la suma de los cuadrados de los errores residuales: $$x_{LS} = \left\{ x \in \mathbb{C}^{n} \colon \big\lVert \mathbf{A} x - b \big\rVert_{2}^{2} \text{ se minimiza} \right\}$$

Explotar SVD - resolver componentes de espacio rango y nulo

Una propiedad útil de las transformaciones unitarias es que son invariantes bajo la norma $2-$. Por ejemplo $$\lVert \mathbf{V} x \rVert_{2} = \lVert x \rVert_{2}.$$ Esto proporciona una libertad para transformar problemas en una forma más fácil de manipular. En este caso, $$\begin{align} r\cdot r &= \lVert \mathbf{A}x - b \rVert_{2}^{2} = \lVert \mathbf{U}^{*}\left(\mathbf{A}x - b \right) \rVert_{2}^{2} = \lVert \mathbf{U}^{*}\left(\mathbf{U} \, \Sigma \, \mathbf{V}^{*}x - b \right) \rVert_{2}^{2} \\ &= \lVert \Sigma \mathbf{V}^{*} x - \mathbf{U}^{*} b \rVert_{2}^{2}. \end{align}$$ Al cambiar a la forma en bloque se separan los términos de espacio rango y espacio nulo, $$\begin{align} r^{2}(x) &= \big\lVert \Sigma \mathbf{V}^{*} x - \mathbf{U}^{*} b \big\rVert_{2}^{2} = \Bigg\lVert % \left[ \begin{array}{c} \mathbf{S} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{array} \right] % \left[ \begin{array}{c} \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} \\ \color{red}{\mathbf{V}_{\mathcal{N}}}^{*} \end{array} \right] % x - \left[ \begin{array}{c} \color{blue}{\mathbf{U}_{\mathcal{R}}}^{*} \\ \color{red}{\mathbf{U}_{\mathcal{N}}}^{*} \end{array} \right] b \Bigg\rVert_{2}^{2} \\ &= \big\lVert \mathbf{S} \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} x - \color{blue}{\mathbf{U}_{\mathcal{R}}}^{*} b \big\rVert_{2}^{2} + \big\lVert \color{red}{\mathbf{U}_{\mathcal{N}}}^{*} b \big\rVert_{2}^{2} \end{align}$$

La separación entre las contribuciones de espacio rango y espacio nulo al error total es también una separación entre componentes que están controladas y no controladas. El vector que controlamos es el vector solución $x$ que aparece solo en el término de espacio rango (azul). Lo que queda es el término de espacio nulo (rojo), un error residual.

Resolver y recuperar la pseudoinversa de Moore-Penrose

Seleccionar el vector $x$ en el término de espacio rango que fuerce a ese término a ser $0$: $$\mathbf{S} \color{blue}{\mathbf{V}_{\mathcal{R}}}^{*} x - \color{blue}{\mathbf{U}_{\mathcal{R}}}^{*} b = 0 \qquad \Rightarrow \qquad x = \color{blue}{\mathbf{V}_{\mathcal{R}}} \mathbf{S}^{-1} \color{blue}{\mathbf{U}_{\mathcal{R}}}^{*}b.$$ Esta es la solución de mínimos cuadrados $$\color{blue}{x_{LS}} = \color{blue}{\mathbf{A}^{\dagger} b} = \color{blue}{\mathbf{V}_{\mathcal{R}}} \mathbf{S}^{-1} \color{blue}{\mathbf{U}_{\mathcal{R}}}^{*}b$$ usando la SVD delgada.

Queda un término de error que no podemos eliminar, un error residual, dado por $$r^{2}\left(x_{LS}\right) = \big\lVert \color{red}{\mathbf{U}_{\mathcal{N}}^{*}} b \big\rVert_{2}^{2},$$ que cuantifica la porción del vector de datos $b$ que reside en el espacio nulo.

Answer 2

Primer condición de orden:

$$\frac{d}{dx}||Ax-b||_2^2=\frac{d}{dx}(Ax-b)'(Ax-b)=A'(Ax-b)=0$$

Por lo tanto,

$$x=(A'A)^{-1}A'b=A^+b$$

Answer 3

Lo primero que hay que saber es que los posibles valores de $Ax$ cubren todos y solo el rango de $A$: $\{Ax: x\in V\}=\mathrm{Range}(A)$. Esto significa que el problema $Ax=b$ es solucionable si y solo si $b\in\mathrm{Range}(A)$.

Ahora consideremos el caso más general en el que $b\in \mathrm{Range}(A)$ no es necesariamente cierto, y por lo tanto el problema no puede resolverse exactamente. Dado que $A$ nos permite generar todo en $\mathrm{Range}(A)$ mediante una elección adecuada de $x$, entonces es natural que el $x$ que brinde la mejor solución (en cuanto a la norma $\|\cdot\|_2$) sea aquel tal que $Ax$ sea igual a la proyección de $b$ en $\mathrm{Range}(A)$. En otras palabras, denotando con $\mathbb P_R$ el proyector en $\mathrm{Range}(A)$, queremos que $x$ cumpla con $$Ax=\mathbb P_R b.\tag A$$ Con esta elección de $x$, la distancia entre $Ax$ y $b$ es entonces $$\|b-\mathbb P_R b\|_2=\|(I-\mathbb P_R)b\|_2.$$

Ahora, ¿por qué entra en la discusión la descomposición en valores singulares (SVD)? Bueno, mayormente porque proporciona una forma fácil de encontrar el $x$ que satisface (A).

Para ver esto, escribamos la SVD en notación dyádica como $$A=\sum_k s_k u_k v_k^*,\tag B$$ donde $s_k>0$ son los valores singulares, y $\{u_k\}$ y $\{v_k\}$ son bases ortonormales para $\mathrm{Range}(A)$ y $A^{-1}(\mathrm{Range}(A)\setminus \{0\})=\mathrm{Ker}(A)_\perp$, respectivamente. La expresión (B) es equivalente a la forma tal vez más común de ver la SVD, $A=UDV^\dagger$, con $v_k$ siendo las columnas de $V$ y $u_k$ las columnas de $U$.

La pseudo-inversa $A^+$ tiene una buena expresión en función de su SVD: $$A^+=\sum_k s_k^{-1} v_k u_k^*. \tag C$$ Con estas expresiones, puedes notar cómo $$A(A^+ y)=\sum_k u_k u_k^* y=\sum_k u_k\langle u_k,y\rangle=\mathbb P_R y.$$ Ahora, recuerda que queremos acercarnos lo más posible a $b$, y por lo tanto estamos buscando algún $y$ tal que $$A(A^+ y)=\mathbb P_R y=\mathbb P_R b.$$ Esto nos dice que queremos que $y=b$, y por lo tanto $x=A^+ b$.

Answer 4

Considera la ecuación ${ Ax = b }$ donde ${ A \in \mathbb{R} ^{m \times n} ,}$ ${ b \in \mathbb{R} ^{m} }.$ El objetivo es encontrar una solución aproximada ${ \hat{x} \in \mathbb{R} ^{n} ,}$ aquella para la cual $${ \lVert A \hat{x} - b \rVert = \min _{x \in \mathbb{R} ^n} \lVert Ax - b \rVert}.$$ Tal ${ \hat{x} }$ no necesita ser única. Por ejemplo, cuando ${ m = 1 }$ y ${ b = 0 }.$

Tenemos dos enfoques, y el enfoque de SVD proporciona una solución explícita.

Sin SVD:

Elige cualquier base ${ \mathcal{B} = (v _1, \ldots, v _r) }$ de ${ \text{colsp}(A) }.$

Una opción natural para la base es: Supongamos que la forma escalonada reducida por filas de ${ A }$ es ${ A ^{'} = E _k \ldots E _1 A = (A ^{'} _1, \ldots, A ^{'} _{n}) ,}$ con columnas pivote en posiciones ${ i _1 \lt \ldots \lt i _r }.$ Entonces ${ (A ^{'} _{i _1}, \ldots, A ^{'} _{i _r}) }$ es una base de ${ \text{colsp}(A^{'}) }.$
Ahora ${ \mathscr{B} = (A _{i _1}, \ldots, A _{i _r}) }$ es una base de ${ \text{colsp}(A) }$. [Ya que para cualquier matriz elemental ${ E },$ ${ (M _{i _1}, \ldots, M _{i _r}) }$ es una base de ${ \text{colsp}(M) }$ si y solo si ${ (EM _{i _1}, \ldots, EM _{i _r}) }$ es una base de ${ \text{colsp}(EM) }$].

Queremos un ${ \tilde{x} \in \mathbb{R} ^{r} }$ para el cual ${ \mathcal{B} \tilde{x} }$ (un punto genérico en ${ \text{colsp}(\mathcal{B}) = \text{colsp}(A) }$) sea el más cercano a ${ b }.$ Para tal ${ \tilde{x} },$ tenemos que ${ b - \mathcal{B}\tilde{x} }$ es ortogonal a ${ \text{colsp}(\mathcal{B}) },$ es decir, ${ \mathcal{B} ^{T} (b - \mathcal{B} \tilde{x}) = 0},$ es decir, ${ \mathcal{B} ^{T} b = \mathcal{B} ^{T} \mathcal{B} \tilde{x} ,}$ es decir, ${ \tilde{x} = (\mathcal{B} ^T \mathcal{B}) ^{-1} \mathcal{B} ^{T} b }.$

Entonces diferentes elecciones de bases ${ \mathcal{B} }$ de ${ \text{colsp}(A) }$ dan diferentes puntos ${ \mathcal{B}\tilde{x} = \mathcal{B}(\mathcal{B} ^T \mathcal{B}) ^{-1} \mathcal{B} ^{T} b \in \text{colsp}(A) }$ más cercanos a ${ b }.$ (Por cierto, la invertibilidad de ${ \mathcal{B} ^{T} \mathcal{B} }$ se sigue, por ejemplo, de aquí).

Con SVD:

[Igual que la respuesta de dantopa arriba]

Considera la notación en aquí. Sustituir la SVD de ${ A }$ en ${ \lVert Ax - b \rVert ^2 }$ da $${ \begin{align*} \lVert A x - b \rVert ^{2} &= \left\lVert \mathscr{U} \begin{pmatrix} \Sigma _r &0 \\ 0 &0 \end{pmatrix} \mathscr{V} ^{*} x - b \right\rVert ^{2} \\ &= \left\lVert \mathscr{U} \left( \begin{pmatrix} \Sigma _r &0 \\ 0 &0 \end{pmatrix} \mathscr{V} ^{*} x - \mathscr{U} ^{*} b \right) \right\rVert ^{2} \\ &= \left\lVert \begin{pmatrix} \Sigma _r &0 \\ 0 &0 \end{pmatrix} \begin{pmatrix} \mathscr{V} _{r} ^{*} \\ \mathscr{V} _{n-r} ^{*} \end{pmatrix} x - \begin{pmatrix} \mathscr{U} _{r} ^{*} \\ \mathscr{U} _{m-r} ^{*} \end{pmatrix} b \right\rVert ^{2} \\ &= \left\lVert \Sigma _{r} \mathscr{V} _{r} ^{*} x - \mathscr{U} _{r} ^{*} b \right\rVert ^{2} + \left\lVert \mathscr{U} _{m-r} ^{*} b \right\rVert ^{2}. \end{align*} }$$

Los minimizadores de esta cantidad son precisamente aquellos ${ x }$ que satisfacen ${ \Sigma _{r} \mathscr{V} _{r} ^{*} x - \mathscr{U} _{r} ^{*} b = 0 }$ es decir, $${ \mathscr{V} _{r} ^{*} x = \Sigma _{r} ^{-1} \mathscr{U} _{r} ^{*} b. }$$

Nota que ${ x = \mathscr{V} _{r} \Sigma _{r} ^{-1} \mathscr{U} _{r} ^{*} b }$ es una solución particular, por lo tanto, el conjunto completo de soluciones es $${ \begin{align*} x &\in \mathscr{V} _{r} \Sigma _{r} ^{-1} \mathscr{U} _{r} ^{*} b + \mathcal{N}(\mathscr{V} _r ^{*}) \\ &= A ^{\dagger} b + \text{span}(\mathscr{V} _{n-r}). \end{align*} }$$

En particular, el conjunto de minimizadores es de dimensión ${ (n-r) }$.

Nota que ${ A ^{\dagger} b = \mathscr{V} _{r} \Sigma _{r} ^{-1} \mathscr{U} _{r} ^{*} b }$ está en ${ \text{span}(\mathscr{V} _r) }$ y por lo tanto es ortogonal a ${ \text{span}(\mathscr{V} _{n-r}) }.$ Así que entre los minimizadores ${ A ^{\dagger} b + \text{span}(\mathscr{V} _{n-r}) },$ el minimizador con la menor norma es ${ x = A ^{\dagger} b }.$