经典的SDR算法： 用半正定松弛法 ( Semidefinite Relaxation) 求解二次优化问题

前言

本文是博主对于 Zhi-quan Luo 老师的经典著作 《Semidefinite Relaxation of Quadratic Optimization Problems》 的读书笔记，希望可作为对全文以中文形式的核心梳理。

单刀直入

首先， Semidefinite Relaxation (SDR) 适用的问题可以写为如下形式：

$$\begin{array}{rl}\underset{x\in {\mathbb{R}}^n}{\min }& x^{T}Cx\\ \text{ s.t. }& x^T{A}_{i}x{⊵}_i{b}_i,i=1,\dots ,m\end{array}\text{\hspace{1em}(1)}$$

有几个注意点：

• 暂时定义在实数域上， 但复数域有类似的结论。
• ${⊵}_i$ 表示 $\ge$, $\le$, $=$ 的 其中之一。也就是说， 只要限制条件是二次约束即可。
• $C$， $A_i$ 都是对称矩阵。 即, $C,A_1,\dots ,A_m\in {\mathbb{S}}^n$。

显然， 目标函数和限制条件中都有两个$x$， 因此这类问题也被称为 二次约束二次规划 (quadratically
constrained quadratic programs, QCQP) 问题。 需要注意的是， 这并不是一个凸问题。 因为 $x^T{A}_{i}x>0$ 和 $x^T{A}_{i}x=0$ 都并非一个凸集。

而SDR方法的目的就是将问题转化为一个凸问题从而求解。

我们引入一个新的变量 $X=x{x}^T$， 并以此为优化变量， 那么问题(1)可以被写为：
$$\begin{array}{rl}\underset{X\in {\mathbb{S}}^n}{\min }& \mathrm{Tr}(CX)\\ \text{ s.t. }& \mathrm{Tr}\left(A_{i}X\right){⊵}_i{b}_i,i=1,\dots ,m,\\ & X⪰0,\mathrm{rank}(X)=1\end{array}\text{\hspace{1em}(2)}$$
详细解释一下这一步的转化：

• 注意到 (1) 的 目标函数 $x^{T}Cx$显然是一个标量。 而标量可以看成一个 $1\times 1$的矩阵， 且迹等于本身。 即 $\mathrm{t}\mathrm{r}(x^{T}Cx)=x^{T}Cx$。 而根据迹的性质， 有 $\mathrm{t}\mathrm{r}(AB)=\mathrm{t}\mathrm{r}(BA)$。 也就是说， $\mathrm{t}\mathrm{r}(x^{T}Cx)=\mathrm{t}\mathrm{r}(Cx{x}^T)=\mathrm{t}\mathrm{r}(CX)$。 这也就变成了 (2)中的目标函数。 那么 (2) 中的限制条件，也就可以根据相同的变换轻易得到。
• 相比于(1)， (2)中多了两个限制条件， 这是由于 $X=x{x}^T$ 引入的。 显然$\mathrm{rank}(X)=1$， 因为$X$的所有列都线性相关。 也因此，只有一个特征值等于 $x^{T}x$，显然是大于0的，所以肯定是半正定矩阵。

这一步转化的意义在哪呢？

事实上， 除了 $\mathrm{rank}(X)=1$ 这个限制以外， 整个问题是一个凸问题。 由于迹函数是个仿射变换， 很容易证明其凸性。 而半正定约束也显然是个凸集， 因为两个半正定矩阵的和一定也是一个半正定矩阵。

因此， 工程上很自然的想法就是， **暂时得忽略 该约束$\mathrm{rank}(X)=1$， 并通过凸优化方法 （事实上以CVX工具包为代表的内点法数值算法）来求解出一个解 $X$. 这一步大家百度下CVX的使用就可以了，没有任何障碍。

而问题的关键点在于，始终我们问题需要得到的是$x$，而非$X$。 如果通过凸优化求出的$X$刚好就是一个$\mathrm{rank}(X)=1$的矩阵，那么$x$可以直接得到（对X做个SVD或EVD啥的都可以）。 但更多时候， $X$没有道理是一个 $\mathrm{rank}(X)=1$的矩阵。 也就是说， 根本不存在一个$x$， 使得 $X=x{x}^T$。

工程实现

那么SDR要解决的另一个问题就是如何从$X$中恢复出一个满足原问题约束的$x$。
出于工程实现角度的思路其实非常清晰， 我们可以直接求解下面的问题以获得 $x$：
$$\underset{x}{\min }\Vert X-x{x}^T{\Vert }_F^2\text{\hspace{1em}(3)}$$
即找一个最接近$X$的 $x{x}^T$。 这个问题很经典了，而他的闭式解就是， $x$是$X$的最大特征向量乘以最大特征值的平方根。 即， 将$X$的特征值分解表示为：
$$X^{\star }=\sum _{i=1}^r{\lambda }_i{q}_i{q}_i^T$$
其中${\lambda }_i$， $q_i$，代表第$i$大的特征值和对应的特征向量。 那么(3)的最优解 $x$就是 $\tilde{x}=\sqrt{{\lambda }_1}{q}_1$。

但是，问题又出现了！ 虽然$X$满足了(2)中的所有的约束， 但由于在这一步中，(3)中求得的$x$是一个近似的结果。 此时得到的$x$，不一定满足原问题(1)中的约束了！ 也就是说， 经过一顿转化之后， 求得的$x$并不是原问题的可行解！此时怎么办呢？ 解决方案也非常符合工程思维：直接从$x$出发找一个离$x$最近的可行解， 这就是SDR的最终结果。

最后提两点：

• 直接从$x$出发找一个离$x$最近的可行解， 这本身就是一个独立的问题， 具体问题具体分析，后面会给一个具体的例子。
• 毫无疑问，在SDR中有许许多多的近似，工程处理。 因此， SDR的解远非全局最优解。

具体实例

为节约篇幅， 本文只讲一个例子。 也不是论文中给出的， 而是HBF混合波束成形的经典论文 《Alternating Minimization Algorithms for Hybrid Precoding in Millimeter Wave MIMO Systems》中的 一个使用。 有兴趣的读者也可以去看 SDR论文的原文， 接触更多的例子。

以HBF为例， 作者将问题转化为了如下的形式， 熟悉HBF的读者应该很亲切了：

$$\begin{array}{ll}\underset{{\mathbf{F}}_{\mathrm{B}\mathrm{B}}}{\mathrm{minimize}}& {\Vert {\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}-{\mathbf{F}}_{\mathrm{R}\mathrm{F}}{\mathbf{F}}_{\mathrm{B}\mathrm{B}}\Vert }_F^2\\ \text{ subject to }& {\Vert {\mathbf{F}}_{\mathrm{B}\mathrm{B}}\Vert }_F^2=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}\end{array}\text{\hspace{1em}(4)}$$

作者使用了 SDR作为一种算法来求解该问题。 那么首先， 为什么会选用 SDR呢？ 因为可以看到， 目标函数和限制条件都是二次形式， 这就是我们之前说到的SDR目标问题：QCQP问题。 接下来，通过这个例子展示如果实战使用SDR。

首先， 刚刚介绍SDR时都是向量形式的变量$x$, 而(4)中都是矩阵，不利于处理， 因此第一步先进行列化：
$$\begin{array}{rl}{\Vert {\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}-{\mathbf{F}}_{\mathrm{R}\mathrm{F}}{\mathbf{F}}_{\mathrm{B}\mathrm{B}}\Vert }_F^2& ={\Vert \mathrm{vec}\left({\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}-{\mathbf{F}}_{\mathrm{R}\mathrm{F}}{\mathbf{F}}_{\mathrm{B}\mathrm{B}}\right)\Vert }_2^2\\ & ={\Vert \mathrm{vec}\left({\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}\right)-\mathrm{vec}\left({\mathbf{F}}_{\mathrm{R}\mathrm{F}}{\mathbf{F}}_{\mathrm{B}\mathrm{B}}\right)\Vert }_2^2\\ & ={\Vert \mathrm{vec}\left({\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}\right)-\left({\mathbf{I}}_{N_s}\otimes {\mathbf{F}}_{\mathrm{R}\mathrm{F}}\right)\mathrm{vec}\left({\mathbf{F}}_{\mathrm{B}\mathrm{B}}\right)\Vert }_2^2.\end{array}$$
进行变量代换后 ($\mathbf{f}=\mathrm{vec}\left({\mathbf{F}}_{\mathrm{o}\mathrm{p}\mathrm{t}}\right)$,$\mathbf{b}=\mathrm{vec}\left({\mathbf{F}}_{\mathrm{B}\mathrm{B}}\right)$ and $\mathbf{E}={\mathbf{I}}_{N_s}\otimes {\mathbf{F}}_{\mathrm{R}\mathrm{F}}.$，原问题可写为：

$$\begin{array}{rlr}& \underset{\mathbf{b}}{\mathrm{minimize}}\Vert t\mathbf{f}-\mathbf{E}\mathbf{b}{\Vert }_2^2& \\ & \text{ subject to }& \{\begin{array}{l}\Vert \mathbf{b}{\Vert }_2^2=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}\\ t^2=1\end{array}\end{array}\text{\hspace{1em}(5)}$$

**这里必须强调的一点是， 作者引入了 $t$这样一个辅助变量。这是SDR的常见技巧之一！**接下来阐释这个$t$的作用。

• 首先， 对于原问题的影响上：以(5)求解之后， 如果得到$t=1$， 那么$\mathbf{f}$就是(4）的最优解。 而如果$t=-1$， 那么$\mathbf{f}$就是(4)的最优解的相反数， 仍等效于得到了最优解。
• 其次， 加入这个$t$辅助变量后， 可以更好地将(5)写成SDR的标准形式。 继续看：

目标函数可以写为：

$$\Vert t\mathbf{f}-\mathbf{E}\mathbf{b}{\Vert }_2^2=\left[\begin{array}{ll}{\mathbf{b}}^H& t\end{array}\right]\left[\begin{array}{cc}{\mathbf{E}}^H\mathbf{E}& -{\mathbf{E}}^H\mathbf{f}\\ -{\mathbf{f}}^H\mathbf{E}& {\mathbf{f}}^H\mathbf{f}\end{array}\right]\left[\begin{array}{l}\mathbf{b}\\ t\end{array}\right]$$
注意到： 此时令 $\mathbf{y}=\left[\begin{array}{l}\mathbf{b}\\ t\end{array}\right]$, $\mathbf{C}=\left[\begin{array}{cc}{\mathbf{E}}^H\mathbf{E}& -{\mathbf{E}}^H\mathbf{f}\\ -{\mathbf{f}}^H\mathbf{E}& {\mathbf{f}}^H\mathbf{f}\end{array}\right]$, 上式就等于 ${\mathbf{y}}^H\mathbf{C}\mathbf{y}$，而$\mathbf{C}$是个显然的半正定矩阵：$\mathbf{C}=[\mathbf{E},-\mathbf{f}{]}^H[\mathbf{E},-\mathbf{f}]$，一个矩阵的共轭转置乘以本身肯定是半正定矩阵（可从特征值角度得证）。

那么(5)中的另外两个约束又可以转换为：
$$\begin{array}{rl}\Vert \mathbf{b}{\Vert }_2^2& =\left[\begin{array}{ll}{\mathbf{b}}^H& t\end{array}\right]\left[\begin{array}{cc}{\mathbf{I}}_{N_{\mathrm{R}\mathrm{F}}^t}{N}_s& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}\mathbf{b}\\ t\end{array}\right]=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}\\ t^2& =\left[\begin{array}{ll}{\mathbf{b}}^H& t\end{array}\right]\left[\begin{array}{cc}{0}_{N_{\mathrm{R}\mathrm{F}}^t}{N}_s& 0\\ 0& 1\end{array}\right]\left[\begin{array}{l}\mathbf{b}\\ t\end{array}\right]=1\end{array}$$

大家可以细致体会下， 为什么SDR可以应用广泛的原因， 因为大部分二次约束都可以通过稍加构造就写成QCQP的形式。

那么，最后经过一步变量代换：
$$\begin{array}{rl}& \mathbf{y}=\left[\begin{array}{c}\mathbf{b}\\ t\end{array}\right],\mathbf{Y}=\mathbf{y}{\mathbf{y}}^H,\mathbf{C}=\left[\begin{array}{cc}{\mathbf{E}}^H\mathbf{E}& -{\mathbf{E}}^H\mathbf{f}\\ -{\mathbf{f}}^H\mathbf{E}& {\mathbf{f}}^H\mathbf{f}\end{array}\right]\\ & {\mathbf{A}}_1=\left[\begin{array}{cc}{\mathbf{I}}_{N_{\mathrm{R}\mathrm{F}}^t{N}_s}& 0\\ 0& 0\end{array}\right],{\mathbf{A}}_2=\left[\begin{array}{cc}{0}_{N_{\mathrm{R}\mathrm{F}}^t}{N}_s& 0\\ 0& 1\end{array}\right]\end{array}$$

（5）终于可写为SDR的标准形式：
$$\begin{array}{rl}& \underset{\mathbf{Y}\in {\mathbb{H}}^n}{\mathrm{minimize}}\underset{\mathrm{Tr}}{(}\mathbf{C}\mathbf{Y})\\ & \text{ subject to }\{\begin{array}{l}\mathrm{Tr}\left({\mathbf{A}}_1\mathbf{Y}\right)=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}\\ \mathrm{Tr}\left({\mathbf{A}}_2\mathbf{Y}\right)=1\\ \mathbf{Y}⪰0,\mathrm{rank}(\mathbf{Y})=1\end{array}\end{array}$$
通过松弛掉最后一个约束， 问题就可以通过凸优化直接求解。
最后必须强调， 如何从$\mathbf{Y}$中恢复出$\mathbf{b}$（即原问题真正的目标）。 首先，用特征值分解得到$\mathbf{y}$。 此时 $\mathbf{y}$的最后一个元素大概率不是1，而按照定义案例$t$应该必须是$1或-1$。此时处理方式很简单， 直接对最后一个元素取正负号， 强行归成$1$或$-1$即可。 而$\mathbf{b}$（ $\mathbf{y}$的前面的元素）也不符合$\Vert \mathbf{b}{\Vert }_2^2=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}$的约束，此时的做法就是将$\mathbf{b}$强行归一化成$\Vert \mathbf{b}{\Vert }_2^2=\frac{N_{\mathrm{R}\mathrm{F}}^t{N}_s}{N_t}$即可。 操作非常简单， 那么其对应的思路就是我们在上文中提到的：当得到了$x$后， 以离$x$最近的可行解作为原问题的解。

实例2

再举个近期的实例好了。 是来自于 智能反射面相关paper 《Intelligent Reflecting Surface Enhanced Wireless
Network: Joint Active and Passive Beamforming Design》中， 对于passive beamforming 的设计， 同样使用了 SDR 算法。

原问题可以转化为这样的形式：
$$\begin{array}{cl}{\max }_{\mathbf{v}}& {\mathbf{v}}^H\Phi {\Phi }^H\mathbf{v}+{\mathbf{v}}^H\Phi {\mathbf{h}}_d+{\mathbf{h}}_d^H{\Phi }^H\mathbf{v}\\ \text{ s.t. }& \left|v_n\right|=1,\forall n=1,\cdots ,N\end{array}$$
简单而言， 就是其余参数都已给定， 要求优化向量$\mathbf{v}$， 其中$\mathbf{v}$的每个元素需要满足恒模约束。
那么同样的， 可以通过引入辅助变量 $t$， 把问题转换为 齐次的 QCQP 问题：

$$\begin{array}{cl}{\max }_{\overline{\mathbf{v}}}& {\overline{\mathbf{v}}}^H\mathbf{R}\overline{\mathbf{v}}\\ \text{ s.t. }& \left|{\bar{v}}_n\right|=1,\forall n=1,\cdots ,N+1\end{array}$$

其中，
$$\mathbf{R}=\left[\begin{array}{cc}\mathbf{\Phi }{\Phi }^H& \mathbf{\Phi }{\mathbf{h}}_d\\ {\mathbf{h}}_d^H{\mathbf{\Phi }}^H& 0\end{array}\right],\overline{\mathbf{v}}=\left[\begin{array}{l}v\\ t\end{array}\right]$$

那么根据SDR的方法， 令 $\mathbf{V}=\bar{\mathbf{v}}{\bar{\mathbf{v}}}^H$, 就可以将问题转化为：

$$\begin{array}{ll}{\max }_{\mathbf{V}}& \mathrm{tr}(\mathbf{R}\mathbf{V})\\ \text{ s.t. }& {\mathbf{V}}_{n,n}=1,\forall n=1,\cdots ,N+1,\\ & \mathbf{V}⪰0\end{array}$$

当然这里同样的放松掉了$\mathrm{rank}(\mathbf{V})=1$这一约束。 可以看到， 上述问题已经是凸问题了， 因为非凸的恒模约束经过转化后变成了 ${\mathbf{V}}_{n,n}=1,\forall n=1,\cdots ,N+1$， 这无疑是一个凸约束 （可以用凸集的定义轻松得证）。

因此， SDR也是处理恒模约束的一种思路。

注意， 这个例子的后续处理，也是从SDR解$\mathbf{V}$获得原问题可行解的一种经典思路：SDR + Gaussian randomization. 这种处理相比于上面提到的直接寻找欧氏距离最近的$\mathbf{v}$，会有通过复杂度提升换来的一定性能增益。

首先， 对$\mathbf{V}$进行 EVD分解， 得到$\mathbf{V}=\mathbf{U}\mathbf{\Sigma }{\mathbf{U}}^H$， 此时， 我们令 $\overline{\mathbf{v}}=\mathbf{U}{\mathbf{\Sigma }}^{1/2}\mathbf{r}$， 其中$\mathbf{r}\in \mathcal{C}\mathcal{N}\left(0,{\mathbf{I}}_{N+1}\right)$。 注意， 如果是上面提到的普通SDR方法， $\mathbf{r}=[1,0,\cdots ,0]$， 即$\bar{\mathbf{v}}$等于最大特征向量。 而这里，引入了一个随机变量， 使得我们尝试更多种的$\bar{\mathbf{v}}$， 再选出使得目标函数值最大的$\bar{\mathbf{v}}$作为最优解。 可以看到， 这种方法其实就是通过遍历更多解来获得更好的性能。 最后， 还是把得到的$\bar{\mathbf{v}}$转换为可行解： $\mathbf{v}=e^{j\arg \left({\left[\frac{\bar{v}}{{\bar{v}}_{N+1}}\right]}_{(1:N)}\right)}$。

代码

关于SDR 更细节一些的实现和相关代码， 在下一篇博文 经典的SDR算法（下）：SDR的具体使用细节与相关代码 进行了更详细的介绍。