概述
首先，您有必要了解一下，为什么我们需要用分布式算法来计算我们的组合优化问题（一个QAP问题）。
在实际生产中，即使我们基于传统的MV理论构建标准二次规划模型，使得二次项系数是正定矩阵，然后约束条件采用线性约束，最后在R中调用quadprog包来求解这个一定有解的二次规划问题。然而，在实际应用中，quadprog内置算法为内点法，从内部向边界迭代穷举计算最优解，这样的方法往往在高维统计或约束条件较多时失效或者寻优时间过长。
针对这个痛点，我们研究了基于ADMM分布式算法的Spark实现来满足实际工程中必然有解的二次规划问题。
什么是对偶理论
线性规划是优化算法中最基础的形式，在线性规划早期发展中最重要的发现是对偶问题，即每一个线性规划问题（称为原始问题）有一个与它对应的对偶线性规划问题（称为对偶问题）。许多原始问题非常复杂的问题，在转化为对偶问题后就变得容易了。比如有100个变量，2个约束的线性规划原问题转化为2个变量，100个约束的对偶问题，利用作图法就可以直观解释问题的求解过程。
对偶上升
对于凸函数的优化问题，对偶上升法核心思想就是引入一个对偶变量，然后利用交替优化的思路，一边优化原问题，一边优化对偶问题，使得两者同时达到最优。一个凸函数的对偶函数其实就是原凸函数的一个下界，因此可以证明一个较好的性质：在强对偶性假设下，即最小化原凸函数（primal）等价于最大化对偶函数（dual），两者会同时达到最优。这种转化可以将原来很多的参数约束条件变得少了很多，以利于做优化。具体表述如下：
$$begin{array}{lc}
min & f(x)\
s.t. & Ax = b \
end{array}$$
$$\Longrightarrow L(x, y) = f(x) + y^T(Ax - b)$$
$${\Longrightarrow} g(y) = \inf_x L(x, y)$$
在强对偶性的假设下，原问题和对偶问题同时达到最优。
$$x^{\star} = \arg\min_x L(x, y^{\star})$$
因此，若对偶函数$g(y)$可导，便可以利用梯度上升法，交替更新参数，使得同时收敛到最优。迭代如下：
begin{split}
x^{k + 1} : & =argmin_x L(x, y^k) quad text{($x$-最小化步)} \
y^{k + 1} : & = y^k + alpha^k nabla g(y) = y^k + alpha^k (Ax^{k + 1} - b) quad text{(对偶变量更新，$alpha^k$是步长)} \
end{split}
当$g$不可微的时候也可以将其转化下，成为一个所谓的次梯度法，虽然看起来不错，简单证明下即可知道 $x^k$ 和 $y^k$ 同时可达到最优，但是上述条件要求很苛刻：$f(x)$要求严格凸，并且要求$α$选择有比较合适。一般应用中都不会满足（比如$f(x)$是一个非零的仿射函数），因此对偶上升不会直接应用。
对偶分解
虽然对偶上升方法有缺陷，要求有些严格，但是他有一个非常好的性质，当目标函数$f$是可分的（separable）时候（参数抑或约束可分），整个问题可以拆解成多个子参数问题，分块优化后汇集起来整体更新。这样非常有利于并行化处理。形式化阐述如下：
begin{array}{lc}
min & f(x) = sum^N_{i = 1} f_i(x_i), x_i in mathbf{R}^{n_i}, x in mathbf{R}^n \
s.t. & Ax = sum^N_{i = 1} A_i x_i = b, quad text{(对$A$矩阵按列切分开)} \
end{array}
$${\Longrightarrow L(x, y) = \sum^N_{i = 1}L_i(x_i, y) = \sum^N_{i = 1}(f_i(x_i) + y^TA_ix_i - \frac{1}{N}y^Tb)}$$
因此可以看到其实下面在迭代优化时，$x$-minimization步即可以拆分为多个子问题的并行优化，对偶变量更新不变这对于约束多时还是很有用的。
begin{split}
x_i^{k + 1} : & =argmin_x L_i(x_i, y^k) quad text{(多个$x_i$并行最小化步)} \
y^{k + 1} : & = y^k + alpha^k nabla g(y) = y^k + alpha^k (Ax^{k + 1} - b) quad text{(汇集整体的$x$，然后对偶变量更新)} \
end{split}
对偶分解是非常经典的优化方法，可追溯到1960年代。但是这种想法对后面的分布式优化方法影响较大，比如近期的图结构优化问题。
增广拉格朗日乘数法
从上面可以看到对偶上升方法对于目标函数要求比较苛刻，为了放松假设条件，同时比较好优化，于是就有了增广拉格朗日方法，目的就是放松对于$f(x)$严格凸的假设和其他一些条件，同时还能使得算法更加稳健。
$$L_{rho}(x, y) = f(x) + y^T(Ax - b) + frac{rho}{2}|Ax - b|^2_2
Longrightarrow
begin{array}{lc}
min & f(x) + frac{rho}{2}|Ax - b|^2_2 \
s.t. & Ax = b \
end{array}$$
从上面可以看到该问题等价于最初的问题，因为只要是可行解对目标函数就没有影响。但是加了后面的$(rho/2)|Ax - b|^2_2$惩罚项的好处是使得对偶函数$g_{rho}(y) = inf_x L_{rho}(x, y)$在更一般的条件下可导。计算过程与之前的对偶上升基本一样，除了最小化$x$时候加了扩增项。
$$begin{split}
x^{k+1} & = argmin_x L_{rho}(x, y^k) \
y^{k+1} & = y^k + rho(Ax^{k+1} - b) \
end{split}$$
这就是乘数法，可能也是因为更新对偶变量$y$时步长由原来变化的 $alpha^k$转为固定的$ρ$了吧。该算法在即使$f(x)$不是严格凸或者取值为$+infty$情况都可以成立，适用面更广。同样可以简单证明primal变量$x$和对偶变量$y$可以同时达到最优。
虽然增广拉格朗日方法有优势，但也破坏了对偶上升方法利用分解参数来并行的优势。当$f$是可分时，对于增广拉格朗日却是不可分的（因为平方项写成矩阵形式无法用之前那种分块形式），因此在 $x -min$步时候无法并行优化多个参数$x_i$。如何改进，继续下面的议题就可以慢慢发现改进思想的来源。
为了整合对偶上升可分解性与乘数法优秀的收敛性质，人们就又提出了改进形式的优化ADMM。目的就是想能分解原函数和扩增函数，以便于在对$f$更一般的假设条件下并行优化。
什么是ADMM算法
参考boyd大神的论文以及这篇lomyal的中文博文的介绍,我们先讲讲什么是ADMM算法。
ADMM的全称是交替方向乘子法（Alternating Direction Method of Multipliers），这是一种求解统计优化问题的计算框架, 适用于求解分布式凸优化问题。 ADMM 通过分解协调过程，将大的全局问题分解为多个较小、较容易求解的局部子问题，并通过协调子问题的解而得到大的全局问题的解。
由于我们的QAP模型就是一个完美的凸优化问题，所以本文将讨论我们如何实现分布式的QAP模型。
众所周知，并行计算是大数据处理/高性能计算的主要手段，然而很多时候并没有什么用。
常见的易于并行的算法有Bootstrap和交叉验证，而大多数优化问题由于迭代过程本身的特性都不容易并行执行。
现在ADMM包提供了以下示例：
Lasso & 弹性网（Elastic Net）
Dantzig Selector
基追踪（Basis Pursuit）
最小绝对偏差（Least Absolute Deviation）
借助于 Eigen包和 RcppEigen封装，ADMM包的核心部分是用高效的c++代码编写的。ADMM的计算性能聘美像 glmnet 这样的前沿R包，比大多数现有的求解器表现都更加优越。
安装
拉取怡轩大神的在GitHub分享的ADMM包
if(!require(devtools)){
install.packages("devtools")
}
if(!require(Rcpp)){
install.packages("Rcpp")
}
devtools::install_github("yixuan/ADMM")ADMM 依赖于 RcppEigen, rARPACK 和 ggplot2 包， install_github() 命令将安装前述依赖。
由于 yixuan大神用C++做底层，开发了这个ADMM包，所以以Mac OS X为例，我们还需要确保安装gcc，否则会遇到如下bug：
clang: error: no such file or directory: '/usr/local/lib/libfontconfig.a'
clang: error: no such file or directory: '/usr/local/lib/libreadline.a'解决方案是添加libfontconfig.a文件到对应目录下。
ADMM算法
ADMM包为许多流行的统计模型提供了接口。ADMM算法通常解决以下形式的问题：
$$begin{aligned}
text{minimize}quad & f(x)+g(z)\
text{subject to}quad & Ax + Bz = c
end{aligned}$$
其中$x$和$z$是向量、$A$和$B$是矩阵、$f$和$g$是凸函数。大多数的统计优化问题可以用这样的形式表达，ADMM算法可以用以下的迭代更新公式表示：
$$begin{align}
x^{k+1} & :=underset{x}{argmin}left(f(x)+frac{rho}{2}Vert Ax+Bz^{k}-c+y^{k}/rhoVert_{2}^{2}right)\
z^{k+1} & :=underset{z}{argmin}left(g(z)+frac{rho}{2}Vert Ax^{k+1}+Bz-c+y^{k}/rhoVert_{2}^{2}right)\
y^{k+1} & :=y^{k}+rho(Ax^{k+1}+Bz^{k+1}-c)
end{align}$$
$rho>0$ 表示步长参数,在参考资料里我们可以找到更多关于详情。
快速上手
本节将演示一些ADMM包的基础用法，我们通过举一些简单的例子来展示ADMM函数的主要使用方法。在下一节，我们将有进一步深入的讨论。
首先，我们为 Lasso, Elastic Net 和 Dantzig Selector 模型生成一些伪随机数据。
set.seed(123)
n = 100
p = 20
nonzero = 5
b = matrix(c(runif(nonzero), rep(0, p - nonzero)))
x = matrix(rnorm(n * p, mean = 1.2, sd = 2), n, p)
y = 5 + x %*% b + rnorm(n)和R中其他大多数建模函数不同，ADMM使用了引用类架构建立并且拟合模型。因此语法也是遵循面向对象编程的（OOP）风格。下面是模型表示的经典例子：
调用一个特定的模型函数来创建一个模型对象。
通过这个模型对象的成员函数设置必要的参数和选项。
调用模型拟合成员函数实际运行评估程序通过。
满足更多的任务需求,比如策划和预测。
首先, 函数调用时简单粗暴的: 只需要提供数据矩阵和响应向量作为参数。以下代码建立了三个类型一致的模型对象。
library(ADMM)
mod1 = admm_lasso(x, y)
mod2 = admm_enet(x, y)
mod3 = admm_dantzig(x, y)注意,在这个阶段不进行真正的计算。模型对象仅仅是描述模型的设置,这些设置可以通过调用的成员函数修改。函数语法类似ggplot2的“链式调用”
mod2$penalty(alpha = 0.3)
mod2$opts(maxit = 1000)
mod3$penalty(lambda_min_ratio = 0.01)以上命令设置Elastic Net模型点 $alpha$ 参数为0.3。在 Dantzig Selector 模型中调整了调优参数序列，设置迭代次数的上限为1000。
在设置了必要参数和选项后，模型可以通过$fit()成员函数来拟合，这就是实际的计算。
fit1 = mod1$fit()
fit2 = mod2$fit()
fit3 = mod3$fit()现在 $beta$ 向量被包含在了我们返回结果的 beta 字段中。我们也可以用结果对象中的 $plot() 成员函数画图。
print(fit1$beta[1:6, 1:6])library(ggplot2)
fit1$plot() %+% ggtitle("Solution Path for Lasso")
fit2$plot() %+% ggtitle("Solution Path for Enet")
fit3$plot() %+% ggtitle("Solution Path for Dantzig Selector")
ADMM包的一个吸引人的特性是,大多数建模功能是“可链接的”,在某种意义上,一个成员函数可以跟着另一个调用。
因此上面的命令可以简化为一些较短的代码:
admm_lasso(x, y)$fit()$plot()
mod2 = admm_enet(x, y)$penalty(alpha = 0.3)$opts(maxit = 1000)
mod2$fit()$plot()
fit3 = admm_dantzig(x, y)$penalty(lambda_min_ratio = 0.01)$fit()
fit3$plot()并行ADMM
将数据分块后，每次迭代所需的时间减少，但收敛所需的迭代次数增加，通信成本同样不可忽略。
收敛速度较慢（通过容忍更大的误差来减少迭代）
实际的并行效果尚待考量（速度不一定比串行更快，但可以解决单机跑不了全部数据的问题）
基于ADMM算法的二次规划
假设$f$是如下（凸）的二次函数
$$f(x) = \frac{1}{2}x^TPx + q^T x + r$$
$P$是对称的半正定矩阵$P in mathbf{S}^n_{+}$。这种形式问题也包含了$f$是线性或者常数的特殊情况。若$P + rho A^TA$可逆，那么$x-update$步求个导即有如下的显示解，是$v$的仿射函数
$$x^{+} = (P + \rho A^TA)^{-1}(\rho A^Tv - q)$$
因此在$x$-minnimiztion步只需要做两个矩阵运算即可，求逆与乘积，选用合适的线性运算库即可以得到不错的计算性能。当然还可以利用一些矩阵分解技巧，这个要看矩阵大小和稀疏程度。因为对于$Fx=g$，可以将$F = F_1F_2cdots F_k$，然后 $F_iz_i = z_{i-1}, z_1 = F_1^{-1}g,x= z_k$，这样会更节省计算时间。其他矩阵计算技巧，基本都是如何对矩阵大规模求解，利用矩阵的稀疏性、缓存分解等来提高性能。此处不赘述，有个很重要的求逆的定理很有用：
$$(P + \rho A^TA)^{-1} = P^{-1} - \rho P^{-1}A^T(I + \rho AP^{-1}A^T)^{-1}AP^{-1}$$
如果对于上述二次函数受限于某仿射集$x$-update步就更复杂些，如
$$f(x) = \frac{1}{2}x^T Px + q^T x + r \quad \textbf{dorm} \,f=\{x \| Fx = g\}$$
$x$-update还有个重要的KKT方程可用：
$$
begin{pmatrix}
P + rho I & F^T \
F & 0 \
end{pmatrix}
begin{pmatrix}
x^{k + 1}\
v \
end{pmatrix}
+
begin{pmatrix}
q - rho(z^k - u^k) \
-g \
end{pmatrix}
= 0
$$
参考文献
作为分享主义者(sharism)，本人所有互联网发布的图文均遵从CC版权，转载请保留作者信息并注明作者 Harry Zhu 的 FinanceR 专栏:https://segmentfault.com/blog...，如果涉及源代码请注明GitHub地址：https://github.com/harryprince。微信号: harryzhustudio
商业使用请联系作者。
**粗体** _斜体_ [链接](http://example.com) `代码` - 列表 > 引用。你还可以使用@来通知其他用户。