原文：C++ Differentiable Programming: Searching For An Optimal Dither Pattern
作者：Alan Wolfe
翻译：无阻我飞扬

摘要：本文作者介绍了梯度下降算法，通过可微编程实现寻找一种最佳的图像抖动模式，详细介绍了其中的五个步骤，并通过结果展示了图像效果。读懂本文，需要有一定的高等数学知识，以下是译文。

实现这篇文章的C++代码在Github上的网址是：https：//github.com/Atrix256/DitherFindGradientDescent

神经网络现在是一个热门话题。围绕它们有许多奥秘和神秘感，但是它们的核心其实只是简单的程序，程序参数通过使用梯度下降进行调整。

（如果对神经网络感兴趣，可能会觉得这篇文章很有趣：如何用反向传播训练神经网络）

虽然梯度下降可以用于很多其它情况，事实上，甚至可以推广神经网络的核心功能来处理其它类型的程序，这正是这篇文章中所做的。

为了能够使用梯度下降来优化程序的参数，程序必须有如下构成：

它具有指定如何处理一些其它数据的参数
可以通过某种方式对它做得如何给出一个分数评价

除了上述两点之外，就像着色器程序或SIMD程序一样，希望程序尽可能无分支（if 语句）。这样做的原因是因为在理想的情况下，整个程序应该由可微操作组成。分支（if语句）导致不连续和不可微。有些方法可以处理分支，有些分支实际上不会影响结果，但这是一个应该牢记在心的好的指导方针。正因为如此，还是希望远离不可微的函数——比如可以试图使用“step”函数而不是if语句。

这篇文章将详细介绍如何在C ++中使用可微编程来实现特定的目标。结果的展示，以及生成结果的简单/无外部依赖的C ++代码位于
https://github.com/Atrix256/DitherFindGradientDescent。

首先，简单介绍一下梯度下降。

一维梯度下降

以一个 f(x)\ f(x) 形式的函数为例，它只需要一个输入，所以是一维的。

可以把这样的函数想象成在数轴上的每个点都有一个值。

可以将这些值看作一个高度，它给了一个函数形式 y=f(x)\ y=f(x) ，仍然称之为一维形式的函数，尽管它现在有两个维度。

下面看一个函数 y=3x+1\ y=3x+1

还记得直线方程式是 y=mx+b\ y=mx+b ,其中m是线（riserun \frac{rise}{run} 或 yx \frac{y}{x} ）的斜率，而b是线与Y轴交叉的位置。

在微积分中，会发现斜率m也是函数的导数：dydx \frac{dy}{dx}

斜率/导数表示x每添加1，y就会对应的增加多少。

假设在这个图上x=1 x=1 ( 这会让y=4 y=4 ）,并且假设希望从当前的位置向下走。可以通过查看这个点的斜率/导数，即3（对于线上的每个点是3）来做到这一点。由于导数是正数，这意味着向右移动会使y值变大（线向上走），向左移动会使y值变小（线向下走）。

所以，如果想向下到更小的y值，需要从x中减去值。

一种更简单的方法是，从x值中减去导数，使的y值更小。

这是一个核心的事实，它会帮助引导你克服困难：减去导数（稍后减去梯度），使值变小。减去的值通常乘以一些标量值，使其移动得更快或更慢。

如果有一个更复杂的函数会发生什么，比如y=(x−2)2 y=(x-2)^2

假设在这个图上的点x=1 x=1 ，那么y=1 y=1 。现在，沿着哪条路向下走？这个函数的导数是y=2x−4y=2x-4 ，可以将x值插入到该点的斜率/导数中：-3。

减去导数向下移，这意味着需要从x中减去一个负值; 又叫做需要添加一个值到X。

正如大家所看到的，给x增加一个值并使它向右移动，实际上就是向下移。

规则生效了，万岁！

二维梯度下降

当有不止一个维度的时候，事情会变得更复杂一点，但是并不是那么复杂，所以坚持下去！

看看这个函数z=xy z=xy

假设（x，y）在点（1,1）—在右上角：z=1z=1 ，并且假设想要向下移。现在有两个变量（x和y），而不是只有一个变量取（x）的导数。如何处理这个问题？

答案是PARTIAL导数（偏导数）。

首先，假设y是一个常数值，而不是一个变量。这会给出x：的偏导数∂ z∂ x \frac{\partial\ z}{\partial\ x} 。由此可知，如果向x加一个值，z会加多少。这是一个特殊的沿着x轴的斜率。

在这种情况下，z相对于x的偏导数就是：y。

对另一个变量做同样的事情，z相对于y的偏导数就是：x。

现在对每个变量都有偏导数，把它们放到一个矢量中。这个矢量被称为梯度，有一些吓人的符号，看起来像这样：

∇z=∇f(x,y)=(∂ z∂ x,∂ z∂ y) { \nabla z} ={ \nabla f(x,y)}=(\frac{\partial\ z}{\partial\ x} ,\frac{\partial\ z}{\partial\ y} )

对于这个函数，梯度是：

∇z=∇f(x,y)=(y,x) { \nabla z} ={ \nabla f(x,y)}=(y,x)

这使得梯度在一个具体点：

∇z=∇f(1,1)=(1,1) { \nabla z} ={ \nabla f(1,1)}=(1,1)

在最后一节中，看到导数/斜率指向函数变大的地方。梯度也是如此，它们也指向函数变大的方向！

所以，如果想向下，需要从x和y中减去值。事实上，从当前的点向下最快的方式是从x和y中减去相同的值。这是因为梯度不只是指向它变大的地方，而是指向最快的地方。所以，逆向梯度也指向最快变小的地方。

很酷吧？

可以通过查看函数图来直观的确认这一点。

在继续之前，关于斜率，导数和梯度的最后一件事情，当它们指向最大增长方向时，仅对于非线性函数的图上的无限小点有效。当沿着梯度的相反方向移动时，这将非常重要，但是要确保用非常小的步骤来帮助找到图表上的最低点。

为什么要使用梯度下降？

为什么要使用梯度下降？假设有一个函数：

w=f(x,y,z) w=f(x,y,z)

当然，可以为x，y和z选择一些随机的起始值，然后使用梯度下降找到最小的w，但是谁在乎呢？

给这些变量一些其它名称，看看这个值是否变得更加明显一点：

DamageTakenMultiplier=CalculateDamageTakenMultiplier(Armor,Dodge,Resist) DamageTakenMultiplier = CalculateDamageTakenMultiplier(Armor, Dodge, Resist)

现在，通过只改变变量的名称，可以看到，可以使用梯度下降来找到多少Armor，Dodge和resist，使得角色受到最少的伤害。现在可以告诉你如何将统计点分配给一个角色以获得最好的结果。

请注意，如果试图找到可能的最高数字，而不是最低数字，则可以将函数乘以-1，并以相同的方式执行其它所有操作。也可以做梯度ASCENT（梯度上升），这相当于乘以-1并做梯度下降。

问题

下面是一些在进行梯度下降时可能遇到的常见问题:

局部极小值—当到达碗底（这个图形就像一个碗状的），但它并不是碗的最深处。
平坦导数 —由于这些导数非常小，所以很难逃离局部地区，这也使得每一次移动都非常小。
不连续性 —问题空间（图形）在没有警告的情况下突然改变，使得梯度下降做错误的事情。

下面是一个局部最小值与全局最小值的例子。可以看到，根据在这张图表上开始的位置，如果唯一的规则是“向下”，可能会进入更深的碗，或更浅的碗。

（图片来自维基百科由KSmrq提供 http://commons.wikimedia.org/wiki/File:Extrema_example.svg, GFDL 1.2,
https://commons.wikimedia.org/w/index.php?curid=6870865）

下面是一个平坦导数的例子。可以想象，如果在x = 1时，可以看到导数会让左边减小y值，但是这个数字非常非常小。这是一个问题，因为通常在减去乘数之前乘以导数或梯度，所以只需要朝目标迈出一小步。

也可以想成一个完全平坦导数，而这个导数恰好为0。在这种情况下，无论数字有多大或多小，都可以乘以导数，根本不会移动！

下面是一个不连续的函数，如果x小于0.5，则值为1，否则为x。这基本上显示了在可微编程中使用if语句时会发生什么。如果从右侧开始，应该向左移动来提高分数。然而，它会一直持续向左移，直到x小于0.5，那么分数会突然变得更糟，导数将变成0，就被卡在这个地方了！

解决这些问题的方法有很多，但都是深层次的话题。如果没有别的，应该知道这些问题是存在的，所以可以知道它们什么时候会影响，和/或为什么应该避免它们，如果有所选择。

如果想避免计算怎么办？

假设没有计算出所有这些偏导数。或者更务实一点，不想坐下来手工计算一些通用的C ++代码的梯度函数！

有一些好消息。

虽然需要梯度的偏导数，不需要做这些计算来获取它们。

以下是其它一些获得偏导数的方法：

有限差分—概念上超级简单，但计算速度慢，并不总是非常精确。更多信息：有限差分
反向传播—神经网络使用什么。也称为自动微分反向模式。快但有点复杂。更多的信息：如何训练反向传播的神经网络
-对偶数—也称为正向模式自动微分。在相同的邻域中，速度不如反向模式快。对于程序员来说超级，超级方便，超赞。更多信息：对偶数和自动微分

仔细猜猜要使用哪一个？是的，对偶数！

一言以蔽之，如果有代码使用了浮点数，则可以将其更改为使用模板类型。然后，在代码中使用对偶数，而不是浮点数。得到的输出将是代码输出的特定值，也是代码在这个值的梯度。更好的是，这不是一个数值方法（它不是一个近似值），它是分析的方法（确切的说）。

对偶数是惊人的！

由于已将代码模板化，因此当不想要或不需要梯度时，仍然可以将它用于浮点数模式。

可微编程/梯度下降梗概

这里是将要使用梯度下降进行可微编程的一般框架：

初始化参数为随机（但有效）的值，将其存储在对偶数中。
运行代码，以对偶数作为输入，说明它是如何工作的。
把结果（这是对偶数）放入一个得分函数给出一个分数。通常情况下，分数越小越好。如果不是，只需乘以-1即可。
因为做了这项工作，并使用对偶数来计算分数，现在有一个梯度来描述如何调整参数使分数更好。
使用梯度调整参数并返回到第2步。重复，直到想要的任何退出条件被选中：也许当一定次数的迭代发生，或者当分数低于某个值时。

这是策略，深入讨论将要探讨的具体问题。

寻找一个理想的图像抖动模式

以下是要解决的问题：

想要找到一个3×3抖动模式，当用它抖动一个图像（通过在图像上反复重复3×3模式），然后以特定的数量模糊结果，它是如此接近被同样数量的模糊的原始图像。

这听起来有点挑战性对不对？这不是真的那么糟糕，不要担心（：

代码的步骤（可微的）是：

抖动源图像
模糊结果
模糊源图像
计算它们有多相似的分数
使用所有这些梯度下降优化抖动模式

再一次，需要用对偶数来微积分这个事情，这样就得到一个梯度来修改抖动模式，以便有更好的得分。

步骤1—抖动源图像

抖动图像是一个非常简单的过程。

要抖动它，以便将灰度图像作为输入，并使用抖动模式将其转换为黑白图像。

（如果开始使用彩色图像，则会显示如何将其转换为灰度：将RGB转换为灰度）

对于源图像中的每个像素（x，y），查看抖动模式中的像素（x％3，y％3），如果抖动模式像素小于源，则会写出黑色像素，否则将写一个白色的像素。

if (sourcePixel(x,y) < ditherPixel(x%3, y%3))
pixelOut(x,y) = 0.0;
else
pixelOut(x,y) = 1.0;

但是有个问题，这是一个分支，它会造成不连续性，这将不能有很好的导数来帮助实现目标。

写上面的抖动操作的另一种方法是这样写的：

difference = ditherPixel(x%3, y%3) - sourcePixel(x,y);
pixelOut(x,y) = step(difference);

这里的“step”是“阶跃函数”（是一种特殊的连续时间函数，是一个从0跳变到1的过程，属于奇异函数），如果x> = 0则为1，否则为0。

（来自维基百科的图像由Omegatron提供（自己的工作）[CC BY-SA 3.0
（https://creativecommons.org/licenses/by-sa/3.0）或GFDL（http://www.gnu.org/copyleft/fdl.HTML）％5D，通过维基共享资源）

去掉了分支（if语句），但仍然有一个不连续的函数。

幸运的是，可以用一个近似阶跃函数的其它函数。比如使用这样的公式0.5+atan(100∗x)/pi0.5+atan(100*x) /pi:

不幸的是，结果不是那么好，所以把它切换到 0.5+atan(1000∗x)/pi0.5+atan(1000*x) /pi，得到了一个更好的结果：

这个函数确实有平坦导数问题，但它的效果却很好。在这种情况下，平坦导数似乎并不是一个大问题。

为了把它们集中在一起，所使用的抖动像素的可微版本看起来像这样：

difference = ditherPixel(x%3, y%3) - sourcePixel(x,y);
pixelOut(x,y) = 0.5+atan(10000.0f * difference) / pi;

作为这个抖动过程的输入，采取：

源图像
一个3×3抖动模式，其中每个像素是一个对偶数

作为输出，这个抖动过程如下：

抖动的图像被转换为黑白（每个像素的值为1.0或0.0）
它与源图像大小相同
每个像素是具有9个导数的对偶数。每个抖动像素有一个导数

步骤2—模糊结果

模糊抖动的结果并不是那么困难。本文中使用了高斯模糊（也叫高斯平滑，是在Adobe Photoshop、GIMP以及Paint.NET等图像处理软件中广泛使用的处理效果，通常用它来减少图像噪声以及降低细节层次），但用其它模糊也很容易。

这里有一些高斯模糊代码（可以从这篇博文：高斯模糊阅读），在适当的地方把它转换为使用模板类型，而不是浮点/像素，也确保没有分支或任何不连续。

幸好这里并没有什么可修正的，所以不是太困难。

这使得抖动的结果是对偶数每像素的，并对它们进行高斯模糊，保留并正确地修改梯度（导数），就如同做高斯模糊。

步骤3 —模糊源图像

模糊源图像很容易，因为最后一步做了一个普通的高斯模糊函数。使用通用的高斯模糊函数来模糊图像。这不需要作为对偶数来完成，所以它是正常像素和常规像素。

大家可能想知道为什么这个部分不需要作为对偶数来完成。

简单点回答是因为这些值不依赖于抖动模式（这是用导数跟踪的）。
更数学的解释是，实际上可以考虑这些对偶数，它们只是有一个零的梯度，因为它们本质上是常数，而且与函数的参数无关。梯度只会隐含地为零，就像可能引入该函数的任何其它常数值一样。

步骤4 —计算相似度

接下来需要计算抖动和模糊结果（由对偶数组成）和源模糊（由正常像素组成）之间的相似度分数。

使用的相似度分数只是MSE或“均方误差”。

为了计算MSE，对每一个像素做如下操作：

error = ditheredBlurredImage(x,y) - blurredImage(x,y);
errorSquared = error * error;

在对每个像素做平方误差之后，只需取它们的平均值就可以得到MSE。

关于MSE的一个有趣的事情是，因为错误是取平方的，所以相比较大的错误，MSE更偏向于较小的错误，这是一个很好的属性。

关于MSE的不太好的属性是，它可能决定某个东西在数学上是一个很小的差别，即使人类会认为这是一个巨大的差异感知。反之亦然。尽管如此，选择用它，因为它很简单，并且最终得到了相当好的结果。

如果想深入了解“感知相似度分数的图像”，看看这些链接：

SSIM - http://www.cns.nyu.edu/~lcv/ssim/
GMSD - https://arxiv.org/ftp/arxiv/papers/1308/1308.3052.pdf
多尺度SSIM -http：//www.cns.nyu.edu/~zwang/files/papers/msssim.pdf

在这个步骤之后，得到一个MSE值，该值表示图像的相似程度。较低的值意味着较低的平均平方误差，所以较低的数字确实更好。

另一个优点是，MSE的值是一个带有9个偏导数梯度的对偶数，用于描述在调整每个参数时MSE的变化量。

该梯度告诉我们如何调整参数（3×3抖动像素！）以降低MSE的值！

步骤5— 把它们放在一起

现在是把所有这些放在一起的时候了，并且使用梯度下降来使抖动模式更好。

下面是程序的运行方式：

将3×3抖动模式初始化为随机值，将梯度中的导数设置为1.0，代表它们所代表的变量。
做这个循环的1000次迭代：
1. 抖动并模糊源图像
2. 与源图像模糊相比，计算该结果的MSE
3. 使用MSE值的梯度，从抖动模式中的每个像素中减去相应的偏导数，但通过“学习率”缩放偏导数
输出最好的结果

在从0处开始循环迭代，学习率从3.0开始，但是在每次迭代时都会衰减，在迭代到999时衰减到0.1。它从1开始以帮助逃离局部最小值，并且在最后使用非常小的速率来尝试并深入到最小值的发现。

在调整抖动模式像素之后，将它们锁定在0和1之间。

还需要提到的一点是，当正在做梯度下降的时候，会跟踪看到的最佳得分的抖动模式。

这样，在经过1000次迭代之后，如果看到了比现在更好的，就用它代替最终的结果。

如果正确地调整参数（学习率，迭代等等），想来，这操作应该不会经常出现，但总的来说，最终状态并不是遇到的最佳状态，所以在大多数情况下，这是一个获得更好结果的好方法。

结果

有没有注意到这个帖子叫做“寻找一种理想的抖动模式”而不是“发现一种理想的抖动模式”？（：

结果很不错，但是还有可能会更好。尽管如此，这里所讨论的技术是一个很好的开端，沿着可微程序设计的道路，以及相类似的主题。

下面是一些能够用代码得到的结果。点击查看完整大小的图像，缩小的图像有锯齿问题。

图像从左到右分别是：原始的图像，使用的抖动模式（重复）的图像，抖动的图像，模糊的抖动图像，最后是模糊的原始图像。目的是使最后两幅图像看起来尽可能接近，使用MSE作为它们多么接近的度量标准。

下图是使用sigma为10的高斯模糊的起始状态：

下图是梯度下降的1000次迭代之后。注意到顶部的黑色块与其开始的地方相比已经消失了。

下图是使用高斯模糊sigma 为1的起始状态：

下图是1000次迭代之后，这是相当不错的结果：

最后，下图没有任何模糊：

经过1000次迭代之后，它实际上看起来更糟！

完全不用模糊会产生一些可怕的结果。模糊赋予了算法在如何成功上更多的自由，而没有模糊，找到解决方案的余地则少得多。

在MSE计算之前使用模糊的另一个好处是，模糊是低通滤波器。这意味着在MSE计算之前，更高的频率消失了。这样做的结果是该算法将有利于更接近蓝色噪声抖动的结果。相当整齐对吧？！

结束语

希望大家通过可微编程和梯度下降享受这个旅程，希望大家能够跟上。

以下是超出这篇文章讨论的一些可能有趣的事情：

让它从一组图像中学习，而不仅仅是这个单一的图像。这应该有助于防止“过度拟合”，并让它找到一个适用于所有图像的抖动模式，而不仅仅是这一个特定的图像。
使用单独的一组图像来衡量结果的准确性，该结果不作为训练的一部分，以帮助证明它确实没有过度拟合训练数据。
尝试在学习中应用“小腐败”，以防止过度拟合或陷入局部最小 —其中一个想法是每个导数有一定的几率，不要将更改应用于抖动模式像素。这会增加梯度下降的随机性，而不是只是沿着最陡峭的方向前进。
可以生成抖动模式的公式，优化该公式的系数/项，而不是优化抖动模式。如果得到了很好的结果，最终会得到一个可以用于抖动的公式，而不是一个模式，对于避免在像素着色器中读取纹理来进行抖动，这可能是一个很好的例子。

我不是一个数据科学家或机器学习专家，所以有很多改进要做。这里所做的和其它算法有很多重叠之处，无论是在机器学习领域还是在机器学习领域之外。

首先，可以使用牛顿的梯度下降法，它可以通过在计算中使用二阶导数来更快地找到最小值。

如果喜欢这篇文章，请查看本书，以获取有关梯度下降（模拟退火，遗传算法等）算法的更多细节。这是一本很好的书，很容易阅读！基础算法

1月13日，SDCC 2017之数据库线上峰会即将强势来袭，秉承干货实料（案例）的内容原则，邀请了来自阿里巴巴、腾讯、微博、网易等多家企业的数据库专家及高校研究学者，围绕Oracle、MySQL、PostgreSQL、Redis等热点数据库技术展开，从核心技术的深挖到高可用实践的剖析，打造精华压缩式分享，举一反三，思辨互搏，报名及更多详情可点击此处查看。

C++可微编程：寻找一种最佳的图像抖动模式相关推荐

5种最佳像素化图像的方法
5种最佳像素化图像的方法 1. 什么是像素化? 2. 像素化有什么用? 3. 如何像素化图像? 参考 Pixelate 像素化这篇博客将讨论像素化及如何以五种最佳方式对图像进行像素化.有时希望在在线 ...
寻找最佳特征维度_寻找5种最佳设计模式书
寻找最佳特征维度 Sometimes back I wrote an article for Best Core Java Books for beginners. Today I am sharin ...
10个最佳价格行动交易模式
这10种最好的价格动作交易模式是我最喜欢的,成功的交易者每天都用这些模式赚钱.要灵活,因为每个都有很多变化.如果你对各种可能性保持开放的心态,你将开始每天.在每个市场.在每个时间段看到它们. 所以,哪 ...
编程最佳字体_网站上使用的10种最佳装饰字体
编程最佳字体我们选择的字体在网站中起着主要作用. 正确的字体组合可以使您的网站脱颖而出. 此外,他们与您的品牌紧密结合,以传递正确的信息. 在各种字体样式中,装饰字体是一个不错的选择,如果您想产生影 ...
键盘输入一个高精度的正整数ｎ，去掉任意ｓ个数字后剩下的数字按原左右次序将组成一个新的正整数。编程对给定的ｎ和ｓ，寻找一种方案，使得剩下的数最小
键盘输入一个高精度的正整数n,去掉任意s个数字后剩下的数字按原左右次序将组成一个新的正整数.编程对给定的n和s,寻找一种方案,使得剩下的数最小. 如输入: 1785434 输出: 13 #includ ...
键盘输入一个高精度的正整数N，去掉其中任意S个数字后剩下的数字按原左右次序将组成一个新的正整数。编程对给定的N和S，寻找一种方案使得剩下的数字组成的新数最小。（C++）（贪心法）
算法经典题型12 键盘输入一个高精度的正整数N,去掉其中任意S个数字后剩下的数字按原左右次序将组成一个新的正整数.编程对给定的N和S,寻找一种方案使得剩下的数字组成的新数最小.输入数据均不需判错.输出 ...
python开源流程图软件_适用于Linux的10种最佳流程图和图表软件
图表是我们联系信息并处理其重要性的好方法. 它们有助于沟通关系和抽象信息,并使我们可视化概念. 流程图和图表工具可用于从基本工作流程图到复杂网络图 ,组织图, BPMN ( 业务流程模型和表示法 ), ...
软件工程师转正申请_这是申请软件工程师工作的4种最佳方法-以及如何使用它们。...
软件工程师转正申请 by YK Sugi 由YK Sugi 这是适用于软件工程师工作的最佳方法,以及确切的使用方法. (Here are the best methods for applying t ...
19种最佳HTML5和JavaScript游戏引擎和模板
您想编写在线游戏吗? 您很幸运,对于任何在线游戏开发人员和有志成为在线游戏开发人员的人来说,这都是一个美好的时光. Flash的时代已经过去,这为数十种(也许数百种!)不同的引擎,库和常见Web语言的 ...

C++可微编程：寻找一种最佳的图像抖动模式