python中的线性代数

Question

鉴于高m by n矩阵X,我需要计算s = 1 + x(X.T X)^{-1} x.T.这里x是行向量并且s是标量.有没有一种有效的(或推荐的)方法来在python中计算它？

不用说,X.T X将是对称肯定的.

我的尝试:

如果我们考虑QR分解X,即,X = QR哪里Q是正交,R则是上三角形,则X.T X= R.T R.

QR可使用可以容易地得到分解numpy.linalg.qr,即

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

但话又说回来,是否有一种特别有效的计算方法inv(R.T R)？

Answer 1

如果您对 进行QR因式分解X，得到X.T X = R.T R，则可以通过使用前向和后向替换（是下三角！）来避免使用np.linalg.inv(和) ：np.linalg.solveR.Tscipy.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

其中@是python3矩阵乘法运算符。