Como o NumPy resolve mínimos quadrados para sistemas sub-determinados?

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

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

Esperaríamos que isso funcionasse apenas se X tivesse a forma (N, 5) onde N> = 5 Mas por que e como?

Recebemos de volta 5 pesos conforme o esperado, mas como esse problema é resolvido?

Não é como se tivéssemos 2 equações e 5 incógnitas?
Como numpy poderia resolver isso?
Ele deve fazer algo como interpolação para criar equações mais artificiais?

Meu entendimento é que o numpy.linalg.lstsq depende da rotina LAPACK dgelsd .

O problema é resolver:

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

Obviamente, isso não tem uma solução exclusiva para uma matriz A cuja classificação é menor que o comprimento do vetor $\mathbf{b}$ . No caso de um sistema indeterminado, dgelsdfornece uma solução $\mathbf{z}$ tal que:

• $A\mathbf{z} = \mathbf{b}$
• $\| \mathbf{z} \|_2 \leq \|\mathbf{x} \|_2$ para todos os$\mathbf{x}$ que satisfazem$A\mathbf{x} = \mathbf{b}$ . (ie$\mathbf{z}$ é a solução de norma mínima para o sistema indeterminado.

Exemplo, se o sistema for $x + y = 1$ , numpy.linalg.lstsq retornaria $x = .5, y = .5$ .

Como o dgelsd funciona?

A rotina dgelsdcalcula a decomposição de valor singular (SVD) de A.

Vou apenas esboçar a idéia por trás do uso de um SVD para resolver um sistema linear. A decomposição do valor singular é uma fatoração $U \Sigma V' = A$ onde $U$ e $V$ são matrizes ortogonais e $\Sigma$ é uma matriz diagonal em que as entradas diagonais são conhecidas como valores singulares.

A classificação efetiva da matriz $A$ será o número de valores singulares efetivamente diferentes de zero (ou seja, suficientemente diferentes de zero em relação à precisão da máquina, etc ...). Seja $S$ uma matriz diagonal dos valores singulares diferentes de zero. O SVD é assim:

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

O pseudo-inverso de $A$ é dado por:

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

Considere a solução $\mathbf{x} = A^\dagger \mathbf{b}$ . Então:

$$\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*}$$

Existem basicamente dois casos aqui:

1. O número de valores singulares diferentes de zero (ou seja, o tamanho da matriz $I$ ) é menor que o comprimento de $\mathbf{b}$ . A solução aqui não será exata; resolveremos o sistema linear no sentido dos mínimos quadrados.
2. $A\mathbf{x} - \mathbf{b} = \mathbf{0}$

$U$

Equivalência de pseudo-inverso

$A$

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

Para um sistema indeterminado, você pode mostrar que o pseudo-inverso fornece a solução de norma mínima.

$A$

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