Math – Twisted Meadows

工作需要，最近在学高性能计算相关的知识。
这个领域，不涉足的时候觉得高深莫测，真正接触起来，不仅不神秘，而且十分有趣。
本文基于我对 LAFF-On Programming for High Performance （UTAustinX UT.PHP.16.01x）这门课的学习，结合个人在机器学习方面的经验总结形成，与大家分享。不当之处请指正。

如何实现高性能计算

从字面上来说，算得快就是高性能计算。
以此为定义，拥有超级计算机是做高性能计算的第一步，因为它从硬件性能上支持了「high performance」。
但显然，「硬件快」只是「high performance」的一个组成部分。在有限的硬件性能下，还有方法可以算得快。
例如分布式和并行化：把任务拆分成可以并行执行的多个部分，用多台机器同步运算，是个好主意。
或者给机器上插很多块显卡（GPU），利用显卡的并行计算能力，来做异构计算（由CPU执行串行部分，GPU执行并行部分）。
——如果只有一台机器，没有显卡，机器的性能也很一般呢？

本文接下来要讨论的，就是在常规架构的普通计算机单机上，如何做高性能计算编程。
特殊地，我们讨论最常用的一种运算：矩阵乘法。我们将会看到，它在时下流行的神经网络中如何被广泛应用，又是如何被最高效地实现为C代码。

机器学习中的回归问题，类似于以往各种工程上的优化问题。定义一个 Cost Function：

(1) $\begin{equation*} J(\theta) = \frac{1}{2m}\sum_{i=1}^m(h_\theta(x^i)-y^i)^2 \end{equation*}$

对这个函数求解最优化问题。

由于 $J$ 是一个关于 $\theta$ 的函数，所以目标是求一组最优的 $\theta$ ，来使 $J$ 取得最小值。

一个符合直觉的做法是Gradient Descent。也就是迭代地用 $J$ 对 $\theta$ 求导，顺着梯度最大的方向移动，直到收敛于最低点。

另一个方法则是Normal Equation，对我来说这是一个非常强的公式，可以把Gradient Descent花了那么多次迭代才解决的事用一次运算就搞定：

(2) $\begin{equation*} best \theta = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top \overrightarrow{y} \end{equation*}$

这里面的 $\mathbf{X}$ 是一个矩阵， $\overrightarrow{y}$ 是列向量， $\theta$ 其实也是列向量。
$\mathbf{X}$ 的每一行是一个样本，第 $i$ 行其实就是 $x_i^1, x_i^2, x_i^3, \dots$ 这样的一组特征(feature)，其中的每个值都可以视作特征空间里的一个维度。
对于第 $i$ 个样本，它对应的（标准）结果就是 $y_i$ ，总共 $m$ 个样本结果组成了向量 $\overrightarrow{y}$ 。
（回顾一下Cost Function：损失函数其实是定义了我们的Hypothesis $h_\theta(\mathbf{X})$ 与真实结果 $\overrightarrow{y}$ 之间的距离。以此衡量训练结果。）
再补充Hypothesis：

(3) $\begin{equation*} h_\theta(x) = \theta^\top \overrightarrow{x} = \theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n \end{equation*}$

Normal Equation 的思路类似于中学数学里求极值：对某个方程求导，然后让式子等于 0，等于 0 的点即为极值点。
但为什么上述思路能被(2)这样一个式子实现呢。

这个问题是我去年初遇到的。当时就想把阳哥哥给我的解答记录下来。拖了这么久，终于给Wordpress装了LaTeX插件，所以来记录一下：

(2)的推导方法有两种，第一种是 cost function $J$ (1)对 $\theta$ 这个向量的每一个分量求偏导

$\begin{cases} \frac{\partial J}{\partial \theta_1} = \frac{1}{m} (h_\theta(x^1)-y^1) x_1\\ \frac{\partial J}{\partial \theta_2} = \frac{1}{m} (h_\theta(x^2)-y^2) x_2\\ \dots \\ \frac{\partial J}{\partial \theta_n} = \frac{1}{m} (h_\theta(x^n)-y^n) x_n \end{cases}$

假设 $\theta$ 是个n维向量，就得到了n个式子，令这n个式子都等于0，我们就得到了一个线性方程组:

$\begin{cases} \frac{1}{m} (h_\theta(x^1)-y^1) x_1 = 0\\ \frac{1}{m} (h_\theta(x^2)-y^2) x_2 = 0\\ \dots \\ \frac{1}{m} (h_\theta(x^n)-y^n) x_n = 0 \end{cases}$

这个线性方程组可以写成矩阵相乘的形式，即

(4) $\begin{equation*} \mathbf{X}^\top\mathbf{X}\theta = \mathbf{X}^\top \overrightarrow{y} \end{equation*}$

(4)这个式子就是所谓的Normal Equation。它跟式(2)是等价的。

第二种方法要稍微高级一点儿，写出来会简洁得多。就是直接把cost function写成矩阵的形式，即

$J = (\overrightarrow{y} - \mathbf{X}\theta)^\top (\overrightarrow{y} - \mathbf{X}\theta)$

然后直接以这个式子对向量 $\theta$ 求导，也会直接得出normal equation。
至于如何直接对向量求导，你可以参考wikipedia的matrix calculus。其实本质上对向量求导就是一种记号，和对向量的每个分量求导没的本质区别。

再回头看一次这个好用的Normal Equation：

$best \theta = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top \overrightarrow{y}$

对于有 $n$ 个特征维度的 $\mathbf{X}$ 来说，用 Normal Equation 求 $\theta$ 的复杂度是 $O(n^3)$ ，而 Gradient Descent 的复杂度为 $O(kn^2)$
所以，当特征维度不多时，都可以用 Normal Equation 快速求解。当 features $n > 10000$ 时，才开始考虑使用 Gradient Descent 。

但是，这个公式显然也不是能够无条件使用的。你可以看到其中有 $(\mathbf{X}^\top\mathbf{X})^{-1}$ ，意味着括号内的矩阵必须是可逆的。
真的是这样吗？——老师说，此处的矩阵并不必须是可逆的。在代码中，对这个矩阵的求逆直接用pinv()函数来做(也就是求伪逆、广义逆)。阳哥哥说，使用广义逆来算 $\theta$ 是 $(\mathbf{X}^\top\mathbf{X})$ 不可逆时的一种传统的处理方式，可以用的原因是广义逆得出的 $\theta$ 符合一些良好的统计性质。（具体参考：高惠璇的《统计计算》）
不过， $(\mathbf{X}^\top\mathbf{X})$ 不可逆往往意味着数据模型有问题（Features数量比样本数量还大，或是非常强的collinearity），遇到不可逆的情况还是返回去检查模型比较好。

The path twists, and the future is uncertain.

分类： Math

高性能计算初步：矩阵乘法

如何实现高性能计算

对 ML 中 Normal Equation 的理解