工程中非凸优化利器： Majorization-Minimization

新闻资讯 | 2024-07-01 15:16

我来给这个专栏除草啦！这次介绍一个算法Majorization-Minimization (MM)。MM可以作为一个理论框架，但是我们这次不涉及收敛性什么的，就说一下在工程当中如何去应用它，包括一些构造技巧以及在机器学习和通信中的应用的例子。

本文分成一下几个部分：

Part I: MM算法原理

Part II:用MM导出经典算法

Part II.A:泰勒展开法 -- DC Programming, Gradient Descent, Cubic Regularized Newton

Part II.B:凸不等式 -- Expectation Maximization

Part III: 实际应用例子

Part III.A: 机器学习 -- 稀疏优化

Part III.B: 通信 -- 智能反射表面、MIMO系统

如果有啥没有说对的欢迎指出来，因为能力有限，之前也没有用过MM，所以如果有哪里错了请一定指出来！

Part I: MM算法原理

MM主要的思想呢就是非凸函数太难最小化了，我就构造几个简单的子问题分而治之。原理就是下面这张图

假设 $x_0$ 在A点，函数f很难最小化，我们就构造一个好优化并且在它上面的函数u，我们找到u的最小值B点，在下一步 $x_1$ 移动到和B横坐标一样的C点，用表达式写出来呢，就是

$x_{k+1}=\ ext{argmin}_{x \\in \\mathcal{D}}u(x_k,x)\\\\$

其中 $\\mathcal{D}$ 是可行域。

为什么u要在f上面呢，因为我们想让函数值下降，如果u在f的上面，我们可以保证

$f(x_{k+1}) \\leq u(x_k,x_{k+1}) \\leq u(x_k,x_k)\\\\$

从图上看呢，就是C点比B点低比A点低。这么一来，目标函数的值就不会上升了，我们期望它能下降，原理就是这么简单。算法也很简单

Majorization-Minimization的名字就是这么来的，找一个在上面的函数u——Majorization，最小化函数u——Minimization。

Part II:从MM中导出经典算法

虽然MM原理简单，有不少经典的算法可以用MM构造出来，我们能构造出很多经典算法，我们借着这个介绍两种MM中构造u函数方法——泰勒展开和凸不等式。

Part II.A:泰勒展开法 -- DC Programming, Gradient Descent, Cubic Regularized Newton

构造出来的函数u必须是好优化的，那么我们就想咯，一次函数和二次函数挺好优化的，再加上u和f需要在A点相切，我们就用泰勒展开咯。这里的阶数是按照[0]里面定义的，可能有点儿奇怪。

一阶展开 -- Difference of Convex Programming

我们只泰勒展开一次项并且保证在函数上面，那就只有凹函数了，凹函数也就是反着的凸函数，它们长这个样子

我们怎么应用这个性质呢，考虑下面一个非凸优化

$\\min_{x \\in \\mathcal{D}}\緻brace{f(x)}_{\ ext{convex}}+ \緻brace{g(x)}_{\ ext{concave}}\\\\$

我们看到f(x)是凸的，如果g(x)也能变成凸的整个问题就是一个凸优化问题。考虑到g(x)的线性展开总是在g(x)的上面，那么我们就可以用一阶展开对g(x)进行Majorization，MM的算法如下，每次迭代解一个凸问题：

$x_{k+1}=\ ext{argmin}_{y \\in \\mathcal{D}}f(x) + g(x_k) + \\langle \ abla g(x_k), y-x_k\\rangle \\\\$

这也就是DC programming[1]以及如果第一项f(x)=0的话就是generalized power method[2]，如果 $\\mathcal{D}$ 是Stiefel manifold的话就是最近马毅老师力推的matching, stretching, and projection[3]。

二阶展开 -- Gradient descent

对于一般的函数，一阶展开并不会始终在函数的上方，那么要怎么办呢，那就加一个正数（二次项），使得一阶展开加上这个二次函数在函数的上方，假设原函数是 $x^3$ ，一阶展开和加了二次函数的图长这个样子

可以看到，在可见范围内，一阶展开始终在原函数的下方，并不满足要求，但是加了二次函数就满足要求了，MM算法如下

$x_{k+1}=\ ext{argmin}_{y \\in \\mathcal{D}}f(x_k) + \\langle \ abla f(x_k), y - x_k \\rangle + \\frac{1}{2\\eta}\\|y-x_k\\|_2^2\\\\$

这个东西又恰好是我们熟悉的梯度下降，推导过程也很简单，令 $g(y)=f(x_k) + \\langle \ abla f(x_k), y - x_k \\rangle + \\frac{1}{2\\eta}\\|y-x_k\\|_2^2$ ，求出g(y)导数等于零的点也就是

$g'(y)=\ abla f(x_k) + \\frac{1}{\\eta}(y-x_k)=0\\\\$

化简出来也就是 $y=f(x_k) - \\eta \ abla f(x_k)$ ，也就是我们熟悉的梯度下降方法。如果在目标函数上加一项不可微的函数，那就是近端梯度下降啦[4]。

三阶展开 -- Cubic Regularized Newton：一阶二阶都是著名算法，那二阶泰勒展开加个三阶项怎么样？这就变成了Cubic Regularized Newton[5]。它长这个样子

$x_{k+1}=\ ext{argmin}_{y \\in \\mathcal{D}}f(x_k) + \\langle \ abla f(x_k), y - x_k \\rangle + \\langle \ abla^2f(x_k)(y-x_k),y-x_k\\rangle + \\frac{L}{6}\\|y-x_k\\|^3 \\\\$

这个方法也是一个最近的热点，在这里推销一下专栏编辑 @Zeap 关于这个优化方法的文章[6.1][6.2]。

Part II.B:凸不等式 -- Expectation Maximization

还有一类就是用一些不等式放缩比如琴声不等式，可以用MM推出著名的EM算法。这个方面感觉文章[0]的ppt讲得特别好，这里就懒一下直接放截图啦。

Part III: 实际应用例子

Part III.A: 机器学习 -- 稀疏优化

我们知道如果要让一个东西变得稀疏，我们去最小化它的0范数，但是零范数又太难优化了，我们就去最小化一个凸函数1范数，其实1范数和零范数差得挺多的

相比较于1范数，这个 $\\phi(x)=\\log\\left(1+\\frac{|x|}{\\epsilon}\\right)$ 更加像0范数一点儿，所以它诱导稀疏的能力会更强一点儿。虽然它既不是凸的也不是凹的，但是可以注意到 $\\psi(x)=\\log\\left(1+\\frac{x}{\\epsilon}\\right)$ 是一个凹函数，一阶展开始终在它上面，我们就可以用上面所说的一阶展开的方法进行Majorization。

在这个例子里面，优化问题是 $\\min_{x \\in \\mathcal{D}}\\sum_{i=1}^n \\psi(|x^{(i)}|)$ ，注意到 $\\psi(x)=\\log\\left(1+\\frac{x}{\\epsilon}\\right)$ 是一个凹函数，按照前面DC Programming的方法，MM每次去最小化一个一阶展开

$\\begin{aligned}x_{k+1}&=\ ext{argmin}_{y \\in \\mathcal{D}}\\sum_{i=1}^n \\left( \\psi(|x_k^{(i)}|) + \\psi'(|x_k|)(|y|-|x_k|) \\right) \\\\&=\ ext{argmin}_{y \\in \\mathcal{D}}\\sum_{i=1}^n \\left( \\psi'(|x_k|)(|y|) \\right)=\ ext{argmin}_{y \\in \\mathcal{D}}\\sum_{i=1}^n \\left( \\frac{1}{|x_k^{(i)}|+\\epsilon}|y| \\right) \\end{aligned}$

其中 $x^{(i)}_k$ 表示第k轮迭代的变量的第i个分量。这个就是著名的reweighted l1 norm minimization[7]。效果上比直接最小化1范数好了很多呢（a是ground-truth，c是1范数，d是reweighted l1）

Part III.B: 通信 -- 智能反射表面

我在写这篇文章的时候，想找一个机器学习以外的比较热应用，然后找了好几天没找到，就去看去年通信和信号处理会的best paper里看有没有，终于找到了一篇用MM的[8](MM算法和详细推导在paper的第四页)。具体背景我也看不懂，大概是要最大化一个类似于SNR的比值，直接看优化问题吧

$\\begin{aligned}&\\max_{v \\in \\mathbb{C}^n}&& \\frac{v^HY_1v}{v^HY_ev}\\\\ &\ ext{subject to}&& |v_k|=1 \\end{aligned}\\\\$

其中 $Y_1, Y_e$ 都是正定矩阵，那个约束条件是因为那个v是个反射角，有 $v_k=e^{i\ heta_k}$ 所以会有那个模长是1的constraint。因为是最大化，所以在MM过程当中想找一个函数从始终在原函数的下面。令 $y=v^H Y_e v$ ，可以看到目标函数对于 $(v,y)$ 是凸的，一次展开式在原函数的下面，对 $(v,y)$ 展开得到

$\\frac{v^HY_1v}{y}\\geq \\frac{v_z^HY_1v_z}{y}+ \\frac{\\partial \\left( \\frac{v^HY_1v}{y}\\right) }{\\partial v}(v - v_z) + \\frac{\\partial \\left( \\frac{v^HY_1v}{y}\\right) }{\\partial y}(y - y_z) \\\\=\\frac{2\ ext{Re}(v^H_zY_1v)}{v_z^HY_ev_z}- \\frac{v_z^HY_1v_z}{(v_z^HY_ev_z)^2}\緻brace{v^HY_ev}_{\\mathcal{P}_1}\\\\$

对于模长是1的constraint来说，第一项是一次项很好优化，第二项是二次项有点儿困难，所以我们继续对 $\\mathcal{P}_1$ 进行Majorization（找一个在它上面的函数）。对于厄米矩阵 $L \\preceq M$ 我们有

$\\begin{aligned}x^HLx &=x_0^H L x_0 + 2\ ext{Re}((x-x_0)^H L x_0) + (x - x_0)^H L(x-x_0) \\\\ &\\leq x_0^H L x_0 + 2\ ext{Re}((x-x_0)^H L x_0) + (x - x_0)^H M(x-x_0) \\\\ &=x^HMx + 2\ ext{Re}(x^H(L-M)x_0) + x_0^H(M-L)x_0 \\end{aligned}\\\\$

取 $L=Y_e, M=\\lambda_{\ ext{max}}(Y_e)$ 我们就可以接着推导

$\\begin{aligned}\\frac{v^HY_1v}{v^HY_ev}&\\geq \\frac{2\ ext{Re}(v^H_zY_1v)}{v_z^HY_ev_z}- \\frac{v_z^HY_1v_z}{(v_z^HY_ev_z)^2}v^HY_ev\\\\ &\\geq \\frac{2\ ext{Re}(v^H_zY_1v)}{v_z^HY_ev_z}- \\frac{v_z^HY_1v_z}{(v_z^HY_ev_z)^2}\\left\\{ v^H \\lambda_{\ ext{max}}(Y_e)v + 2\ ext{Re}(v_z^H[Y_e - \\lambda_{\ ext{max}}(Y_e)I]v) + v_z^H[\\lambda_{\ ext{max}}(Y_e)I -Y_e]v_z\\right\\}\\\\ &=\\frac{2\ ext{Re}(v^H_zY_1v)}{v_z^HY_ev_z}- \\frac{v_z^HY_1v_z}{(v_z^HY_ev_z)^2}\\left\\{\\lambda_{\ ext{max}}(Y_e)n + 2\ ext{Re}(v_z^H[Y_e - \\lambda_{\ ext{max}}(Y_e)I]v) + v_z^H[\\lambda_{\ ext{max}}(Y_e)I -Y_e]v_z\\right\\}=g(v,v_z) \\end{aligned}$