联立方程的最小二乘解答案

【问题标题】：Least Squares solution to simultaneous equations联立方程的最小二乘解
【发布时间】：2010-11-17 18:34:05
【问题描述】：

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

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 的内置功能可能没有帮助！

编辑：我有一大套 old(x,y) new(x,y) 适合。问题是超定的（数据点多于未知数），所以简单的矩阵求逆是不够的 - 正如我所说，我真的需要拟合误差。

【问题讨论】：

你有一组 (x_i,y_i,x'_i,y'_i)s 或 x+/-dx, y+/-dy ... 还是什么？如果 x, y, x', y' 各有一个，您可以只做一个精确的解决方案，并且没有办法提取错误估计......
没有一个最小化类型函数 LMA/Gauss-Newton 给出直接错误。我想我可以计算出最佳拟合，然后计算出每个点的误差。我原以为它比这简单得多（即线性 LSquers 的简单模式），我只是愚蠢
有趣的问题！我发布了一些应该可以解决问题的代码，但有趣的部分是再次打破我生锈的旧数学技能:)

标签： algorithm language-agnostic math

【解决方案1】：

感谢 eJames，这几乎就是我所拥有的。我从一本旧的军队测量手册中对其进行编码，该手册基于早期的“测量员说明”笔记，该笔记必须有 100 年的历史！（它使用 N 和 E 表示 North 和 East 而不是 x/y ）

拟合优度参数将非常有用 - 如果选定的点使拟合更差，我可以交互地丢弃它们。

FindTransformation(vector<Point2D> known,vector<Point2D> unknown) {
{
    // sums
    for (unsigned int ii=0;ii<known.size();ii++) {
       sum_e += unknown[ii].x;
       sum_n += unknown[ii].y;
       sum_E += known[ii].x;
       sum_N += known[ii].y;                            
       ++n;         
    }

    // mean position
    me = sum_e/(double)n;
    mn = sum_n/(double)n;
    mE = sum_E/(double)n;
    mN = sum_N/(double)n;

    // differences
    for (unsigned int ii=0;ii<known.size();ii++) {

       de = unknown[ii].x - me;
       dn = unknown[ii].y - mn;

       // for P
       sum_deE += (de*known[ii].x);
       sum_dnN += (dn*known[ii].y);
       sum_dee += (de*unknown[ii].x);
       sum_dnn += (dn*unknown[ii].y);

       // for Q
       sum_dnE += (dn*known[ii].x);
       sum_deN += (de*known[ii].y);                     
   }

double P = (sum_deE + sum_dnN) / (sum_dee + sum_dnn);
double Q = (sum_dnE - sum_deN) / (sum_dee + sum_dnn);

double R = mE - (P*me) - (Q*mn);
double S = mN + (Q*me) - (P*mn);
}

【讨论】：

酷。想想他们 100 年前是如何计算的，我不寒而栗 :)

【解决方案2】：

下面的代码应该可以解决问题。我使用以下公式计算残差：

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;
}

【讨论】：

我通常会使用更具描述性的变量名称，但在这种情况下，我认为像sumOfNewXValues 这样的长名称实际上会使所有内容更难阅读。数学公式似乎是一个特例。

【解决方案3】：

定义 3x3 矩阵 T(P,Q,R,S) 使得(x',y',1) = T (x,y,1)。然后计算

A = \sum_i |(T (x_i,y_i,1)) - (x'_i,y'_i,1)|^2

并针对 (P,Q,R,S) 最小化 A。

自己编写代码是一个大中型项目，除非您可以保证数据条件良好，尤其是当您希望从程序中获得良好的错误估计时。您最好使用支持误差估计的现有最小化器。

粒子物理类型将使用minuit 直接来自CERNLIB（最容易在fortran77 中完成编码），或来自ROOT（使用c++ 编码，或者应该可以通过python 绑定访问） .但是，如果您还没有这些工具中的任何一个，那就是一个很大的安装。

我相信其他人可以推荐其他最小化器。

【讨论】：

【解决方案4】：

您可以使用levmar 程序进行计算。它经过测试并集成到包括我的在内的多种产品中。它在 GPL 下获得许可，但如果这是一个非开源项目，他将为您更改许可（收费）

【讨论】：

谢谢，有很多免费的 LMA 实现。我没有意识到这是一种看起来很简单的必要方法。
它的误差估计很难。这需要解决方案矩阵（IIRC）的雅可比行列式。请记住，LS 错误不是世界考虑错误的意义上的错误。它们是答案的数值稳定性的量度。所以，有些东西可能是正确的解决方案，但不是特别稳定（即稍微改变一个值会导致目标函数发生很大变化）。
是的，我认为从用户的角度来看，解决误差的最佳方法是计算最佳位置，然后计算每个点的残差，然后取它们的标准差。
@mgb - 请注意，成熟的用户不会期待（并且可能不想要它）。 LM 错误通常由数学包报告

【解决方案5】：

要找到 P、Q、R 和 S，您可以使用最小二乘法。我认为令人困惑的是，最小二乘的通常描述使用 x 和 y，但它们与您的问题中的 x 和 y 不匹配。您只需要将您的问题仔细转换为最小二乘框架。在您的情况下，自变量是未变换的坐标 x 和 y，因变量是变换后的坐标 x' 和 y'，可调参数是 P、Q、R 和 S。（如果这还不够清楚，让我知道，我会发布更多详细信息。）

一旦你找到了 P、Q、R 和 S，然后 scale = sqrt(P^2 + Q^2) 然后你可以从 sin rotation = Q / scale 和 cos rotation = P / 找到旋转规模。

【讨论】：

是的，线性最小二乘法的通常表示是为了拟合标量方程 y = F(x) = \sum_i c_i f_i(x);这个问题是拟合向量方程 r' = F(r) = \sum_i c_i f_i(r)。
问题是线性LSF只有一个变量+2个未知数。我有 2 个变量和 4 个未知数（如果您假设比例相同，则为三个）
线性最小二乘法（在通常的表示中）实际上适用于一个变量和 n 个未知数。请参阅en.wikipedia.org/wiki/Linear_least_squares，其中第一个图形示例用于拟合二阶多项式 (n = 3)。
最小二乘可以处理 m > 1 个自变量和 n > 1 个因变量。在我的数字食谱版本中，它在第 15.4 节一般线性最小二乘中有所介绍，或者查看上面引用的维基百科文章。
谢谢，我认为这一章现在更有意义了！当我之前这样做时，总是有一些技巧可以将其转换为线性方程。手动解决方案只是一个简单的差异平方和过程，我认为它很简单，但我遗漏了一些东西。

【解决方案6】：

一个问题是，像这样的数字内容通常很棘手。即使算法很简单，在实际计算中也经常会出现问题。

因此，如果您可以轻松获得具有内置功能的系统，那么最好使用它。

【讨论】：