算法分析中的欧拉方程：基于三数中值法的快速排序

摘要: 三数中值法的快速排序，欧拉方程的求解

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

---

在大学高等数学或考研高等数学中，我们学过几类容易求解的常微分方程的解法，特别地对于常系数线性微分方程，有一套标准化的解法。但如果是变系数的线性微分方程，也基本上很难求解了。

变系数的线性微分方程中有一类比较特殊的方程称为欧拉方程，这是有一套标准解法可以解决的。欧拉方程是指具有以下形式的变系数线性微分方程：

$$x^{n}y^{(n)} + P_{1}x^{n-1}y^{(n-1)} + \cdots + P_{n-1}xy' + P_{n}y = f(x)$$

解决这类方程的过程涉及到变量代换和微分算子，当年我们为了期末考试，为了考研，刷了很多解这种方程的题，了解到这个方程在热的传导、圆膜的振动、电磁波的传播等问题中有所应用。

本文我们介绍一个在算法分析中与欧拉方程有关的例子，还是个比较常见的算法：基于三数中值法的快速排序。

我们首先讨论一下在快速排序中引入三数中值法的动机，并给出其递归式，引入生成函数，将递归式做一些处理得到生成函数满足的微分方程，正是欧拉方程。

然后我们用微分算子法解出生成函数，得到三数中值法的算法分析结果，并与普通快速排序对比，看看引入三数中值法之后的优化效果到底怎么样。

---

快速排序的递归式

快速排序算法是算法与数据结构中的入门必学算法，其算法原理与实现可以参考相应的资料。这里为了分析方便，大致介绍一下算法流程。

对于数组 a[l:r]，将其划分为 2 个独立的更小的部分，将这两个部分分别排序，从而达到将 a[l:r] 排序的目的。

划分数组的过程中把最后位置上的元素放到它排序后的正确位置上，称其为分割元，使得所有较小元素在它前面，所有较大元素在它后面。

用两个指针完成这种划分，其中一个从左向右扫描数组，直到发现比分割元大的元素，另一个从右向左扫描数组，直到发现比分割元小的元素。然后将这两个元素交换，然后继续进行以上过程，直至两个指针相遇。此时两个指针相遇的位置就是分割元的位置。

记 $C_{N}$ 为上面程序将 $N$ 个元素排序所用的平均比较次数。$C_{0} = 0$，根据以上算法的描述，写出其递归式：

$$C_{N} = N + 1 + \frac{1}{N}\sum\limits_{j=1}\limits^{N}(C_{j-1} + C_{N - j})$$

上式的含义是，对于长 $N$ 的数组，划分子数组的过程需要比较 $N + 1$ 次；当分割元为第 $j$ 个最大元素时，划分后的子数组大小为 $j - 1$ 和 $N - j$。

在快速排序中引入三数中值法

上述快速排序算法中，如果分割元始终是当前子数组的最大值，那么会陷入到最坏的情况。因此我们希望分割元尽量是当前子数组的中位数，这样长为 $N$ 的数组在划分后的两个子数组的长度都尽可能接近 $\frac{N}{2}$。这是引入三数中值法的动机。

三数中值法是从长为 N 的数组中选取 3 个数作为样本，取其中间的那个值作为分割元。此时快速排序的递归式如下：

$$C_{N} = N + 1 + \sum\limits_{k=1}\limits^{N}\frac{(N - k)(k - 1)}{\binom{N}{3}}(C_{k-1} + C_{N - k})$$

上式的含义是，其中 $\binom{N}{3}$ 为 $N$ 个数中选取 3 项的方法数。由初等计数中的乘法原理，选出 3 个数，使得第 $k$ 个最小元是分割元的选法是 $(N - k)(k - 1)$。

于是第 $k$ 个最小元素为分割元的概率，从一般快速排序的 $\frac{1}{N}$ 变成了 $\frac{(N - k)(k - 1)}{\binom{N}{3}}$，这就是引入三数中值法带来的改变。

下面我们根据上述递归式对 $C_{N}$ 进行渐近估阶。

欧拉方程的导出

用 $N(N-1)(N-2)$ 乘以递归式两边，然后消去和式中对称的部分，得：

$$N(N-1)(N-2)C_{N} = (N + 1)N(N-1)(N-2) + 12\sum\limits_{k=1}\limits^{N}(N - k)(k - 1)C_{k-1}$$

再用 $z^{N-3}$ 乘以两边，再对 $N$ 求和，记 $C_{N}$ 的生成函数为 $C(z)$，有生成函数的性质，得到生成函数满足的微分方程：

$$C'''(z) = \frac{24}{(1-z)^{5}} + 12 \frac{C'(z)}{(1 - z)^{2}}$$

这是一个高阶方程，我们基本上不用指望得到显式的解，但是细一看，用 $(1 - z)^{3}$ 乘以两边，得：

$$(1 - z)^{3}C'''(z) = \frac{24}{(1-z)^{2}} + 12 (1 - z)C'(z)$$

每个项的 $(1 - z)$ 的次数都等于该项求导的阶数，因此这是一个欧拉方程。可以通过常微分方程相关理论得到显式的解。

微分算子法求解欧拉方程

按照解欧拉方程的常规方法，定义微分算子：

$$\psi C(z) = (1 - z)\frac{\mathrm{d}}{\mathrm{d}z} C(z)$$

于是微分方程转化为：

$$\psi(\psi + 1)(\psi + 2)C(z) = 12\psi C(z) + \frac{24}{(1-z)^{2}}$$

分解因式后，得：

$$\psi(\psi + 5)(\psi - 2)C(z) = \frac{24}{(1-z)^{2}}$$

根据上式，求解该欧拉方程可以转换为以下三个关于 $C(z)$ 的一阶微分方程：

$$\psi U(z) = \frac{24}{(1 - z)^{2}}$$ $$(\psi + 5)T(z) = U(z)$$ $$(\psi - 2)C(z) = T(z)$$

一阶微分方程是容易解出显式解的，根据微分方程相关理论，求解上述方程，最终得到：

$$U(z) = \frac{12}{(1-z)^{2}} - 12$$ $$T(z) = \frac{12}{7}\frac{1}{(1-z)^{2}} - \frac{12}{5} + \frac{24}{35}(1 - z)^{5}$$ $$C(z) = \frac{12}{7}\frac{1}{(1 - z)^{2}}\ln \frac{1}{1 - z} - \frac{54}{49}\frac{1}{(1 - z)^{2}} + \frac{6}{5} - \frac{24}{245}(1 - z)^{5}$$

算法分析结果

将解得的 $C(z)$ 展开，取其 $z^{n}$ 的系数，得：

$$C_{N} = \frac{12}{7}(N + 1)(H_{N+1} - \frac{23}{14})$$

在文章 渐近估阶的威力：欧拉-麦克劳林公式的应用 中我们推导出了普通快速排序的平均比较次数：

$$C_{N} = 2(N+1)(H_{N+1} - 1)$$

其中 $H_{N+1}$ 为调和数，同样在文章 渐近估阶的威力：欧拉-麦克劳林公式的应用 中我们也推导了其渐近估阶：

$$H_{N} = \sum\limits_{k=1}\limits^{N}\frac{1}{k} = \ln N + \gamma + o(1)$$

对比发现，可以发现三数中值的快速排序算法在平均情况下的比较次数与普通快速排序相比仅相差一个常数因子，均为 $O(N\ln N)$。

对比其常数因子可得：利用三数中值快速排序算法，可以节省大约 $14 \%$ 的用于比较的开销。

总结

本文我们分析了快排中的一个常见的优化思路：三数中值法。其动机主要是为了解决分割元取字数最小值时会造成过多无效比较的问题。

写出优化后算法的递归式，利用生成函数性质转化为微分方程，然后我们发现这竟然是我们之前考试刷过的欧拉方程，可以用换元和微分算子法解出显式解得到生成函数。

然后将生成函数展开，取其 $z^{n}$ 的系数就得到了平均情况的比较次数，将调和数的渐近估阶代入，得到比较次数的渐近估阶。

上述步骤中，求解方程得显式解，以及将显式解展开取 $z^{n}$ 系数，计算量非常大，可以用符号计算工具算。

经过以上分析，我们算是知道了三数中值法的引入对比较次数带来的优化效果，具体就是可以节省 $14 \%$ 的比较次数。

当然快排算法的开销除了比较还有交换，分割元更接近中间后，比较次数虽然少了，交换的次数有可能增加，如果要更细致的分析，需要将其考虑进去。

---