计算科学 - 高维欠定线性方程的数值解 - 吾爱随笔录

我需要解决线性回归问题

A x = y

$Ax=y$ 这是非常不确定的。我身边有

10^{6}

$10^6$ 功能，但只有

10^{3}

$10^3$ 方程。所以

A

$A$ 是一个

1, 000 \times 1, 000, 000

$1,000\times 1,000,000$ 矩阵和

y

$y$ 长度向量

1, 000

$1,000$ ，两者都给出，我需要找到

x

$x$ （长度

1, 000, 000

$1,000,000$ ）。

解决方案当然是 $x^*=A^\dagger y$ 在哪里 $A^\dagger$ 是伪逆。 $x^*$ 最小化最小二乘表达式 $(Ax-y)^T(Ax-y)$ .

所有这一切都是众所周知的，我可以使用标准库（例如）非常有效地做到这lstsq一点numpy。问题是我想添加形式的正则化项

min_{x} {(A x - y)^{T} (A x - y) + x^{T} R x}

$\min_x\Big\{ (Ax-y)^T(Ax-y)+x^T R x\Big\}$ 在哪里

R

$R$ 是一些

10^{6} \times 10^{6}

$10^6\times 10^6$ 我拥有的矩阵（它当然非常稀疏）。这个方程的解析解是

x = {(A^{T} A + R)}^{†} y,

$x=\left(A^TA+R\right)^\dagger y\ ,$ 但即使实例化矩阵也是完全不切实际的

A^{T} A + R

$A^TA+R$ （一个

10^{6} \times 10^{6}

$10^6\times10^6$ 矩阵）。以数值稳定的方式获得解的最佳方法是什么？如果这有帮助，我知道如何表达

R

$R$ 作为

R = B^{T} B

$R=B^TB$ 其中 b 是一个形状相同的矩阵

A

$A$ .

--- Preemtive apology:: 我意识到这可能是一个非常基本的问题，但我在我使用的任何包（scikit 等）中都找不到标准的方法。我也知道当 $R$ 是恒等式这是通常的 Tikhonov 正则化，并且有一个评估技巧 $A A^T$ 代替 $A^TA$ ，但我不确定该技巧是否适用于任意 $R$ .