用Python中的QR分解解决超定系统

Question

我正在尝试使用 QR 分解和 linalg.solve 解决超定系统，但我得到的错误是

LinAlgError：数组的最后 2 个维度必须是正方形。

当 R 数组不是正方形时就会发生这种情况，对吧？代码看起来像这样

import numpy as np
import math as ma

A = np.random.rand(2,3)
b = np.random.rand(2,1) 
Q, R = np.linalg.qr(A)
Qb = np.matmul(Q.T,b)
x_qr = np.linalg.solve(R,Qb)

对于任意 A 维度，有没有一种更有效的方法来编写它？如果没有，我该如何让这个代码片段工作？

Answer 1

原因确实是矩阵

R

不是方阵，可能是因为系统是超定的。您可以尝试

np.linalg.lstsq

，找到最小化平方误差的解决方案（如果存在，应该产生精确的解决方案）。

import numpy as np

A = np.random.rand(2, 3)
b = np.random.rand(2, 1) 
x_qr = np.linalg.lstsq(A, b)[0]

Answer 2

您需要使用标志 mode='reduced' 来调用 QR。默认 Q R 矩阵返回为 M x M 和 M x N，因此如果 M 大于 N，则矩阵 R 将是非方矩阵。如果您选择简化（经济）模式，您的矩阵将是 M x N 和 N x N，在这种情况下，求解例程将正常工作。

但是，对于超定系统，您也有向后的方程/未知数。你的代码片段应该是

import numpy as np 

A = np.random.rand(3,2)
b = np.random.rand(3,1) 
Q, R = np.linalg.qr(A, mode='reduced')
#print(Q.shape, R.shape)
Qb = np.matmul(Q.T,b)
x_qr = np.linalg.solve(R,Qb)

正如其他贡献者所指出的，您也可以直接调用 lstsq，但有时直接使用 Q 和 R 更方便（例如，如果您还计划计算投影矩阵）。

Answer 3

如

numpy.linalg.solve

的文档所示：

计算已确定的线性矩阵方程 ax = b 的“精确”解 x。

你的方程组是欠定的不是超定的。请注意，其中有 3 个变量和 2 个方程，因此方程比未知数少。

另请注意，它还提到在

numpy.linalg.solve(a,b)

中，

a

必须是

MxM

矩阵。这背后的原因是，求解方程组

Ax=b

需要计算

A

的逆矩阵，并且只有方阵是可逆的。

在这些情况下，常见的方法是采用 Moore-Penrose 伪逆，它将计算系统的最佳拟合（最小二乘）解。因此，不要尝试求解精确的解决方案，而是使用

numpy.linalg.lstsq

:

x_qr = np.linalg.lstsq(R,Qb)

Answer 4

对于解决同一非方矩阵的多个右侧的情况的更新解决方案：

import numpy as np
import scipy as sp

numRows = 500
numCols = 100
numIn   = 1000 #<! Number of inputs

mA = np.random.randn(numRows, numCols)
mB = np.random.randn(numRows, numIn)

mX = np.zeros(shape = (numCols, numIn))

mQ, mR = sp.linalg.qr(mA, mode = 'economic')

for ii in range(numIn):
    mX[:, ii] = sp.linalg.solve_triangular(mR, mQ.T @ mB[:, ii], check_finite = False)

这应该比调用更高级别的函数（如

solve()

或

lstsq()

）更快。