linear algebra in python

Question

Given a tall m by n matrix X, I need to calculate s = 1 + x(X.T X)^{-1} x.T. Here, x is a row vector and s is scalar. 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 numpy.linalg.qr, that is

Q,R = numpy.linalg.qr(X)

But then again, is there a particularly efficient way to calculate inv(R.T R)?

Answer 1

If you are doing the QR factorization of X, resulting in X.T X = R.T R, you may avoid using np.linalg.inv (and np.linalg.solve) by using forward and backward substitution instead (R.T is lower triangular!) with scipy.linalg.solve_triangular:

import numpy as np
import scipy.linalg as LA
Q,R = np.linalg.qr(X)
# solve R.T R z = x such that R z = y 
# with step (a) then (b)
# step (a) solve  R.T y = x
y   = LA.solve_triangular(R,x,trans='T') 
# step (b) solve R z = y
z   = LA.solve_triangular(R,x)
s   = 1 + x @ z

where @ is the python3 matrix multiplication operator.