• 概率论教育程序（入门文章，可选）
  
• 线性方程组介绍
  
• 有限元法
  
• 最小化二次形式和OLS问题的示例
  
• 从最小二乘到神经网络
  

因此，另一篇文章来自“指间的数学 ”。 今天，我们继续讨论最小二乘法，但是这次是从程序员的角度出发。 这是该系列的另一篇文章，但与众不同，因为它通常不需要任何数学知识。 本文被认为是对理论的介绍，因此从基本技能出发，它需要具备打开计算机并编写五行代码的能力。 当然，我不会在这篇文章上作详细介绍，并且在不久的将来，我将出版续集。 如果我有足够的时间，那么我将用这种材料写一本书。 目标受众是程序员，因此Habr是闯入的正确位置。 总的来说，我不喜欢写公式，而且我真的很喜欢从示例中学习，在我看来这很重要-不仅要看学校董事会上的花式，而且要尝试一切。

因此，让我们开始吧。 让我们想象一下，我有一个三角表面，上面有我的脸部扫描图像（在左图中）。 我需要做些什么来增强功能，将其表面变成怪诞的面具？



在这种特殊情况下，我求解了一个名为Simeon Demi Poisson的椭圆型微分方程。 程序员程序员，让我们一起玩一个游戏：猜猜决定它的C ++代码中有多少行？ 不能调用第三方库，我们只有一个裸编译器可供使用。 答案是削减。

实际上，对于求解器而言，二十行代码就足够了。 如果您将所有内容都包括在内，则包括3D模型文件的解析器在内的所有内容都保持在两百行之内-随便吐痰。

示例1：数据平滑

让我们谈谈它是如何工作的。 让我们从头开始，假设我们有一个正则数组f，例如32个元素，其初始化如下：



然后，我们将执行以下过程一千次：对于每个单元f [i]，我们将相邻单元的平均值写入其中：f [i] =（f [i-1] + f [i + 1]）/ 2。 为了更加清楚，这是完整的代码：

import matplotlib.pyplot as plt f = [0.] * 32 f[0] = f[-1] = 1. f[18] = 2. plt.plot(f, drawstyle='steps-mid') for iter in range(1000): f[0] = f[1] for i in range(1, len(f)-1): f[i] = (f[i-1]+f[i+1])/2. f[-1] = f[-2] plt.plot(f, drawstyle='steps-mid') plt.show() 

每次迭代都会使数组的数据变得平滑，经过一千次迭代，我们将在所有单元格中获得一个恒定值。 这是前一百五十次迭代的动画：



如果您不清楚为什么会出现抗锯齿，请立即停止，握笔并尝试绘制示例，否则继续阅读毫无意义。 三角化的表面与此示例基本相同。 想象一下，对于每个顶点，我们将找到其邻居，计算其质心，然后将顶点移动到该质心，如此十次。 结果将是这样的：



当然，如果您不停止十次迭代，那么一段时间后，整个表面将以与上一个示例完全相同的方式压缩到一个点，整个数组将被填充为相同的值。

示例2：增强/减弱功能

完整的代码可在github上找到 ，但是在这里，我将给出最重要的部分，仅省略3D模型的读写。 因此，对我而言，三角模型由两个数组表示：顶点和面。 verts数组只是一组三维点，它们是多边形网格的顶点。 faces数组是一组三角形（三角形的数量等于faces.size（）），对于每个三角形，顶点数组的索引都存储在该数组中。 我在计算机图形学课程中详细介绍了数据格式和使用它。 还有第三个数组，我从前两个数组（更确切地说，仅从faces数组）开始重新计数-vvadj。 这是一个数组，对于每个顶点（二维数组的第一个索引）存储其邻居的索引（第二个索引）。

 std::vector<Vec3f> verts; std::vector<std::vector<int> > faces; std::vector<std::vector<int> > vvadj; 

我要做的第一件事是为曲面的每个顶点考虑曲率向量。 让我们说明一下：对于当前顶点v，我遍历其所有邻居n1-n4； 然后我计算它们的质心b =（n1 + n2 + n3 + n4）/ 4。 好了，可以将最终曲率矢量计算为c = vb，这与二阶导数的普通有限差分完全不同 。



直接在代码中，它看起来像这样：

  std::vector<Vec3f> curvature(verts.size(), Vec3f(0,0,0)); for (int i=0; i<(int)verts.size(); i++) { for (int j=0; j<(int)vvadj[i].size(); j++) { curvature[i] = curvature[i] - verts[vvadj[i][j]]; } curvature[i] = verts[i] + curvature[i] / (float)vvadj[i].size(); } 

好了，然后我们多次执行以下操作（请参见上图）：将顶点v移至v：= b + const * c。 请注意，如果常数等于1，那么我们的顶点将不会移动到任何地方！ 如果常数为零，则顶点将替换为相邻顶点的质心，这将使我们的表面平滑。 如果常数大于1（使用const = 2.1绘制标题图片），则顶点将沿表面曲率矢量的方向移动，从而增强了细节。 这就是它在代码中的样子：

  for (int it=0; it<100; it++) { for (int i=0; i<(int)verts.size(); i++) { Vec3f bary(0,0,0); for (int j=0; j<(int)vvadj[i].size(); j++) { bary = bary + verts[vvadj[i][j]]; } bary = bary / (float)vvadj[i].size(); verts[i] = bary + curvature[i]*2.1; // play with the coeff here } } 

顺便说一句，如果它小于1，则细节会相反地变弱（const = 0.5），但这不等于简单的平滑处理，“图像对比度”将保持不变：



请注意，我的代码生成了Wavefront .obj格式的3D模型文件，该文件是在第三方程序中渲染的。 例如，您可以在联机查看器中查看生成的模型。 如果您对渲染方法感兴趣，而不是生成模型，请阅读有关计算机图形学的课程 。

示例3：添加约束

让我们回到第一个示例，并做完全相同的事情，只是不会重写数字0、18和31下的数组元素：

 import matplotlib.pyplot as plt x = [0.] * 32 x[0] = x[31] = 1. x[18] = 2. plt.plot(x, drawstyle='steps-mid') for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1])/2. plt.plot(x, drawstyle='steps-mid') plt.show() 

我用零初始化了数组的其余“空闲”元素，但仍然用相邻元素的平均值迭代地替换它们。 这就是数组的演化在前一百五十次迭代中的样子：



很明显，这次解决方案将不会收敛到填充数组的常数元素，而是收敛到两个线性斜坡。 顺便说一句，对每个人来说真的很明显吗？ 如果没有，请尝试使用此代码，我专门给出了一个非常简短的代码示例，以便您可以彻底了解正在发生的事情。

抒情离题：线性方程组的数值解。

让我们得到通常的线性方程组：



可以将其重写，将每个方程式留在等号x_i的一侧：



让我们得到一个任意向量 近似方程组的解（例如零）。

然后，将其粘贴在系统的右侧，我们可以获取更新的解近似向量 。

更清楚地说，从x0获得x1，如下所示：



重复该过程k次，解将由向量近似

为了以防万一，编写一个递归公式：



在关于系统系数的一些假设下（例如，很明显对角元素不应为零，因为我们将它们除以了零），该过程收敛到真实解。 这种体操在文献中被称为Jacobi方法 。 当然，还有其他方法可以对线性方程组进行数值求解，例如更强大的方法，例如共轭梯度法 ，但也许Jacobi方法是最简单的方法之一。

再次是示例3，但另一方面

现在，让我们仔细看一下示例3中的主循环：

 for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1])/2. 

我从一些初始向量x开始，并对其进行了1000次更新，可疑的更新过程与Jacobi方法类似！ 让我们明确地写出这个方程组：



花一点时间，确保我的Python代码中的每次迭代严格来说都是该方程组的Jacobi方法的一个更新。 值x [0]，x [18]和x [31]对我来说分别是固定的，它们不包含在变量集中，因此它们被转移到右侧。

总的来说，我们系统中的所有方程看起来都像-x [i-1] + 2 x [i]-x [i + 1] =0。这个表达式只不过是二阶导数的普通有限差分 。 就是说，我们的方程组只是简单地规定，在任何地方（除了在点[[18]处），二阶导数都必须等于零。 记得，我说过很明显迭代应该收敛到线性斜坡吗？ 这就是为什么二阶导数的线性导数为零的原因。

您知道我们刚刚解决了Laplace方程的Dirichlet问题吗？

顺便说一句，专心的读者应该注意到，严格来讲，在我的代码中，线性方程组的求解不是通过Jacobi方法求解，而是通过Gauss-Seidel方法求解 ，这是Jacobi方法的一种优化方法：

示例4：泊松方程

让我们稍微改变一下第三个示例：每个单元不仅放置在相邻单元的质心中，而且还放置在质心中加上一些（任意）常数：

 import matplotlib.pyplot as plt x = [0.] * 32 x[0] = x[31] = 1. x[18] = 2. plt.plot(x, drawstyle='steps-mid') for iter in range(1000): for i in range(len(x)): if i in [0,18,31]: continue x[i] = (x[i-1]+x[i+1] +11./32**2 )/2. plt.plot(x, drawstyle='steps-mid') plt.show() 

在前面的示例中，我们发现放置在质心中是Laplace算子的离散化（在我们的例子中是二阶导数）。 也就是说，现在我们正在解决一个方程式系统，该系统说我们的信号必须具有一定的常数二阶导数。 二阶导数是表面的曲率； 因此，分段二次函数应该是我们系统的解决方案。 让我们检查32个样本中的样本：



数组长度为32个元素，我们的系统经过数百次迭代收敛到一个解决方案。 但是，如果我们尝试128个元素的数组呢？ 这里的一切都非常可悲，必须已经计算了数千次迭代：



Gauss-Seidel方法非常易于编程，但不适用于大型方程组。 您可以尝试使用例如多重网格方法来加快速度。 换句话说，这听起来很麻烦，但是这个想法非常原始：如果我们想要一个具有一千个元素的分辨率的解决方案，我们可以首先以十个元素进行求解，得到一个近似的近似值，然后将分辨率加倍，再次求解，依此类推，直到达到预期的结果。 实际上，它看起来像这样：



而且您不能炫耀，而使用方程组的真实求解器。 因此，我通过构造矩阵A和列b来解决同一示例，然后求解矩阵方程Ax = b：

 import numpy as np import matplotlib.pyplot as plt n=1000 x = [0.] * n x[0] = x[-1] = 1. m = n*57//100 x[m] = 2. A = np.matrix(np.zeros((n, n))) for i in range(1,n-2): A[i, i-1] = -1. A[i, i] = 2. A[i, i+1] = -1. A = A[1:-2,1:-2] A[m-2,m-1] = 0 A[m-1,m-2] = 0 b = np.matrix(np.zeros((n-3, 1))) b[0,0] = x[0] b[m-2,0] = x[m] b[m-1,0] = x[m] b[-1,0] = x[-1] for i in range(n-3): b[i,0] += 11./n**2 x2 = ((np.linalg.inv(A)*b).transpose().tolist()[0]) x2.insert(0, x[0]) x2.insert(m, x[m]) x2.append(x[-1]) plt.plot(x2, drawstyle='steps-mid') plt.show() 

这是该程序工作的结果，请注意，该解决方案立即生效：



因此，实际上，我们的函数是分段二次的（二阶导数是常数）。 在第一个示例中，我们将第二个二阶导数设置为零，第三个非零，但是到处都相同。 第二个例子是什么？ 我们通过指定曲面的曲率来求解离散泊松方程 。 让我提醒您发生了什么：我们计算了传入曲面的曲率。 如果我们通过将输出处的曲面曲率设置为等于输入处的曲面曲率（常量= 1）来解决泊松问题，则什么都不会改变。 当我们简单地增加曲率（常量= 2.1）时，就会增强面部特征。 如果const <1，则结果曲面的曲率会减小。

更新：另一个玩具作为作业

开发SquareRootOfZero提出的想法，请尝试以下代码：



这是默认结果，红色列宁是初始数据，蓝色曲线是它们的演变，在无限远处，结果收敛到一个点：



这是系数为2的结果：



作业：为什么在第二种情况下，列宁首先变成捷尔任斯基，然后又收敛到列宁，但规模更大？

结论

可以将许多数据处理任务（尤其是几何形状）表述为线性方程组的解决方案。 在本文中，我没有讲述如何构建这些系统，我的目标只是表明这是可能的。 下一篇文章的主题将不再是“为什么”，而是“如何”以及以后要使用的求解器。

顺便说一句，毕竟在文章标题中，最小的正方形。 您在文字中看到了吗？ 如果不是，那么这绝对不是吓人的，这正是对“如何？”问题的答案。 保持在线，在下一篇文章中，我将确切显示它们的隐藏位置，以及如何修改它们以访问功能非常强大的数据处理工具。 例如，在十行代码中，您可以获得：



敬请期待更多！