Given a tall m×n matrix X, I need to calculate s=1+x(XtX)−1xt. Here, x is a row vector. Is there an efficient (or, recommended) way to compute this in python?
Needless to say, XtX will be symmetric positive definite.
If we consider the QR decomposition of X, i.e., X=QR, where Q is orthogonal, R is upper triangular, then XtX=RtR.
QR decomposition can be easily obtained using
Q, R = numpy.linalg.qr(X)
But then again, is there a particularly efficient way to calculate (RtR)−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 XTX=ei for i=1,…,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 XTX, 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.