python - 当一个变量的值已知时，在 numpy 中解决超定系统

Question

我正在尝试使用 numpy.solve 在 Python 中求解一个超定系统功能。我知道其中一个变量的值，并且我知道如果我能以某种方式插入该已知值，理论上我可以为系统找到一个独特的解决方案。

我的系统的形式是 AxC=B .变量分为两组，一组 N 个变量和一组 T 个变量(尽管这对数学无关紧要)。 A 是 (T*N x T+N)矩阵，C是变量向量，长度为(T+N) , B 是长度为 (T*N) 的向量.

我怎么知道numpy.solve (或 Python 中的另一个函数，但请不要推荐最小二乘法，我需要唯一、精确的解决方案，我知道它存在)使用其中一个变量的已知值？

我的系统的一个简单示例是:

|1 0 0 1 0|     |n1|     |B1|
|1 0 0 0 1|     |n2|     |B2|
|0 1 0 1 0|  X  |n3|  =  |B3|
|0 1 0 0 1|     |t1|     |B4|
|0 0 1 1 0|     |t2|     |B5|
|0 0 1 0 1|              |B6|  

B 的元素的值当然是已知的，还有其中一个变量的值，假设我知道 t1=1 .这些点没有任何意义，我只是把它们放在那里，这样角色就不会聚在一起。

Answer 1

正如@Foon 指出的那样，执行此操作的规范方法是减去一列。

但是，在旁注中，由于您的问题被过度确定，您必须使用最小二乘法等方法。根据定义，如果它是一个多定问题，则不存在“唯一、精确的解决方案”。 (否则它将是偶数确定的 - 一个方阵。)

除此之外，这是您的处理方式:

让我们以您的示例等式为例:

|1 0 0 1 0|     |n1|     |B1|
|1 0 0 0 1|     |n2|     |B2|
|0 1 0 1 0|  X  |n3|  =  |B3|
|0 1 0 0 1|     |t1|     |B4|
|0 0 1 1 0|     |t2|     |B5|
|0 0 1 0 1|              |B6|  

正如您所指出的，这是过度确定的。如果我们知道我们的“模型”变量之一(在这种情况下我们说 n1)，它会更加多定。这不是问题，但这意味着我们需要使用最小二乘法，并且没有完全唯一的解决方案。

那么，假设我们知道 n1 应该是什么。

在那种情况下，我们将通过从我们的观察向量中减去 n1 乘以解矩阵中的第一列来重述问题(这是@Foon 的建议):

|0 0 1 0|     |n2|     |B1 - n1|
|0 0 0 1|     |n3|     |B2 - n1|
|1 0 1 0|  X  |t1|  =  |B3 - 0 |
|1 0 0 1|     |t2|     |B4 - 0 |
|0 1 1 0|              |B5 - 0 |
|0 1 0 1|              |B6 - 0 | 

让我们用 numpy 术语使用一个更具体的例子。让我们求解方程 y = Ax^2 + Bx + C。首先，让我们生成数据和“真实”模型参数:

import numpy as np

# Randomly generate two of our model variables
a, c = np.random.rand(2)
b = 1
x = np.linspace(0, 2, 6)

y = a * x**2 + b * x + c
noise = np.random.normal(0, 0.1, y.size)
y += noise

首先，我们将在不了解 B = 1 的情况下解决它。我们可以为此使用 np.polyfit，但为了进入下一步，我们将使用较低级别的方法:

# I'm a geophysist, so I tend to use Gm=d instead of Ax=b
G = np.column_stack([x**2, x, np.ones_like(x)])
d = y

m, _, _, _ = np.linalg.lstsq(G, d)

print "Ideally, this would be 1: ", m[1]

如您所见，我们会得到接近但不完全是 1 的值。在这种情况下(我没有设置种子，所以这会有所不同)，返回的模型参数是

[ 0.13392633,  0.97217035,  0.33645734]

而真正的参数是:

[ 0.14592752,  1.        ,  0.31349185]

现在让我们考虑一下我们知道 b 的事实。我们将创建一个少一列的新 G 并从我们的观察值中减去该列时间 b (d/y):

G = np.column_stack([x**2, np.ones_like(x)])
d = y - b * x
m, _, _, _ = np.linalg.lstsq(G, d)

现在 m 是 [a, c] 并且我们已经使用我们对 b 的知识求解了这两个变量。