斜风细雨作小寒，淡烟疏柳媚晴滩。入淮清洛渐漫漫。
雪沫乳花浮午盏，蓼茸蒿笋试春盘。人间有味是清欢。
---- 苏轼
更多精彩内容请关注微信公众号 “优化与算法”
迭代硬阈值（Iterative Hard Thresholding）算法是求解基于 \({\ell _0}\) 范数非凸优化问题的重要方法之一，在稀疏估计和压缩感知重构等领域应用较多。IHT最初由Blumensath, Thomas等人提出，后来许多学者在IHT算法的基础上不断发展出一些改进算法，如正规化迭代硬阈值算法（Normalized Iterative Hard Thresholding, NIHT)、共轭梯度硬阈值迭代算法（Conjugate Gradient Iterative Hard Thresholding, CGIHT）、基于回溯的硬阈值迭代算法（Backtracking based Iterative Hard Thresholding）等，以及后来与贪婪算法结合产生的一些算法，如硬阈值追踪算法（Hard Thresholding Pursuit, HTP）就是结合子空间追踪算法（SP）和IHT算法的优势得到的。
笔者在此将上述这些算法原理做简要分析，并附上算法和测试代码供参考。
以稀疏近似问题为例，介绍IHT算法。对于以下优化问题：
\[\mathop {\min }\limits_{\bf{x}} \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2{\rm{ }} ~~~~~~s.t.~~~~~{\rm{ }}{\left\| {\bf{x}} \right\|_0} \le S~~~~~~~~~~~~(1)\]
其中 \({\bf{y}} \in {{\bf{R}}^M}\), \({\bf{A}} \in {{\bf{R}}^{M \times N}}\)，且 \(M < N\)， \({\left\| {\bf{x}} \right\|_0}\) 为 \(\bf{x}\) 的 \({\ell _0}\) 范数，表示 \(\bf{x}\) 中非零元素的个数。由于 \(M < N\)， 显然（1）式是一个欠定问题，加上 \({\ell _0}\) 约束，式（1）同时是个组合优化问题，无法在多项式时间复杂度内求解。
为了求解优化问题（1），Blumensath 等人提出了著名的 IHT 算法，其迭代格式为：
\[{{\bf{x}}_{k + 1}} = {H_k}({{\bf{x}}_k} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }({\bf{y}} - {\bf{A}}{{\bf{x}}_k}))~~~~~~~~~~~~~~~~(2)\]
其中 \({H_S}( \cdot )\) 表示取 \(\cdot\) 中最大的 \(S\) 个非零元素。
式(2)看起来并不复杂，直观上来看就是对（1）式用常数步长的梯度下降法来求解，其中在每次迭代中对更新的结果进行一次硬阈值操作。为什么硬阈值操作+梯度下降可以求解(1)式的 \({\ell _0}\) 范数优化问题呢？文献[1]中作者利用Majorization-Minimization优化框架推导出了（2）式所示迭代公式，并证明了IHT算法在 \({\left\| {\bf{A}} \right\|_2} < 1\) 时，算法能收敛到局部最小值。以下对IHT算法的原理做简要分析。
1.2 Majorization-Minimization优化框架首先简单介绍一下MM优化框架，实际上MM优化框架是处理非凸优化问题中常见的一种方法，其形式简单，思想优雅。后期关于MM思想再以专题形式进行详细的总结和讨论。
(1) 式的目标函数为
\[{F}({\bf{x}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 ~~~~~~~~~~~~~~~~~~~~~~~~~~~(3)\]
(1) 无法得到一个闭式解，只能用迭代方法来求得一个近似解。
MM优化思想：当目标函数 \(F({\bf{x}})\) 比较难优化时，通常可以选择更容易优化的替代目标函数 \(G({\bf{x}})\)，当 \(G({\bf{x}})\) 满足一定的条件时，\(G({\bf{x}})\) 的最优解能够无限逼近原目标函数 \(F({\bf{x}})\) 的最优解。这样的 \(G({\bf{x}})\) 需要满足3个条件：
容易优化
对所有 \(\bf{x}\)，\(G({\bf{x}}) \ge F({\bf{x_{k}}})\)
\({G_k}({{\bf{x}}_k}) = F({{\bf{x}}_k})\)
在满足 \({\left\| {\bf{A}} \right\|_2} < 1\) 条件时，目标函数 (4) 可用如下替代目标函数来代替：
\[G({\bf{x}},{\bf{z}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 - \left\| {{\bf{Ax}} - {\bf{Az}}} \right\|_2^2 + \left\| {{\bf{x}} - {\bf{z}}} \right\|_2^2~~~~~~~(4)\]
从 (4) 式可以看出，替代目标函数 \(G({\bf{x}},{\bf{z}})\) 是满足MM优化框架中所列条件的。
对式 (4) 展开并化简：
\[\eqalign{
& G({\bf{x}},{\bf{z}}) = \left\| {{\bf{y}} - {\bf{Ax}}} \right\|_2^2 - \left\| {{\bf{Ax}} - {\bf{Az}}} \right\|_2^2 + \left\| {{\bf{x}} - {\bf{z}}} \right\|_2^2 \cr
& = \left( {\left\| {\bf{y}} \right\|_2^2 - 2{{\bf{x}}^{\mathop{\rm T}\nolimits} }{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} + \left\| {{\bf{Ax}}} \right\|_2^2} \right) - \left( {\left\| {{\bf{Ax}}} \right\|_2^2 - 2{{\bf{x}}^{\mathop{\rm T}\nolimits} }{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}} + \left\| {{\bf{Az}}} \right\|_2^2} \right) + \left( {\left\| {\bf{x}} \right\|_2^2 - 2{{\bf{x}}^T}{\bf{z}} + \left\| {\bf{z}} \right\|_2^2} \right) \cr
& = \left[ {\left\| {\bf{x}} \right\|_2^2 - \left( {2{{\bf{x}}^T}{\bf{z}} + 2{{\bf{x}}^T}{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} - 2{{\bf{x}}^T}{{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}}} \right)} \right] + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr
& = \left[ {\left\| {\bf{x}} \right\|_2^2 - 2{{\bf{x}}^T}\left( {{\bf{z}} + {{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{y}} - {{\bf{A}}^{\mathop{\rm T}\nolimits} }{\bf{Az}}} \right)} \right] + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr
& = \sum\limits_i {\left[ {{\bf{x}}_i^2 - 2{{\bf{x}}_i}\left( {{{\bf{z}}_i} + {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{y}} - {\bf{A}}_i^{\mathop{\rm T}\nolimits} {\bf{Az}}} \right)} \right]} + \left\| {\bf{y}} \right\|_2^2 + \left\| {\bf{z}} \right\|_2^2 + \left\| {{\bf{Az}}} \right\|_2^2 \cr} \]