数值解法入门

<pre />

• 下载 SolverDemoCode - 615.82 KB

引言

最常见的问题之一是计算方程的答案。如果方程是一个简单方程，只有一个未知变量，并且可以分离变量，那么求解很容易，例如计算 y

y = a*x +b, y = a0 + a1*x + a2*x^2, y = a0*e^x, y = x^(a/x + b)

当根据测量到的 x、y 值求解 a、b 时，问题就变得困难得多。常见示例包括从一组测量值中求解 GPS 接收器的位置，将 NMR 回波拟合到洛伦兹分布，从记录的数据中计算热系统的表观时间常数，或计算超越方程的解。

数值求解的经典技术是计算一个误差函数，并搜索（猜测）一个能最小化误差的解。不同的问题适用于不同的方法。但它通常是一个迭代过程。

本文介绍了一些经典方法，并提出了一种新颖的方法，该方法结合了所介绍的一些技术的优点，通常很有用。通常，“最佳”技术取决于可用的计算能力和时间、所需的精度以及误差函数的平滑度和形状。寻找最佳技术在一定程度上是主观的，并且相对于问题而言，它是一门艺术。提供的这种新方法旨在提供一个有趣的例子。

提供的源代码相对容易使用，调用者定义一个误差函数（委托/回调）和函数系数的初始估计。作为 API，无需理解其背后的数学原理，即可通过几次调用解决大量问题。

背景

这些技术非常灵活。如果您能够将您想要（或不想要）的东西描述为一个要最小化的误差函数，并且该函数最好具有平滑的单调斜率，那么该解将快速收敛。

示例 1

假设您有一组数据点，需要计算“最佳”拟合该线的线。

y = mx + b

在这种情况下，y 和 x 是已知的，但给定任意两个 y 和 x 对，由于测量噪声，将计算出不同的 m 和 b 对。目标是最小化点的误差。

误差函数可以是，给定一个“b”和一个“m”，为每个已知的 x 计算一个预测的 y，然后

Error = sum( (y-ypred)*(y-ypred)) 对于所有 y 和 yPredicted。

这是一个用最小二乘法解决的简单问题，提供的源代码演示了如何将直线或二次曲线拟合到数据集。提供了线性求解器和二次求解器（Solve_FirstOrder、Solve_SecondOrder）的调用。请注意，这最小化的是误差的平方而不是 delta 误差，当存在 x、y 测量噪声时，delta 误差通常更可取。

LinearLeastSquares lls = new LinearLeastSquares();
double [] x = {1, 2, 3, 4, 5};
double [] y = {3, 5, 7, 9, 11};
lls.Solve_FirstOrder(x, y);
if ((1 = lls.Solution[0])&& (2 == lls.Solution[1]))
{
    // Found the correct line that passes through y at a given x.
    // y = lls.Solution[1] * x + lls.Solution[0];</p>
}
 

示例 2

根据多个有噪声的距离测量值或已知点（GPS / 三角测量）的角度测量值来计算您所在的位置。在这种情况下，如果您位于您猜测（系数）所在的位置，那么实际测量值将具有某个值，预测值和实际值之间的差异就是误差。在这种情况下，误差通过值相对于参数的导数进行转换，将误差转换为对估计值的更新。这就是非线性最小二乘法或高斯-牛顿方法的基础。其最简单的形式是

J*dA = R

dA 是系数的调整（作为向量），

R 是每个 x 位置的误差（作为向量）。

J 是一个称为雅可比矩阵的矩阵，它将误差坐标转换为系数。

1. A = 初始估计

2. dA = J^-1 * R

3. A = A + dA

4. 如果 R 的幅度足够小，则停止，否则转到步骤 2。

矩阵在测量数量和系数数量上都会增长，而且通常不是方阵，当进行平方运算并乘以 J 的转置时，有时它是不可逆的。在这些情况下，会使用伪逆。伪逆的问题在于，由于舍入误差，它有时不准确，不同的矩阵包使用不同的方法来最小化问题。有时这种方法不收敛到解而不进行修改。

因此，不提供迭代非线性代码示例，因为它会太大太复杂，而且通常，具有创造性的猜测方案会收敛得更快。

示例 3

假设热量从 A 传递到 B，再从 B 传递到 C。我们可能知道 A 处的加热器电流以及 C 处的温度，但无法直接了解热量如何通过 B 传递。

B = (B-A)*e^(-t/tba) - (C-B)*e^(-t/tbc)

您拥有 C 和 A 的值，但没有 B，也没有时间常数 tba 和 tbc，指数函数远非线性。如果存在一个唯一的解可以最小化误差（我们测量的值减去猜测的常数所预测的值），那么问题就可以解决。在这种情况下，非线性最小二乘法也可以工作，但随着系数的增加，该解可能会收敛到不同的局部最小值，具体取决于不同的初始估计，而不会收敛到理想的答案。

// Given many x and y, solve for a, b, c.
// These are the "real" parameters we're trying to guess at.
double a = 1.001, b = 2.002, c = 3.003;

// Generate the fake data, y = a*x^2 + b*x*e^x + c*e^(x + 1)
for (int i = 0; i < data.Length; i++) // data.Length is 100.
{
    double di = (double)(i * 0.1);
    data[i] = a*di*di + b*di*System.Math.Exp(di) + c*System.Math.Exp(di + 1);
}

// Set up each parameter's search criteria.
Solver.Parameter pa = new Solver.Parameter(
                "a",     // Name.
                1.05,       // Center
                1.2,     // Max
                0.8,     // Min
                0.1,     // Initial step size
                0.0001);  // Final step size
Solver.Parameter pb = new Solver.Parameter("b", 2.05, 2.1, 1.9, 0.1, 0.0001);
Solver.Parameter pc = new Solver.Parameter("c", 3.05, 3.1, 2.9, 0.01, 0.0001);
s.Parameters.Add(pa);
s.Parameters.Add(pb);
s.Parameters.Add(pc);

// Pre-Solve
//   ComputeError    The error function.
//   200             200 generations, each generation is only 10 calculations
//   10              10 walks per variable before moving on.
//   1.25            How much over testing is done per generation.
//   0.00005         After reaching this error stop trying.
Solver s = new Solver();
s.Solve_DirectSearch(ComputeError, 10, 10, 1.25, 0.00005);
pa.CenterValue = pa.Solution;
pb.CenterValue = pb.Solution;
pc.CenterValue = pc.Solution;

// Now re-solve with slope variate, slope variate assumes your close enough
// to the solution that the errors can more or less fit a quadratic and the
// trench of a sadle point can be "walked" to improve the solution.
double err = s.Solve_SlopeVariate(ComputeError, 20, 5, 1.25, 0.00000000005);

// Our provided "error" we can convert to stdev error.
err = System.Math.Sqrt(err) / data.Length;
label1.Text = 
    "Actual A, B, C = " + a + ", " + b + ", " + c + "\r\n” + 
    “Solved A, B, C = " + pa.Solution + ", " + pb.Solution + ", " + pc.Solution + 
    "\r\nError = " + err; 

理论

作者使用的大多数方法或文档中看到的都需要一个误差函数，然后进行猜测，改进猜测，然后重复，大致的草图如下：

1. 初始猜测
2. 计算误差
3. 改进猜测
4. 计算误差
5. 如果误差不够小，转到步骤 3，否则重复。

每种算法的变化在于如何进行改进步骤。

遗传算法：

创建一组可能的系数，然后用一个“最佳”或几个“最佳”值重复。

直接搜索：

一次更改一个参数，同时减小步长重复。

最小二乘法：

将误差视为从解到解决方案的可翻译向量空间。

斜率变化

拟合（最小二乘法）一个系数对另一个系数的影响，然后成对进行直接搜索。

关键特性比较

将误差与任意两个系数作图会产生一个表面图。想象一下试图找到地图上的最低点，但只知道去过的地方的海拔。

· 最小二乘法将误差“压平”成一个解，但如果地图太复杂，那就不行了。

· 迭代最小二乘法利用斜率来确定下一步尝试的方向和大小，但如果存在二阶导数（斜率的变化率，或 Hessian）或更高阶的东西，那么方向和/或大小可能会导致相对相同的误差出现在新位置，解可能不会收敛。 反转 大型矩阵耗时耗力。

· 使用高阶导数（Hessian 矩阵）可以在曲线非常平滑（无噪声）且表面上的凹陷（局部最小值）很少的情况下有所帮助，算法可能会锁定在其中。

· 检查几何模式以计算搜索区域（单纯形法）或网格模式，然后迭代地细化搜索区域（遗传算法）以围绕可能的候选点，当存在噪声和局部最小值时效果很好，但这需要大量的计算。对于 5 个系数，一次 20 个测试点搜索是 20^5 = 3,200,000 次误差计算。如果每次误差计算都涉及计算 10k 个数据点的误差之和，并进行 20 次迭代，这可能需要 6400 亿次计算。以 10 Hz 解决相同的问题可能需要超过 20 GHz 的处理器。

· 一次只更改一个参数称为模式搜索或直接搜索，或 Hook-Jeeves 算法。在这种情况下，从一个位置开始，沿着一条直线一次探索解决方案，并改变步长。通常，如果设置得当，Hook-Jeeves 算法会过冲/欠冲并探索以找到解决方案。缺点是它可能会锁定在局部最小值，或者可能会跳过局部最小值。对于大量的系数，按顺序更改每个系数也可能导致跳过局部最小值。测试 5 个系数和 10k 个数据点可能需要 1000 次迭代或 5000 万次计算，而不是 6400 亿次。在 10Hz 的情况下，这可能需要使用双核处理器、GPU 或 ASIC 来完成。

新算法的理论

算法可能无法收敛到正确解的情况有很多。假设两个系数的误差图有一个狭窄的深沟，与系数轴成一定角度。这通常是其他擅长初步收敛到正确解的方法停止收敛的原因。在这种情况下，大多数搜索不会锁定到沟槽并沿着它移动，而是锁定在某个横截面上，并减小步长而无法收敛到正确解。噪声和局部最小值会使情况变得更糟。

这是两个系数在复杂方程中的误差图（越浅表示误差越小，中心为理想解）。

误差有一个从左到右的占主导地位的沟槽，大多数算法都会锁定到这个沟槽，有些会锁定到中心沟槽，但大多数无法很好地在其中收敛。在中心对角线井中优化的解需要算法检测、改变并沿着该点进行跟踪。

有几种方法可以通过 Hessian 或斜率“旋转”到沟槽，但它们容易受到噪声的影响，并且在沟槽不弯曲且表面没有涟漪时效果最好。

对于两个系数 a 和 b，在给定 b 的情况下，通过图的切片，以一定步长采样几次，得到一条直线。当这条直线可以通过最小二乘法拟合曲线来近似时，并且该曲线有一个最低点，那么这个最低点就是该切片中沟槽的最低点。连接多个切片的“最低”点，然后再次进行最小二乘法拟合，可以得到代表 a 和 b 系数沟槽底部的直线。由于噪声等原因，当步长足够大时，它会跳过局部最小值和其他小的误差表面波动。该算法也会跟踪弯曲的沟槽。随着步长的减小，算法将继续收敛。

为确保改进是真实的，而不是由于噪声，改进是成对进行的，并且仅当从“a”开始或从“b”开始可以进行相同的改进时才应用。

该算法有几个很好的特性

• · 轴成对选择（所有组合），因此如果任何组合可以改进，算法将快速收敛。
• · 以一定间隔进行子采样，这可以锁定误差函数的主要特征，并倾向于忽略局部最小值。
• · 当算法锁定在局部最小值时，改进另一对轴的影响通常会“弹出”系统，使其脱离该最小值。
• · 算法收敛速度极快。
• · 使用最小二乘法拟合二次曲线，并假设另一个算法已提供相对较好的初始估计。

该算法有一个明显的弱点，但可以通过将其与其他算法结合来解决。如果要在检查的值范围内不存在占主导地位的最小值，该算法将不起作用。当使用遗传网格搜索提供可能的“估计”位置时，这不是问题。

性能

性能与所选方程有关。有大量的方程集最适合每种算法，以下演示了一个典型案例，其中“斜率变化”方法比其他方法效果更好。

y = (1.001)*x^2 + (2.002)*x*e^x + (3.003)*e^(x + 1)

Paris Genetic 0.9, 2, 3.00999 Error 0.4099, 10k error function evaluations.

Paris Genetic Multiple 0.9, 2, 3.00999 Error 0.4099, 40k error function evaluations.

Direct Search 1.13, 2.00259, 3.00055, Error of 0.2157, 2k error function evaluations

Slope Variation 0.96177, 2.00163, 3.00439, Error of 0.047, 200 error function evaluations.

Slope Variation 1.000994, 2.0019999, 3.00300, Error of 6.35 x10^-6, 19.7k error function evaluations.

Paris 和 Direct search 技术只能略微改善误差，并且收敛缓慢，甚至在后期几乎没有收敛。Slope Variation 方法收敛迅速，并且持续快速收敛，而其他技术则不能。

结论

当与一种初始方法配对以生成精炼估计时，新方法通常会收敛得更快，并且几乎可以以任意精度继续收敛。这些特性使其成为快速精炼可能解的良好方法。

参考文献

zip 文件包含一个 pdf 文件，其中包含许多方法的详细数学推导。大多数使用与维基百科中所示相同的系数名称或相似的名称，以避免混淆。

http://en.wikipedia.org/wiki/Non-linear_least_squares

http://en.wikipedia.org/wiki/Genetic_algorithm

http://en.wikipedia.org/wiki/Pattern_search_%28optimization%29