Given a tall $m \times n$ matrix $X$, I need to calculate $s = 1 + x(X^t X)^{-1}x^t$. Here, $x$ is a row vector. Is there an efficient (or, recommended) way to compute this in python?
Needless to say, $X^t X$ will be symmetric positive definite.
My attempt:
If we consider the QR decomposition of $X$, i.e., $X = QR$, where $Q$ is orthogonal, $R$ is upper triangular, then $X^t X = R^t R.$
QR decomposition can be easily obtained using
Q, R = numpy.linalg.qr(X)
But then again, is there a particularly efficient way to calculate $(R^t R)^{-1}$?
The way I see it, any preprocessing such as QR decomposition, will make the computation substantionally more time and memory expensive, espeially for large $m$. For sure, QR decomposition works as a normalization in matrix computations and you can avoid some numerical instability problems with it, but I don't think it is of much use in this situtation. Hence, I would probably perform the matrix-matrix product $X^TX
Regarding the inverse, if you do need the matrix itself for some reason, you can either solve $n$ linear equations $X^TX=e_i$ for $i=1,\dotsc ,n$ or you can use the SVD implemented in python for amtrix $X$. That is stable, but rather costly once again. However it provides complete knowledge of the operator $X^TX$, since you obtain the spectral decomposition in this way and at the same time you also have a complete knowledge about $X$, which may come handy, depending on the algorithm.