深入分析：Lasso问题和原子范数问题研究

写在前面

本文将主要围绕Lasso问题和原子范数等经典问题进行对偶问题的推导、分析，由于笔者的数理基础浅薄，下面的证明过程若存在错误，欢迎评论指正。

LASSO问题

推导

​ 问题定义：\(\underset{x}{\min}||y-X\beta||_2^2+\tau||\beta||_\mathcal{A}\)

​ 问题推导：

​ 0、上述问题是典型的无约束问题，可以通过变量替换的思想进行处理。

​ 1、令\(z=X\beta\)，上述问题更新为\(g(u)=\min_{\beta,z}=\frac{1}{2}\Vert y-z\Vert_2^2+\lambda\Vert\beta\Vert_1+u^T(z-X\beta)\).

​ 2、可以观察到\(g(u)\)中关于\(\beta\)和\(z\)的元素项不存在耦合关系，因此可进步将\(g(u)\)问题拆解为独立的最小项\(g(\beta)\)和\(g(z)\)，其中\(g(\beta)=\min_\beta~\lambda\Vert\beta\Vert_1-u^TX\beta\), \(g(z)=\min_\beta~\frac{1}{2} \Vert y-z\Vert_2^2+u^Tz\)

​ 3、\(g(\beta)=\min_\beta~\lambda\Vert\beta\Vert_1-u^TX\beta=-\max_\beta~\lambda(-\Vert\beta\Vert_1+\frac{u^TX\beta}{\lambda})=-\lambda~g^*(\frac{X^Tu}{\lambda})=-\lambda I_{\{v_i|\Vert v\Vert_\infty \leq1\}}(\frac{X^Tu}{\lambda})\),这个最小项可以表征为示性函数形式，示性函数\(f^*(y)=\begin{cases}0&||y||_*\leq1\\ \infty&\text{otherwise}\end{cases}=I_{\{z:||z||_*\leq1\}}(y)\).

​ 4、对\(g(z)\)求极值，可以得到\(-(y-z)+u=0\)，即\(z=y-u\).

​ 5、将上述约束代入\(g(z)\)，可以得到下式：

​ \(g(z)=\frac{1}{2}\Vert u\Vert^2_2+u^T(y-u)=u^Ty-\frac{1}{2}\Vert u\Vert_2^2=-\frac{1}{2}[\Vert y\Vert_2^2+\Vert u\Vert_2^2-2u^Ty]+\frac{1}{2}\Vert y\Vert_2^2=-\frac{1}{2}\Vert y-u\Vert_2^2+\frac{1}{2} \Vert y\Vert_2^2\)

​ 那么对偶问题可以表示为如下形式：

​ \(\max_u -\frac{1}{2}\Vert y-u\Vert_2^2+\frac{1}{2} \Vert y\Vert_2^2 ~s.t.X^Tu\leq\lambda\)

原子范数对偶问题

推导

​ 有噪声情况下，原子范数的原问题可以抽象为：\(\operatorname*{min}_{x}\|x\|_{\mathcal{A}}\mathrm{subject~to}\left\{\begin{array}{l}{{\mathbf{y}={\mathcal{F}}_{M}x+\mathbf{n},}}\\ {{\|\mathbf{n}\|_{2}\leq\epsilon.}}\end{array}\right.\)

​ 对偶函数可以写为\(g(x,u)\)在\(C\)上的下确界，即\(g(\mathbf{c},\xi)=\inf~L(x,\mathbf{c},\xi)\)

​ 下面对原问题的对偶问题进行推导：

​ 1、原问题的增广拉格朗日目标函数可以表示为：\(L(x,\mathbf{c},\xi)=\|x\|_{\mathcal{A}}+\mathrm{Re}\left[\mathbf{c}^{H}\left(\mathbf{y}-\mathcal{F}_{M}x-\mathbf{n}\right)\right]+\xi\left(\mathbf{n}^{H}\mathbf{n}-\epsilon^{2}\right)\)

将拉格朗日方程进行重写，\(\inf~L(x,\mathbf{c},\xi)=\mathrm{Re}[c^Hy-c^Hn]+\xi[n^Hn-\epsilon^2]+\inf[\Vert x\Vert_A-\mathrm{Re}[\lambda^HF_Mx]]\\\)

​ 2、下确界的求解是关于\(x\)的最小化，因此对原拉格朗日增广函数的最小化可以转换为对\([\Vert x\Vert_A-\mathrm{Re}[\lambda^HF_Mx]]\)求下确界。在求这项下确界时，需要对式中的噪声功率参数\(n\)和对偶变量\(\xi\)求偏导寻找极值点。

​ 当对噪声功率参数\(n\)求偏导时[目的是为了使噪声功率最小化]，有\(\frac{\partial g(\mathbf{c},\xi)}{\partial c}=-c+2\xi n=0\)，可以得到最佳极值点\(n_o=\frac{c}{2\xi}\)，此时对应的对偶函数为\(\begin{array}{c}g(\mathbf{c},\xi)|_{\mathbf{n_\circ}}=\mathrm{Re}\left[\mathbf{c}^{H}\mathbf{y}\right]-\dfrac{\mathbf{c}^{H}\mathbf{c}}{2\xi}+\xi\left(\dfrac{\mathbf{c}^{H}\mathbf{c}}{4\xi^{2}}-\epsilon^{2}\right)+\+\inf_{x}\left(\|x\|_{\mathcal{A}}-\mathrm{Re}\left[\mathbf{c}^{H}\mathcal{F}_{M}x\right]\right).\end{array}\)

​ 当对对偶变量\(\xi\)求偏导时，有\(\frac{\partial g(\mathbf{c},\xi)_{n_o}}{\partial \xi}=\frac{c^Hc}{4\xi^2}-\epsilon^2=0\),可以得到最佳极值点\(\xi_0=\frac{\Vert c\Vert}{2\epsilon}\).

​ 最后，基于最优极值点对偶函数可以表示为\(g(c)|_{n_o,\xi_o}=\mathrm{Re}[c^Hy-\epsilon\Vert c\Vert_2+\inf_x(\Vert x\Vert_A-\mathrm{Re}[c^H\mathcal{F}_{M}]x)]\).

​ 对于下确界项，对每个\(x_i\)，有\(\mathrm{Re}[(c^H\mathcal{F}_{M})_ix_i]=\mathrm{Re}[(\mathcal{F}_{M}^Hc)^H_ix_i]=|(\mathcal{F}_{M}^Hc)_i||x_i|\cos\phi_i\)，\(\phi_i\)表示\(x_i\)和\({F}_{M}^Hc\)间的角度，基于此可以得到以下结论：\(\begin{array}{c}|x_i|-\mathrm{Re}\left[\left(\mathcal{F}_M^H\mathbf{c}\right)_i^Hx_i\right]=|x_i|\left[1-|\left(\mathcal{F}_M^H\mathbf{c}\right)_i|\cos\phi_i\right]\geq|x_i|\left[1-|\left(\mathcal{F}_M^H\mathbf{c}\right)_i|\right].\end{array}\)

​ 当\(|\mathcal{F}_{M}^Hc|\leq1\)时下确界项为0；当\(|x_i|\left[1-|\left(\mathcal{F}_M^H\mathbf{c}\right)_i|\right]<0\)时下确界可以达到\(-\infty\)。

​ 3、整理上述讨论，有噪声下的原子范数的对偶问题可以表征为:

​ \(g(\mathbf{c})=\left\{\begin{array}{lcl}\operatorname{Re}\left[\mathbf{c}^H\mathbf{y}-\epsilon\Vert c\Vert_2\right],&\|\mathcal{F}_M^H\mathbf{c}\|_{\infty}\leq1\\ -\infty,&\quad\mathrm{otherwise.}\end{array}\right.\)

​ 在上式中，\(\mathcal{F}_M^H\mathbf{c}\)中\(\mathcal{F}_M^H\)表示逆FFT算子，对偶多项式可以表示为\(H(z)=\mathcal{F}_M^H\mathbf{c}=\sum\limits_{m=0}^{M-1}c_m z^m=\sum\limits_{m=0}^{M-1}c_m e^{-j\left(2\pi\frac{d}{\lambda}t\right)m}\)，其中\(z(t)=e^{-j\left(2\pi\frac{d}{\lambda}t\right)}\).

​ 4、为了进一步抽象\(g(\mathbf{c})\)，我们可以作以下表示：

​ 令\(a(\omega)=[1,e^{j\omega},...,e^{j\omega(L-1)}]^T\)为\(L-1\)次的三角多项式向量，那么因果三角多项式可以表征为：\(H(\omega)=\sum_{l=0}^{L-1}h_l e^{-j\omega l}=\mathbf{a}(\omega)^H\mathbf{h}\),其中\(\textbf{h}=\left[h_0,\cdots,h_{L-1}\right]^T\in\mathbb{C}^L\)表示多项式系数向量.

​ 对于非负三角多项式，可以有Hermitian矩阵\(R(\omega)=|H(\omega)|^2=H(\omega)H(\omega)^H=a(\omega)^Hhh^Ha(\omega)=\sum_{k=-(L-1)}^{L-1}r_k e^{-j\omega k}\),其中\(r_k=\sum_{l=0}^{L-1-k}h_l h_{l+k}^*\)，\(k\geq0\)并且\(r_{-k}=r_k^*\)，稀疏\(r_k\)可以通过自相关矩阵\(Q_{L\times L}=hh^H\)的第\(k\)条对角线元素进行计算\(r_k=\sum_{i=1}^{L-k}Q_{i,i+k}\).

​ 令两个多项式\(H(\omega)\)和\(B(\omega)\)满足以下不等关系：\(|H(\omega)|\leq|B(\omega)|,\forall\omega\in[-\pi,\pi]\)

​ 这意味着\(|H(\omega)|^2\leq|B(\omega)|^2,\forall\omega\in[-\pi,\pi]\),定义\(R_H(\omega)=|H(\omega)|^2\)和\(R_B(\omega)=|B(\omega)|^2\)，那么有\(R_H(\omega)\leq R_B(\omega)\)，即\(Q_H \leq Q_B\)，其中\(Q_H=hh^H\)和\(Q_B=bb^H\)为自相关向量\(\mathbf{h}=[h_{0},\cdots,h_{L-1}]^{T}\)和\(\mathbf{b}=[b_{0},\cdots,b_{L-1}]^{T}\)的自相关矩阵.

​ 根据Schur补条件有\(Q_B-hh^H\geq0\)，即\(\begin{bmatrix}\mathbf{Q}_B&\mathbf{h}_{L\times1}\\ \mathbf{h}_{1\times L}^{H}&1\end{bmatrix}\succeq0.\)

​ 令多项式\(H(\omega)\)的振幅均匀有界(对所有\(\omega\in[-\pi,\pi]\)有\(H(\omega)\leq\gamma\),其中\(\gamma \in R_+\)为给定正实数.作为有界三角多项式的特例，令\(|B(\omega)|=\gamma\)，那么\(H(\omega)\leq\gamma\)可以用两个线性不等式抽象，如下：(其中\(R_B(\omega)=\gamma^2\))

​ \(\begin{bmatrix}\mathbf{Q}_{L\times L}&\mathbf{h}_{L\times1}\\ \mathbf{h}_{1\times L}&1\end{bmatrix}\succeq0,\\ \sum_{i=1}^{L-j}Q_{i,i+j}=\left\{\begin{array}{c}\gamma^2,&j=0\\ 0,&j=1,\cdots,L-1.\end{array}\right.\)

​ 有界三角多项式的结果可以用于\(\infty\)范数，因为多项式的最大振幅设置上界意味着多项式对所有\(\omega\in[-\pi,\pi]\)具有一致有界的振幅,即\(\begin{array}{l}\|H\|_{\infty}=\max\limits_{\omega\in[-\pi,\pi]}|H(\omega)|\leq\gamma,\\ |H(\omega)|\leq\gamma,\forall\omega\in[-\pi,\pi].\end{array}\)

​ 回到本节开始处，基于振幅一致有界条件和Schur补条件,对偶问题可以表征为以下凸优化问题：

​ \(\begin{array}{c}\max\operatorname{Re}\left(\mathbf{c}^{H}\mathbf{y}-\epsilon\Vert c\Vert_2\right)~~\operatorname{subject to}\left[\begin{array}{cc}\mathbf{Q}_{M\times M}&\mathbf{c}_{M\times1}\\ \mathbf{c}_{1\times M}&1\end{array}\right]\succeq0,\\ \sum_{i=1}^{M-j}\mathbf{Q}_{i,i+j}=\left\{\begin{array}{cc}1,&j=0\\ 0,&j=1,\cdots,M-1.\end{array}\right.\end{array}\)

代码

% 本处仅给出上述凸优化问题的核心代码
if noise_flag == 0 % 无噪声版本
    cvx_begin sdp quiet
    cvx_solver sdpt3
        variable S(M+1,M+1) hermitian 
        subject to
        S >= 0;
        S(M+1,M+1) == 1;
        trace(S) == 2; % 主对角元素迹为2
        for j = 1 : M-1
            sum(diag(S,j)) == S(M+1-j,M+1); % 非主对角线元素求和为0.
        end
        maximize (real(S(1:M,M+1)'* Y)) %  - 0.5 * norm(c)
    cvx_end
else % noise version
    regular_param = 0.2; % 有噪声需要引入正则化参数
    cvx_begin sdp quiet
    cvx_solver sdpt3
        variable S(M+1,M+1) hermitian
        subject to
        S >= 0;
        S(M+1,M+1) == 1;
        trace(S) == 2;
        for j = 1 : M - 1
            sum(diag(S,j)) == S(M+1-j,M+1);
        end
         maximize (real(S(1:M,M+1)'* Y) - regular_param * norm(c));
    cvx_end
end

原子范数软阈值问题的推导

推导

​ 原子集合由各个正弦曲线的样本组成，\(a_{f,\phi}\in C^n\)，表示为\(a_{f,\phi}=e^{i2\pi\phi}\left[1~e^{i2\pi f}~\cdots e^{i2\pi(n-1)f}\right]^T\)

​ 无限原子集\(\mathcal{A}=\{a_{f,\phi}:f\in[0,1],\phi\in[0,1]\}\)组成了\(x^*\)适当的原子集合，\(x^*\)在对偶问题中可以写成一个稀疏的非负的原子组合。\(x^\star=\sum_{l=1}^k c_l^\star a_{f_l^\star,0}=\sum_{l=1}^k|c_l^\star|a_{f_l^\star,\phi_l}\)，\(c_{l}^{\star}=|c_{l}^{\star}|e^{i2\pi\phi_{l}}\).

​ 相应的对偶范数采用直观的形式：\(\|v\|_\mathcal{A}^*=\sup\limits_{f,\phi}\langle v,a_{f,\phi}\rangle=\sup\limits_{f\in[0,1]}\sup\limits_{\phi\in[0,1]}e^{i2\pi\phi}\sum\limits_{l=0}^{n-1}v_l e^{-2\pi i l f}=\sup\limits_{|z|\le1}\left|\sum\limits_{l=0}^{n-1}v_l z^l\right|\)，\(\|v\|_\mathcal{A}^*\)可以理解为在单位圆上获得的最大绝对值\(\zeta\mapsto\sum_{l=0}^{n-1}v_l{\zeta^l}\),\({\cal A}=\{a_{f,\phi}|f\in[0,1],\phi\in[0,1]\}\)为与线谱原子集相关的原子范数的半正定规划.

​ 根据上式可知向量\(v\in C^n\)的对偶原子范数是复数三角多项式\(V(f)=\sum_{l=0}^{n-1}v_l e^{j2\pi lf}\)的最大绝对值；因此，对对偶原子范数的约束等价于对\(V(f)\)大小的限制：\(||v||_A^\text{t}\leq\tau\Leftrightarrow|V(f)|^2\leq\tau^2,\forall f\in[0,1]\).函数\(q(f)=\tau^{2}-|V(f)|^{2}\)是一个三角多项式，\(q(f)\)非负的充要条件是可以写成三角多项式的平方和.

​ 定义映射\(T：\mathbb{C}^{n}\to\mathbb{C}^{n\times n}\),从输入创建一个Hermitian Toeplitz 矩阵，即\(T(x)=\begin{bmatrix}x_1&&x_2&&...&&x_n\\ x_2^*&&x_1&&...&&x_{n-1}\\ \vdots&&\vdots&&\ddots&&\vdots\\ x_n^*&&x_{n-1}^*&&...&&x_1\end{bmatrix}\).

​ 对于给定的因果三角多项式\(V(f)=\sum_{l=0}^{n-1}v_{l}e^{-2\pi i l f},\)如果有且仅有复Hermitian矩阵\(Q\)存在时有\(|V(f)|\leq\tau\)，这与原子范数对偶问题中第4节证明类似，即有\(T^*(Q)=\tau^2e_1~\mathrm{and}~\begin{bmatrix}Q&v\\ v^*&1\end{bmatrix}\succeq0.\)其中\(e_1=[1,0,0,....,0]^T\)，\(v^*\)表示\(v\)的Hermitian转置.

​ 重写原子范数\(\Vert x \Vert_A=\sup_{\|v\|_\mathcal{A}^*\leq1}<x,v>\)为下列形式：

​ \(\begin{array}{ll}\text{maximize}_{v,Q}&\langle x,v\rangle\\ \text{subject to}&T^*(Q)=e_1\\ &\begin{bmatrix}Q&v\\ v^*&1\end{bmatrix}\succeq0.\end{array}\)

​ 下面对上述问题进行对偶推导：

​ 1、首先需要将上述问题转化为无约束的拉格朗日方程形式，可以表示如下：

​ $L(Q,v,u,\Gamma)=\langle x,v\rangle +\langle u, T^* (Q)\rangle-\langle u,e_1\rangle-\langle\Gamma, \begin{bmatrix}Q&v\ v^*&1\end{bmatrix}\rangle $

​ 2、关于\(v\)的项为$\langle x,v\rangle-\langle\Gamma, \begin{bmatrix}Q&v\ v^*&1\end{bmatrix}\rangle \(，下面对变量\)v$求解极值，则有

​ \(x- Tr(\Gamma\begin{bmatrix}0&I\\ I^*&0\end{bmatrix})=x- Tr(\begin{bmatrix}\Gamma_{12}I^*&\Gamma_{11}I\\ \Gamma_{22}I^*&\Gamma_{21}I\end{bmatrix})=x-2\Gamma_{12}=0\)，那么可以得到\(\Gamma_{12}=\frac{x}{2};\)\(\Gamma_{21}=\frac{x^*}{2}\)

​ 3、关于\(Q\)的项为\(\langle u, T^* (Q)\rangle-\langle\Gamma, \begin{bmatrix}Q&v\\ v^*&1\end{bmatrix}\rangle\)，对变量\(Q\)求解极值前，先将\(\langle u, T^* (Q)\rangle\)进步抽象为\(Tr(T(u)\cdot Q)\),那么关于\(Q\)的偏导可表示为\(T(u)-\langle\Gamma,\begin{bmatrix}I&0\\ 0&0\end{bmatrix}\rangle=0\)，那么则有\(\Gamma_{11}=T(u)\),其中\(F_{22}=t\)，用于半正定约束\(\Gamma=\begin{bmatrix}T(u)&x/2\\ x^*/2&t\end{bmatrix}\geq0\).

​ 4、将\(\Gamma\)结果代入到\(L\)中，那么有如下证明：

​ \(L=Re(v^*x)+Tr(T(u)\cdot Q)-u^*e_1-Tr(\begin{bmatrix}T(u)Q+Re(v^*x)/2&x/2+T(u)v\\ tv^*+Re(Qx^*)/2&Re(x^*v)/2+t\end{bmatrix} ) =-u^*e_1-t=-u_1-t\).

​ 根据半正定约束条件\(T(u)t-xx^*/4\geq0\)，通过对\(u\)和\(t\)缩放则有\(2t\cdot T(2u)-xx^*\geq0\)

​ 这等价于将对应目标函数缩放为\(-u_1/2-t/2\)，那么原问题的对偶形式可以表示如下:

​ \(\begin{array}{l l}{{\mathrm{minimize}_{t,u}}}&{{\frac{1}{2}(t+u_{1})}}\\ {{\mathrm{subject~to}}}&{{\begin{bmatrix}{T(u)}&{x}\\ {x^{*}}&{t}\end{bmatrix}\succeq0.}}\end{array}\)

​ 那么对应有噪声版本下的原问题对偶函数可以表示如下：[\(\tau\)表示正则参数]

​ \(\begin{array}{ll}\text{minimize}_{t,u,x}&\frac{1}{2}\|x-y\|_2^2+\frac{\tau}{2}(t+u_1)\\ \text{subject to}&\begin{bmatrix}T(u)&x\\ x^*&t\end{bmatrix}\succeq0.\end{array}\)

​ 上述问题可以通过凸优化中的SDP解释器求解，但是计算复杂度较高，可以通过交替方向投影算子加速求解,这将在后续的章节进一步讨论。

代码

% 在上述推导过程中讨论了单快拍下有噪声和无噪声版本的原子范数模型
% 在本代码中笔者给出了单快拍和多快拍版本，后续将补充多快拍版本的理论
if noise_flag == 0 % 无噪声情况下的原子范数AST模型
    if snap == 1 % 单快拍模型
        cvx_begin sdp quiet
        cvx_solver sdpt3
            variable T(M, M) hermitian toeplitz
            variable x
            minimize (0.5 * x + 0.5 * T(1,1))
            [x Y'; Y T] >= 0;
        cvx_end
        [Phi, Val] = rootmusic(T, P, 'corr');
        Phis = Phi / 2 / pi ;
        estimated_theta = asind(-Phis * lambda / d);
    else % 多快拍模型
        cvx_begin sdp quiet
        cvx_solver sdpt3
            variable T(M, M) hermitian toeplitz
            variable X(snap, snap) hermitian
            minimize (trace(X)+trace(T))
            [X Y'; Y T] >= 0;
        cvx_end   
        [Phi, Val] = rootmusic(T, P, 'corr');
        Phis = Phi / 2 / pi ;
        estimated_theta = asind(-Phis * lambda / d);
    end
    
else % 有噪声情况下的原子范数AST模型
    if snap == 1 % 单快拍模型
        sigma = 1;
        regular_param = sqrt(M * log(M * sigma));
        cvx_begin sdp quiet
        cvx_solver sdpt3
            variable T(M, M) hermitian toeplitz
            variable x
            variable z(M,1) complex
            minimize (regular_param * 0.5 *(x + T(1,1)) + 0.5 * norm(Y-z))
            [x Y'; Y T] >= 0;
        cvx_end
        [Phi, Val] = rootmusic(T, P, 'corr');
        Phis = Phi / 2 / pi ;
        estimated_theta = asind(-Phis * lambda / d);
        
    else % 多快拍模型
        regular_param = sqrt(M * (snap + log(M) + sqrt(2 * snap * log(M))));
        cvx_begin sdp quiet
        cvx_solver sdpt3
            variable T(M,M) hermitian toeplitz
            variable X(snap, snap) hermitian
            variable Z(M, snap) complex
            minimize (regular_param * (trace(X) + trace(T)) + 1 / 2 * sum_square_abs(vec(Y - Z)));
            [X Y';Y T] >= 0;
        cvx_end
        [Phi, Val] = rootmusic(T, P, 'corr');
        Phis = Phi / 2 / pi ;
        estimated_theta = asind(-Phis * lambda / d);
    end
end

参考文献

[1] Atomic norm denoising with applications to line spectral estimation. https://arxiv.org/abs/1204.0562

[2] Grid-free compressive beamforming. https://arxiv.org/abs/1504.01662

[3] Positive Trigonometric Polynomials and Signal Processing Applications.

[4] Regularized Matrix Factorization for Multilabel Learning With Missing Labels. https://ieeexplore.ieee.org/abstract/document/9198894