1比特压缩感知


1比特压缩感知(One-bit compressive sensing,1BCS)旨在从以下系统中恢复出稀疏信号 $\mathbf{x}\in\mathbb{R}^{n}$,

\begin{equation} \mathbf{b} = \mathrm{Diag}(\mathbf{h}) \mathrm{sign}(\mathbf{A}\mathbf{x} + \boldsymbol{\varepsilon}) \tag{1bCS} \end{equation}

其中,采样矩阵 $\mathbf{A}\in\mathbb{R}^{m\times n}$,观测向量 $\mathbf{b}$ 和符号翻转向量 $\mathbf{h}\in\{-1,1\}^{m}$,噪声 $\boldsymbol{\varepsilon}\in\mathbb{R}^{n}$,对角矩阵 $\mathrm{Diag}(\mathbf{h})$ 的对角元由 $\mathbf{h}$ 组成,符号函数定义为 $\mathrm{sign}(t)=1$ 当 $t>0$,否则 $\mathrm{sign}(t)=-1$。向量情形下,$\mathrm{sign}(\mathbf{x})=$ $(\mathrm{sign}(x_1)$ $\cdots$ $\mathrm{sign}(x_n))^\top$。值得一提的是,观测 $\mathbf{b}$ 是符号翻转后得到的,因而,问题更具挑战性。在该模型中,假设最多有 $k$ 个符号被翻转,即 $\mathbf{h}$ 满足 $\parallel\mathbf{h}-1\parallel_0\leq k$,其中 $k$ 为给定整数,零范数 $\parallel\mathbf{x}\parallel_0$ 表示 $\mathbf{x}$ 中非零元个数。为恢复信号,可求解以下优化模型。

◻️ 双稀疏约束优化模型 \begin{equation} \min_{\mathbf{x}\in\mathbb{R}^{n},\mathbf{z}\in\mathbb{R}^{m}}~ \parallel \mathrm{Diag}(\mathbf{b}) \mathbf{A} \mathbf{x}+\mathbf{z} -\epsilon \parallel^2 + \eta \parallel \mathbf{x} \parallel^2,~~~\textrm{s.t.}~ \parallel\mathbf{x} \parallel_0\leq s,~ \parallel \mathbf{z}_+\parallel_0\leq k \tag{DSCO} \end{equation}
◻️ 阶梯函数正则优化模型 \begin{equation} \min_{\mathbf{x}\in\mathbb{R}^{n}}~ \sum_{i=1}^n (x_i^2+\varepsilon)^{q/2}+\lambda \parallel (\epsilon- \mathrm{Diag}(\mathbf{b}) \mathbf{A} \mathbf{x})_+ \parallel_0 \tag{SFRO} \end{equation}
其中, 正整数 $s\ll n$ 和 $k\ll m$, 参数 $(\epsilon,\eta,\varepsilon,\lambda)>0$ 且 $\mathbf{z}_+=$ $(\max\{0,z_1\}$ $\cdots$ $\max\{0,z_m\})^\top$。度量 $\parallel \mathbf{z}_+\parallel_0$ 计算 $\mathbf{z}$ 中正元素的个数,它与阶梯(又称 0/1 损失)函数相关,其定义为 $\mathrm{step}(t)=1$ 当 $t>0$,否则 $\mathrm{step}(t)=0$。因此,$\|\mathbf{z}_+\|_0=$ $\mathrm{step}(z_1)+$ $\cdots$ $+\mathrm{step}(z_m)$。

程序包 - OBCSpack-Matlab(点击下载)提供了 2 个求解器。其核心算法分别来自以下 2 篇文章,其中, $\texttt{GPSP}$ 和 $\texttt{NM01}$ 分别用来求解模型 (DSCO) 和 (SFRO)。

GPSP - S Zhou, Z Luo, N Xiu, and G Li, Computing one-bit compressive sensing via double-sparsity constrained optimization, IEEE TSP, 70:1593-1608, 2022.
NM01 - S Zhou, L Pan, N Xiu, and H Qi, Quadratic convergence of smoothing Newton’s method for 0/1 loss optimization, SIOPT, 31:3184–3211, 2021.


下面 Matlab 代码(Python 代码类似)展示了如何使用 $\texttt{OBCSpack}$ 求解 1BCS 问题。用户需输入数据 ($\texttt{A}$,$\texttt{b}$,$\texttt{s}$,$\texttt{k}$),然后从 {'$\texttt{GPSP}$','$\texttt{NM01}$'} 中选择一个求解器进行求解。

clc; close all; clear; addpath(genpath(pwd));

n      = 2000;                             % Signal dimension 
m      = ceil(0.5*n);                      % Number of measurements
s      = ceil(0.01*n);                     % Sparsity level
nf     = 0.05;                             % Noisy ratio
r      = 0.02;                             % Flipping ratio
k      = ceil(r*m);

A      = randn(m,n);
T      = randperm(n,s);
xo     = zeros(n,1);                      
xo(T)  = (1+rand(s,1)).*sign(randn(s,1));  
xo(T)  = xo(T)/norm(xo(T));                % True sparse solution
h      = ones(m,1);                        % Flipping vector
T      = randperm(m,k); 
h(T)   = -h(T);
b      = h.*sign(A(:,T)*xo(T)+nf*randn(m,1));; 

solver = {'GPSP','NM01'};
out    = OBCSpack(A,b,s,k,solver{1});  
fprintf(' Time:                  %6.3f sec\n',out.time);
fprintf(' Absolue error:         %6.2f %%\n', norm(xo-out.sol)*100);
fprintf(' Signal-to-noise ratio: %6.2f\n',-10*log10(norm(xo-out.sol)^2));
fprintf(' Hamming distence:      %6.3f\n',nnz(sign(A*out.sol)-b)/m)
Matlab 版程序包 $\texttt{OBCSpack}$ 的输入和输出(Python 版的输入和输出类似)如下所示,其中输入 ($\texttt{A}$,$\texttt{b}$,$\texttt{s}$,$\texttt{k}$,$\texttt{solver}$) 为必需项。因为 $\texttt{NM01}$ 求解模型 (SFRO),所以无需参数 $\texttt{s}$ 和 $\texttt{k}$。因此,如果选择 $\texttt{solver}$='$\texttt{NM01}$',则当 $\texttt{s}$ 和 $\texttt{k}$ 未知时,可以将它们设置为 $\texttt{[]}$。参数 $\texttt{pars}$ 是可选的,但设置其中的一些参数可以提升求解器的性能和解的质量。

function out = OBCSpack(A,b,s,k,solver,pars)
% -------------------------------------------------------------------------
% One-bit compressed sensing problem aims to recover a sparse signal x from
%
%                b = Diag(h).*sign( A*x + noise )
%
% 1) The double sparsity constrained optimization (DSCO)
%
%    min  ||Diag(b)*A*x+y-epsilon||^2 + eta||x||^2
%    s.t. ||x||_0<=s, ||y_+||_0<=k
%
% where (epsilon, eta)>0, s\in[1,n], k\in[0,m] are given.
%
% 2) The step function regularized optimization (SFRO)
%
%    min ||x.*x+vareps||^{q/2}_{q/2} + lam*||(epsilon - Diag(b)*A*x)_+||_0
%
% where (vareps, lam, epsilon)> 0, q\in(0,1).  
% -------------------------------------------------------------------------
% Inputs:
%  A:       The sensing matrix \in R^{m-by-n},                   (REQUIRED)
%  b:       The binary observation \in R^m, b_i\in{-1,1}         (REQUIRED)
%  s:       Sparsity level of x, an integer \in[1,n]             (REQUIRED)      
%  k:       An integer in [0,m], e.g., k = ceil(0.01m)           (REQUIRED)       
%  solver:  A text string, can be one of {'GPSP','NM01'}         (REQUIRED)            
%  pars:    Parameters are optional                              (OPTIONAL) 
%           ------------- For GPSP solving (DSCO)--------------------------
%           pars.eps     - The parameter in the model        (default 1e-4)
%           pars.eta     - The penalty parameter       (default 0.01/ln(n))
%           pars.acc     - Acceleration is used if acc=1        (default 0)
%           pars.big     - Start with a bigger s if big=1       (default 1)
%           pars.maxit   - Maximum number of iterations       (default 1e3) 
%           pars.tol     - Tolerance of halting condition    (default 1e-8)
%           -------------  For NM01 solving (SFRO)-------------------------
%           pars.x0      - The initial point           (default zeros(n,1))
%           pars.q       - Parameter in the objective         (default 0.5)
%           pars.vareps  - Parameter in the objective         (default 0.5)
%           pars.epsilon - Parameter in the objective        (default 0.15)
%           pars.lam     - The penalty parameter                (default 1)
%           pars.tau     - A useful parameter                   (default 1) 
%           pars.maxit   - Maximum number of iterations       (default 1e3)  
% -------------------------------------------------------------------------
% Outputs:
%     out.sol:   The sparse solution x
%     out.time:  CPU time
%     out.iter:  Number of iterations
%     out.obj:   Objective function value at out.sol 
% ------------------------------------------------------------------------