algorithm - 联立方程的最小二乘解

我正在尝试适应从一组坐标到另一组坐标的转换。

x' = R + Px + Qy
y' = S - Qx + Py
Where P,Q,R,S are constants, P = scale*cos(rotation). Q=scale*sin(rotation)

有一个众所周知的“手工”公式，可将P，Q，R，S拟合到一组相应的点。
但是我需要对拟合度进行误差估计-所以我需要一个最小二乘解。

阅读“数字食谱”，但我很难解决如何对其中包含x和y的数据集执行此操作。

任何人都可以指出如何执行此操作的示例/教程/代码示例吗？
不用太在意这种语言。
但是-仅使用Matlab/Lapack/numpy/R的内置功能可能无济于事!

编辑:
我有一大堆旧的(x，y)新的(x，y)适合。这个问题已经确定得太多了(数据点多于未知数)，所以简单的矩阵求逆是不够的-正如我说的那样，我真的需要拟合的错误。

最佳答案

下面的代码应该可以解决问题。我对残差使用了以下公式:

residual[i] =   (computed_x[i] - actual_x[i])^2
              + (computed_y[i] - actual_y[i])^2

然后根据Wolfram MathWorld中描述的general procedure得出最小二乘公式。

我在Excel中测试了该算法，并按预期执行。我使用了十个随机点的集合，然后通过随机生成的转换矩阵对其进行旋转，平移和缩放。

在没有随机噪声施加到输出数据的情况下，此程序将产生与输入参数相同的四个参数(P，Q，R和S)，并且rSquared值为零。

随着越来越多的随机噪声被施加到输出点，这些常数开始偏离正确的值，并且rSquared值相应地增加。

这是代码:

// test data
const int N = 1000;
float oldPoints_x[N] = { ... };
float oldPoints_y[N] = { ... };
float newPoints_x[N] = { ... };
float newPoints_y[N] = { ... };

// compute various sums and sums of products
// across the entire set of test data
float Ex =  Sum(oldPoints_x, N);
float Ey =  Sum(oldPoints_y, N);
float Exn = Sum(newPoints_x, N);
float Eyn = Sum(newPoints_y, N);
float Ex2 = SumProduct(oldPoints_x, oldPoints_x, N);
float Ey2 = SumProduct(oldPoints_y, oldPoints_y, N);
float Exxn = SumProduct(oldPoints_x, newPoints_x, N);
float Exyn = SumProduct(oldPoints_x, newPoints_y, N);
float Eyxn = SumProduct(oldPoints_y, newPoints_x, N);
float Eyyn = SumProduct(oldPoints_y, newPoints_y, N);

// compute the transformation constants
// using least-squares regression
float divisor = Ex*Ex + Ey*Ey - N*(Ex2 + Ey2);
float P = (Exn*Ex + Eyn*Ey - N*(Exxn + Eyyn))/divisor;
float Q = (Exn*Ey + Eyn*Ex + N*(Exyn - Eyxn))/divisor;
float R = (Exn - P*Ex - Q*Ey)/N;
float S = (Eyn - P*Ey + Q*Ex)/N;

// compute the rSquared error value
// low values represent a good fit
float rSquared = 0;
float x;
float y;
for (int i = 0; i < N; i++)
{
    x = R + P*oldPoints_x[i] + Q*oldPoints_y[i];
    y = S - Q*oldPoints_x[i] + P*oldPoints_y[i];
    rSquared += (x - newPoints_x[i])^2;
    rSquared += (y - newPoints_y[i])^2;
}

关于algorithm - 联立方程的最小二乘解，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/1245968/