Cálculo numérico $(X^tX)^{-1}X^t y$ y $(X^tX)^{-1}z$

Question

En algún algoritmo necesito calcular $a=(X^tX)^{-1}X^ty$ y $b=(X^tX)^{-1} z$ en cada paso, donde $X$ es un $n \times p$ matriz no cuadrada ( $n \geq p$ , $p$ se incrementa en cada paso) y $y,z$ son algunos vectores apropiados, pero me gustaría hacerlo de manera eficiente.

Lo que se me ocurrió hasta ahora:

1. Para calcular sólo $a$ cabe destacar que $a$ es una solución de mínimos cuadrados, es decir $a$ minimiza $\Vert Xa-y\Vert_2^2$ , por lo que habría calculado $a$ utilizando la descomposición QR de $X$ (el coste es $O(np^2)$ pero ahora también tengo que calcular $b$ .
2. Así que otro enfoque ingenuo sería calcular $W := X^tX$ y $v = X^ty$ y, a continuación, calcular $a= W^{-1} v$ , $b= W^{-1} z$ . Este enfoque tiene la ventaja de que podríamos utilizar, por ejemplo, una descomposición Cholesky de $W$ (coste $O(p^3)$ ), pero los inconvenientes son que tenemos que calcular $W$ primero (coste $O(np^2)$ ) y que $X^tX$ tiene un número de condición peor que $X$ .

Pero, ¿existe una forma más eficiente de calcular $a$ y $b$ ?

Answer 1

Con $X=QR$ en la variante pequeña donde $R$ es cuadrado, se obtiene $X^tX=R^tR$ para que $a=R^{-1}Q^ty$ y $b=R^{-1}(R^t)^{-1}z$ donde sólo tienes sistemas triangulares que resolver.