$\newcommand{\pinv}{A^+}$ $\newcommand{\inv}{A^{-1}}$
Para terminar, y de manera informal, si $A$ admite inversa, entonces $\pinv = \inv$. Cuando $A$ no tiene inversa, entonces $\pinv$ hace el "mejor trabajo" en invertir $A$, en el sentido de encontrar la proyección del vector $\vec{b}$ al espacio generado por las columnas de A.
$\pinv = V \Sigma^+ U^T$, donde $\Sigma^+$ tiene $1/\sigma_i$ en el espacio dado por el rango de $A$ y "$0$" en el resto:
\begin{equation} \Sigma^+ = \begin{bmatrix} 1/\sigma_1 & & & & \\ & \ddots & & & & \\ & & 1/\sigma_r & & \\ & & & 0 & \\ & & & & \ddots & \\ & & & & & 0 \\ \end{bmatrix} \end{equation}$\newcommand{\vec}[1]{\mathbf{#1}}$
Es el caso que más veces observaremos en métodos potenciales. Aquí $A$ ($n\times m$) con $n>m$. Tenemos más observaciones ($n$) que incógnitas ($m$). En este caso no existen soluciones para $\vec{x}$ dado un $\vec{b}$ genérico, salvo que $\vec{b}\in\text{Col}(A)$, es decir cuando el vector $\vec{b}$ es combinación lineal de los vectores columna en $A$.
Tomemos el siguiente ejemplo ilustrativo donde utilizamos la pseudoinversa para que arroje una "solución" al sistema $A\vec{x}=\vec{b}$,
con $A=\begin{bmatrix}a_1 \\ a_2 \\ \vdots \\ a_n \end{bmatrix}$, $\vec{x}=[x]$, y $\vec{b}=\begin{bmatrix}b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}$.
$A^+ = V\Sigma^{-1}U^T$. En este caso SVD para $A$ es: $U=\vec{a}/||\vec{a}||_2$, $\Sigma=[||a||_2]$ y $V^T=[1]$. Por lo tanto, la "solución" al problema viene expresada por
$$x = A^+ \vec{b} = (V \Sigma^{-1} U^T)\vec{b} = 1\dfrac{1}{||\vec{a}||_2}\dfrac{\vec{a}^T}{||\vec{a}||_2} \vec{b}= \dfrac{\vec{a}^T\vec{b}}{\phantom{+}||\vec{a}||_2^2}.$$La "solución" arrojada por la SVD nos dice que es esencialmente la proyección de $\vec{b}$ sobre $\vec{a}$.
Analicemos el caso donde $A=[\vec{a}]=\begin{bmatrix}1\\1\\1\\1\\1\end{bmatrix}$, $\vec{b}=\begin{bmatrix}2\\3\\4\\3\\2\end{bmatrix}$ y $\vec{x}=[c]$. Queremos ajustar una constante $c$ a una serie de puntos u observaciones. ¿Qué resultado arroja la pseudoinversa?
$$c = \dfrac{\vec{a}^T\vec{b}}{||\vec{a}||_2^2}= \dfrac{2+3+4+3+2}{1+1+1+1+1}=14/5.$$Si observamos el resultado, vemos que $c=\dfrac{1}{5}\Sigma_{j=1}^{5} b_j = <\vec{b}>$ es el promedio de las entradas de $\vec{b}$.
Claramente $A=[14/5]=\begin{bmatrix}14/5\\14/5\\14/5\\14/5\\14/5\end{bmatrix}$, no es $\vec{b}=\begin{bmatrix}2\\3\\4\\3\\2\end{bmatrix}$.
Tomando la misma $A$ y $\vec{x}$ que el ejemplo anterior, construimos $\vec{b} \in \text{Col}(A)$ para forzar una solución al sistema sobredeterminado. Consideramos $\vec{b}=\begin{bmatrix}7\\7\\7\\7\\7\end{bmatrix}$. La solución por medio de la pseudoinversa arroja $c=7$, la solución al sistema.
Basado en video sobre sistemas de ecuaciones
A continuación vemos como podemos implementar este ejemplo en Python:
import numpy as np
A = ...
b = ...
U,S,VT = np.linalg.svd(A,full_matrices=False) # descomposición SVD
c = (VT.T @ np.linalg.inv(np.diag(S)) @ U.T) @ b # resultado
En la práctica de métodos potenciales podemos utilizar la pseudoinversa para obtener la densidad de un rasgo topográfico por medio de la anomalías de aire libre, obtener superficies de tendencia, resolver el problema de la capa equivalente, etc.
Promedio de los valores en b: c = 2.8 Utilizando la expresión analítica dada por la pseudoinversa: c = [[2.8]] Utilizando la pseudoinversa numéricamente: c = [[2.8]]
Volvemos al ejemplo resuelto en otra notebook por gradientes conjugados (CG):
\begin{equation} A\mathbf{x}=\mathbf{b} \end{equation}\begin{equation} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \end{equation}con minimizador $\mathbf{x}^*=(0,1)^T$.
La solución por medio de la pseudoinversa obtiene el valor esperado $\mathbf{x}=(0,1)^T$. En Python escribimos:
A = np.array([2, 1, 1,2],dtype=float).reshape(2,2)
b = np.array([1,2],dtype=float).reshape(2,1)
x = np.linalg.pinv(A)@b # numpy.linalg.pinv obtiene la pseudoinversa
x = [[0. 1.]]^T
Eso es todo por hoy.