最小二乘法

Mosi_62

5.00/5 (15投票s)

2016年3月17日

CPOL

7分钟阅读

32167

1432

在超定系统中近似多项式函数。使用法方程或正交变换。

引言

最小二乘法是一种标准程序，用于将多项式函数近似到一组参考点。多项式的阶数 n 低于参考点的数量。这会导致超定方程组。例如，如果我们有一组 7 个参考点（x 和 y）并想近似 3 次多项式函数，我们可以建立 7 个如下方程：

co + c1 * x + c1 * x^2 = d

对于 3 个独立值。从而得到一个矩阵方程

Cx = d

最小二乘法寻找一个多项式函数，该函数使得每个方程的误差平方和最小。解决此问题的一种方法是使用法方程。

法方程

法方程由卡尔·弗里德里希·高斯用于通过观测点计算恒星的椭圆轨道。它将一个误差向量 r（残差向量）引入上述矩阵方程：

Cx – d = r

经过复杂的矩阵转置、前向替换和后向替换计算后，解决方案是当误差向量最小时：

这意味着上述超定矩阵方程的两边都必须乘以 C 的转置矩阵。这会得到一个 n * n 的矩阵方程，其中 n 是所搜索多项式的阶数。这个新的矩阵方程可以被求解，从而得到向量 x 中的正确结果。
乘以转置矩阵可以将超定矩阵方程转换为一个确定的系统。在我上面的例子中，有 7 个参考点和 3 次多项式，C 有 7 行 3 列。乘以其转置后，它变成一个 3x3 的矩阵，d 也变成一个 3 维向量。可以使用高斯消元法轻松求解此矩阵方程。
因此，第一步是填充矩阵方程。在 C 中，我们必须放入参考点的 x 值。C 的第一列全为 1，因为我们要找的第一个参数是 c0，乘以 1。第二列是 x1...xn，第三列是 (x1)^2...(xn)^2，以此类推。
我定义了一个 TMatrix 结构来包含 C、x 和 d，并创建了一个实例 m 来进行操作。所以

for (i = 0; i < samples; i++)
{
  m.d[i] = y_k[i];
  for (k = 0; k < order; k++)
  {
    m.c[i, k] = Math.Pow(x_k[i], k);
  }
}

填充矩阵方程

现在我们必须将方程乘以 C 的转置。C 的转置矩阵是 C 沿其主对角线翻转。在示例中，3x7 的矩阵变成 7x3 的矩阵，沿红线标记的主对角线翻转。

现在我们必须将方程乘以 C 的转置。C 的转置矩阵是 C 沿其主对角线翻转。在示例中，3x7 的矩阵变成 7x3 的矩阵，沿红线标记的主对角线翻转。

要计算 C^T 乘以 C，C^T 的一行中的值必须乘以 C 的一列中的值，然后相加。所以示例中结果矩阵中的每个值都是 7 次乘法的总和，我们得到一个 3x3 的矩阵。d 也是一样：C^T 中每一行的值乘以 d 中的值然后相加。为此，我们需要一个临时矩阵 a。

for (i = 0; i < order; i++)
{
  a.d[i] = 0.0;
  for (k = 0; k < samples; k++)
  a.d[i] = a.d[i] + m.c[k, i] * m.d[k];
  for (j = 0; j < order; j++)
  {
    a.c[j,i] = 0.0;
    for (k = 0; k < samples; k++)
    {
      a.c[j, i] = a.c[j, i] + m.c[k,j] * m.c[k, i];
    }
  }
}

计算 transposed(C) * C 和 transposed(C) * d，并将结果放入 a。矩阵方程现在可以在 a 中通过高斯法求解，我们得到 a.x[1] .. a.x[n] 中的结果。

使用正交变换

文献中提到，使用法方程可能会导致不准确的结果。这主要是由于乘以 A * A^T 可能会导致结果矩阵的条件变差，从而在后续运算中产生舍入误差。好吧，我目前还没有遇到这种情况，但无论如何：使用正交变换是一个有趣的方法 :-)

使用正交变换基本上是使用 Givens 的旋转矩阵（参见 http://www.mosismath.com/RotationMatrix/RotationMatrix.html），基本思想是向量乘以正交矩阵不会改变向量的长度，因此矩阵方程与残差向量

Cx – d = r

可以通过这些正交矩阵进行变换，而不会影响残差平方和。

如果我们扩展主矩阵 C，它是一个非方阵 n*m，使其成为一个方阵，并将新元素全部设为 0，如下所示：

整个方程变为：

C’x – d = r

可以通过旋转矩阵进行乘法运算，将 C' 转换为上三角矩阵，而不会影响残差平方和。这一点很重要。

变换是迭代进行的。从元素 c₂₁ 开始，构建消除 C' 中该元素的旋转矩阵，进行乘法运算，然后构建下一个旋转矩阵以消除 C' 中下一个元素 c₃₁，以此类推，直到我们得到上三角矩阵，如下所示：

这在一个函数中完成

public void RunGivens(int maxOrder)
{
      int i, k;
      for (i = 0; i < maxOrder; i++)
      {
            for (k = i + 1; k < maxOrder; k++)
            {
                  if (m.c[k, i] != 0.0)
                  {
                        CalcRMatrix(i, k);
                        m.MultMatrix(r, i, k);
                  }
            }
      }
}

用

public void CalcRMatrix(int p, int q)
{
      int i, j;
      double radius;
      radius = Math.Sqrt(m.c[p, p] * m.c[p, p] + m.c[q, p] * m.c[q, p]);
      for (i = 0; i < m.size; i++)
      {
            for (j = 0; j < m.size; j++)
            {
                  r.c[i, j] = 0.0;
            }
            r.c[i, i] = 1.0;
      }
      if (radius != 0.0)
      {
            if (m.c[p, p] == 0.0)
            {
                  r.c[p, p] = 0.0;
                  r.c[p, q] = 1.0;
            }
            else
            {
                  r.c[p, p] = m.c[p, p] / radius;
                  r.c[p, q] = m.c[q, p] / radius;
            }
            r.c[q, q] = r.c[p, p];
            r.c[q, p] = -r.c[p, q];
      }
}

该函数在 r.c 中构建上三角矩阵。

我标记了所有在上三角矩阵中通常不为 0 的元素，但我们对此不感兴趣，标记为“..”，因为这些元素在这种情况下都会变为 0。

要构建这个上三角矩阵，我们基本上必须让旋转矩阵的 p 和 q 从 1 运行到 N。向量 d 也必须乘以相同的旋转矩阵，原则上 r 也是如此。但由于我们并不真正关心 r，所以我们可以省略它。最终（理论上）我们得到：

C’Qx – dQ = rQ

或者

C’’x – d’ = r’

其中 Q 是所有旋转矩阵的乘积。如果我们看一下变换后的 C，我们可以看到第 4 到第 7 行的元素全为 0。这意味着 d' 中第 4 到第 7 个元素相当固定，就像 d' 的第 4 到第 7 个元素一样。

我们寻找矩阵方程的解，其中 r' 的元素平方和最小。当 r' 的第 1 到第 3 个元素为 0 时，这种情况发生。

这意味着我们只需要求解这个剩余的矩阵方程：

通过后向替换，如下所示：

public void SolveRMatrix(int maxOrder, TMatrix mi)
{
      int k, l;
      for (k = maxOrder - 1; k >= 0; k--)
      {
            for (l = maxOrder - 1; l >= k; l--)
            {
                  mi.d[k] = mi.d[k] - mi.x[l] * mi.c[k, l];
            }
            if (mi.c[k, k] != 0)
                  mi.x[k] = mi.d[k] / mi.c[k, k];
            else
                  mi.x[k] = 0;
      }
}

然后我们就完成了 :-)

速度提升

我在这里做了一个小的速度提升：在上面的示例中，在变换过程中，我们最终只关心第 1 列到第 m 列的元素，因为我们正在寻找 m 阶的多项式。如果 p 大于 m，乘以旋转矩阵只会影响大于 m 的行中的元素。因此，我们可以跳过这些乘法，让 p 只从 1 运行到 m，也就是所搜索多项式的阶数。这节省了相当多的计算时间。有了这个，将矩阵 C' 转换为上三角矩阵的函数如下所示：

public void RunGivens(int maxP, int maxQ)
{
      int i, k;
      for (i = 0; i < maxP; i++)
      {
            for (k = i + 1; k < maxQ; k++)
            {
                  if (m.c[k, i] != 0.0)
                  {
                        CalcRMatrix(i, k);
                        m.MultMatrix(r, i, k);
                  }
            }
      }
}

将 maxP 设置为阶数，将 maxQ 设置为样本数（参见我的示例项目 Smallest_Squares_Givens.zip）。

与最小二乘法相比，这是一种完全不同的方法。难道不令人着迷吗？:-)

使用代码

我的 C# 代码示例只有一个主窗口。我实现了一个用于高斯消元法求解矩阵方程的类和一个函数。在 MainWin 类的顶部定义了公共常量 samples 来固定样本数量。样本数量可以在 4 到 40 之间设置。在窗口中有 10 个用于 x 和 y 的文本框，最多可编辑 10 个样本。如果样本数量大于 10，这些文本框将变为不可见。在这种情况下，必须在应用程序启动时初始化样本值。我实现了 7 个样本点，并在 MainWin 函数中初始化它们。
所搜索多项式的阶数也在 MainWin 类的顶部定义。但 order 不是常量。它可以在应用程序窗口右下角的 1 和应用程序中的样本数量之间更改。
应用程序计算近似函数并在 pGraph 面板中的 DrawGraph 函数中绘制它。为此，它搜索 x 和 y 方向的最大值并将图形缩放到适合面板。参考点绘制为红点。多项式在 a.x[0] 到 a.x[n] 中。