¿Cómo resuelve NumPy mínimos cuadrados para sistemas subdeterminados?

Digamos que tenemos X de forma (2, 5)
e y de forma (2,)

Esto funciona: np.linalg.lstsq(X, y)

Es de esperar que esto funcione solo si X tenía forma (N, 5) donde N> = 5 Pero, ¿por qué y cómo?

Volvemos 5 pesos como se esperaba, pero ¿cómo se resuelve este problema?

¿No es que tenemos 2 ecuaciones y 5 incógnitas?
¿Cómo podría numpy resolver esto?
Debe hacer algo como la interpolación para crear más ecuaciones artificiales? ..

Tengo entendido que numpy.linalg.lstsq se basa en la rutina LAPACK dgelsd .

El problema es resolver:

$$\text{minimize} (\text{over} \; \mathbf{x}) \quad \| A\mathbf{x} - \mathbf{b} \|_2$$

Por supuesto, esto no tiene una solución única para una matriz A cuyo rango es menor que la longitud del vector $\mathbf{b}$ . En el caso de un sistema indeterminado, dgelsdproporciona una solución $\mathbf{z}$ tal que:

• $A\mathbf{z} = \mathbf{b}$
• $\| \mathbf{z} \|_2 \leq \|\mathbf{x} \|_2$ para todas las$\mathbf{x}$ que satisfacen$A\mathbf{x} = \mathbf{b}$ . (es decir,$\mathbf{z}$ es la solución de norma mínima para el sistema indeterminado.

Ejemplo, si el sistema es $x + y = 1$ , numpy.linalg.lstsq devolvería $x = .5, y = .5$ .

¿Cómo funciona dgelsd?

La rutina dgelsdcalcula la descomposición del valor singular (SVD) de A.

Esbozaré la idea detrás de usar un SVD para resolver un sistema lineal. La descomposición del valor singular es una factorización $U \Sigma V' = A$ donde $U$ y $V$ son matrices ortogonales y $\Sigma$ es una matriz diagonal donde las entradas diagonales se conocen como valores singulares.

El rango efectivo de la matriz $A$ será el número de valores singulares que efectivamente no son cero (es decir, suficientemente diferentes de cero en relación con la precisión de la máquina, etc.). Sea $S$ una matriz diagonal de los valores singulares distintos de cero. La SVD es así:

$$A = U \begin{bmatrix} S & 0 \\ 0 & 0 \end{bmatrix} V'$$

El pseudoinverso de $A$ viene dado por:

$$A^\dagger = V \begin{bmatrix} S^{-1} & 0 \\ 0 & 0 \end{bmatrix} U'$$

Considere la solución $\mathbf{x} = A^\dagger \mathbf{b}$ . Luego:

$$\begin{align*} A\mathbf{x} - \mathbf{b} &= U \begin{bmatrix} S & 0 \\ 0 & 0 \end{bmatrix} V' V \begin{bmatrix} S^{-1} & 0 \\ 0 & 0 \end{bmatrix} U' \mathbf{b} - \mathbf{b} \\ &= U \begin{bmatrix} I & 0 \\ 0 & 0 \end{bmatrix} U' \mathbf{b} - \mathbf{b}\\ \end{align*}$$

Básicamente hay dos casos aquí:

1. $I$$\mathbf{b}$
2. $A\mathbf{x} - \mathbf{b} = \mathbf{0}$

$U$

Equivalencia de pseudoinverso

$A$

$$A^\dagger = A'\left(AA' \right)^{-1}$$

Para un sistema indeterminado, puede demostrar que el pseudo-inverso le brinda la solución de norma mínima.

Cuando tiene columnas linealmente independientes (por ejemplo, tenemos una matriz delgada), entonces: $A$

$$A^\dagger = \left(A'A \right)^{-1}A'$$