scipy中的线性最小二乘 - QR分解与其他方法的准确性

Question

我尝试使用以下方法在 scipy 中解决线性最小二乘问题 Ax = b：

x = numpy.linalg.inv(A.T.dot(A)).dot(A.T).dot(b) #Usually not recommended

和

x = numpy.linalg.lstsq(A, b)

两者都给出几乎相同的结果。我还尝试使用 QR 算法手动执行此操作，即：

Qmat, Rmat = la.qr(A)

bpr = dot(Qmat.T,b)
n=len(bpr)
x = np.zeros(n)
for i in xrange(n-1, -1,-1):
    x[i] = bpr[i]
    for j in xrange(i+1, n):
        x[i] -= Rmat[i, j]*x[j]
    x[i] /= Rmat[i,i]

然而，这种方法给出的结果非常不准确（大约 1e-2 的错误）。我在代码或数学上犯了 n00b 错误吗？或者，是方法的问题，还是 scipy 本身的问题？

我的 numpy 版本是 1.6.1（来自 http://www.lfd.uci.edu/~gohlke/pythonlibs/ 的 mkl 编译版本），在 x86_64 上使用 Python 2.7.3。

Answer 1

如果您使用这些二进制文件，则 QR 分解由英特尔 MKL 计算，并且可能是正确的。

对我来说，上面的代码在1e-12 内给出了随机矩阵的正确结果的解决方案。你用什么矩阵来测试它，你如何测量误差？

在某些情况下，最小二乘问题是病态的。例如， 对于具有较大零空间的矩阵，舍入误差会影响结果。考虑 rank-1 矩阵：

np.random.seed(1234)
v = np.random.rand(100)
A = v[:,None] * v[None,:]
b = np.random.randn(100)

x = scipy.linalg.lstsq(A, b)[0]
print(np.linalg.norm(A.dot(x) - b))
# -> 9.63612833771

# xp obtained using your above code
print(np.linalg.norm(A.dot(xp) - b))
# -> 3262.61161684

您自制的三角求解比 lstsq 中使用的更仔细编写的 LAPACK 例程更容易出现舍入误差，因此它的准确度会有所降低。

Answer 2

您也可以尝试使用截断特征分解。这意味着使用前 k 个特征值。我使用以下代码对最小二乘回归进行正则化，y=kc

其中u 是特征向量，lambda 是特征值

线性模型的核矩阵：

k=np.dot(X,X.T)

然后特征分解：

 w, v = scipy.linalg.eigh(k, eigvals=(lo, hi))

然后

temp= np.dot(np.dot(v,np.diag(1.0/w)),v.T)
c=np.dot(temp,y)

对于正则化，您应该在内核矩阵的对角线上添加一个小值（如0.001），否则您将拥有负特征值，这会阻止您的内核矩阵不是正定的。