基于KL不等式求解压缩感知问题
Solving Compressed Sensing Problem Based on Kurdyka-?ojasiewicz Inequality
DOI: 10.12677/AAM.2022.111054, PDF, HTML, XML, 下载: 400  浏览: 669  科研立项经费支持
作者: 孙静坤, 许 荻, 任咏红*:辽宁师范大学,辽宁 大连
关键词: KL性质压缩感知线搜索KL Properties Compressed Sensing Line Search
摘要: 压缩感知问题在雷达探测、信号与图象处理、组合优化、金融等诸多领域中有着非常广泛的应用。本文将压缩感知问题在一定的条件下转化为等价的无约束问题。通过研究这类无约束问题的KL (Kurdyka-Łojasiewicz)性质,本文采用BB步长,应用相应的非单调线搜索的方法求解这类非线性函数的收敛性和线性收敛速率。
Abstract: Compression sensing problem is widely used in radar detection, signal and image processing, combination optimization, finance and many other fields. In this paper, the compressed sensing problem is transformed into an equivalent unconstrained problem under certain conditions. By studying the Kurdyka-Łojasiewicz (KL) properties of such unconstrained problems, we use BB steps to solve the linear convergence and linear convergence rate of such nonlinear functions.
文章引用:孙静坤, 许荻, 任咏红. 基于KL不等式求解压缩感知问题[J]. 应用数学进展, 2022, 11(1): 462-472. https://doi.org/10.12677/AAM.2022.111054

1. 引言

Donho、Candes和Tao等人 [1] 在研究高度欠定的核磁共振成像问题时,正式提出了压缩感知的概念。压缩感知,也被称为压缩传感,它的基本思想是利用特定的矩阵,将K-稀疏的高维信号空间变换为低维信号空间,利用稀疏性信号的先验条件,通过某个线性或者非线性重建模型重建原始信号。压缩感知理论广泛应用于雷达探测、图像处理、无线传感网络等很多领域。压缩感知问题的求解备受关注。

压缩感知问题的模型如下:

min x 0 s .t . y = A x (1)

其中,A代表 m × n ( m n )的观测矩阵,x表示原始信号,且 x n ,y代表观测值,且 y m x 0 表示向量x的零范数(即非零元的个数)。求解问题(1),需要列出x中所有非零位置的线性组合,有 C n k 种。问题(1)的数值计算很不稳定,并且是一种NP (Nondeterministic Polynomially)难的非凸优化问题。具有代表性的求解方法有匹配追踪(MP) [2],正交匹配追踪(OMP) [3],稀疏自适应匹配追踪(SAMP) [4] 等。Donohoh、Candes和Tao等人在文献 [5] 中将其转化成与之等价的 l 1 范数进行求解,证明了若观测矩阵A满足RIP (Restricted Isometry Property),那么问题(1)与 l 1 范数问题等价,并且建立了问题(1)的等价形式如下 [6]:

min x 1 s .t . y = A x (2)

本文用罚函数法,将问题(2)转化为与其同解的无约束问题:

min x R n { f ( x ) : = λ 2 A x y 2 + x 1 }

其中, λ > 1

KL不等式作为收敛性分析中的一个重要研究工具,并且应用广泛。像半代数函数、tame函数、向量范数等等都是KL函数。早在1963年S. Łojasiewicz [7] 在实分析函数中,研究一般的最速下降曲线问题解的轨迹是否为有限长度时,提出Łojasiewicz不等式:

F ( x ) c | F ( x ) | θ x U ( x ¯ )

其中, F : n 的实解析函数, x ¯ F 1 ( 0 ) 为临界点且 θ [ 1 2 , 1 ) c > 0 。并由此推导出最速下降曲线

的解的轨迹是有限长度,有界轨迹收敛到临界点。随后,Kurdyka [8] 将这一结果扩展到可在o-minimal结构中定义的可微函数中,相应的广义不等式称之为KL不等式。本文主要建立与其同解的无约束优化问题模型的KL性质,从而设计基于BB (Barzilai-Borwein)步长线搜索算法。

2. 预备知识

对于给定的广义实值函数 f : X ( , + ] ,若

dom f : = { x X | f ( x ) < }

则称 f ( x ) 是正常的。对给定的 α β ,记

[ α f β ] : = { x X | α f ( x ) β }

对于给定的 x ¯ X δ > 0 B ( x ¯ , δ ) 表示以 x ¯ 为中心, δ 为半径的闭球。x到S的距离定义为

dist ( x , S ) : = inf { y x : y S }

定义1 [9] (Kurdyka-Łojasiewicz性质)设 f : X ( , + ] 是正常的下半连续的函数。考虑任意 x ¯ dom f 。若存在常数 η ( 0 , + ] ,满足如下条件的连续凹函数 φ : [ 0 , η ) +

1) φ ( 0 ) = 0 φ 在开区间 ( 0 , η ) 上连续可微;

2) 对所有的 s ( 0 , η ) φ ( s ) > 0 ,以及点 x ¯ 的邻域U使得对所有 x U [ f ( x ¯ ) < f < f ( x ¯ ) + η ] ,都有下述不等式成立

φ ( f ( x ) f ( x ¯ ) ) dist ( 0 , f ( x ) ) 1

则称函数f在点 x ¯ 处具有KL性质。若f在集合 dom f 中的每点都具有KL性质,则称f是KL函数。

定义2 [10] (Lipschitz连续)设F为从集合 D n m 上的单值映射。设 X D 。F在X上是Lipschitz连续的,若存在 κ + = [ 0 , + ) ,满足

F ( x ) F ( x ) κ x x , x , x X

其中 κ 称为F在X上的Lipschitz常数。

定义3 [10] (次微分)函数 f ( x ) 在x处的所有次梯度的集合称为 f ( x ) 在x处的次微分,即

f ( x ) = { x X : f ( y ) f ( x ) x , y x , y X }

dom F ( x ) = { x | F ( x ) }

定义4 [11] 设存在 f 0 R ,S是一非空子集,满足 f ( x ) = f 0 x S ,令 γ 是一正的常数。称 γ 阶增长条件在S上成立,若存在常数 c > 0 ,存在S的邻域V,满足对所有的 x n V ,下述不等式成立:

f ( x ) f 0 + c [ dist ( x , S ) ] γ

尤其,若 γ = 1 ,称一阶增长条件成立,若 γ = 2 ,称第二阶(或二阶)增长条件成立。

定义5 [11] (凸函数)设C是一凸集合, f : C ( , + ] 是一增广实值函数。对任意 x C y C ,有

f ( ( 1 t ) x + t y ) ( 1 t ) f ( x ) + t f ( y ) 0 < t < 1

则函数f是C上的凸函数。

引理1 对任意给定的 x n

1 ( x ) = { μ R n | μ i = { 1 x i > 0 1 x i < 0 [ 1 , 1 ] x i = 0 , i = 1 , 2 , , n }

注记1 根据文献 [12] 中的引理2.1,一个正常函数在所有非稳定点处具有KL性质。因此为了证明它是KL函数,只需证明函数在任意稳定点处是否具有KL性质。

注记2 若 0 f ( x ¯ ) ,则称 x ¯ 是函数f的稳定点,记 crit f 为f的稳定点集。

3. 主要结果及其应用

考虑如下问题模型

min x R n { f ( x ) : = λ 2 A x b 2 + x 1 } (3)

其中, A m × n ( m < n ), x n b m λ > 1 为参数, x = x 1 2 + x 2 2 + + x n 2

x 1 = | x 1 | + | x 2 | + + | x n |

f 1 ( x ) : = λ 2 A x b 2 f 2 ( x ) : = x 1 ,则 f ( x ) = f 1 ( x ) + f 2 ( x )

引理2 [13] 对于任意的 x n ,给定 p > 0 ,则p范数,即 x p = ( i = 1 n x i p ) 1 p ,具有KL性质。

由引理2知, f 2 ( x ) 具有KL性质。

引理3 [14] 假设 f 1 dom f 2 上具有 L > 0 的Lipschitz连续梯度,即

f 1 ( x ) f 1 ( y ) L x y x , y dom f 2

f ( x ) 是下有界的。则

f ( x ) = { f 1 ( x ) } + f 2 ( x ) x dom f

定理1 假设 f 1 ( x ) dom f 上局部Lipschitz连续,对 x ¯ dom f ,若 f 1 ( x ) x ¯ 点处满足二阶增长条件,即

f 1 ( x ) f 1 ( x ¯ ) + c x x ¯ 2 c > 0

f ( x ) 为KL函数。

证明:由 f 1 ( x ) = λ 2 A x b 2 x n ,可知 f 1 ( x ) = λ * A T ( A x b ) 2 f 1 ( x ) = A T A 。因此 f 1 ( x ) 是凸的,且在 x ¯ 点满足二阶增长条件。由次微分定义,有下式成立

f 1 ( x ¯ ) f 1 ( x ) λ A T ( A x b ) , x ¯ x x n

f 1 ( x ) f 1 ( x ¯ ) λ A T ( A x b ) , x x ¯ x n

因此有

f 1 ( x ) f 1 ( x ¯ ) λ A T ( A x b ) x x ¯ x n (4)

由于 f 1 ( x ) x ¯ 点处二阶增长条件,可得

f 1 ( x ) f 1 ( x ¯ ) c x x ¯ 2 x n (5)

将(4)和(5)结合,可得

c x x ¯ 2 λ A T ( A x b ) x x ¯ x n

c λ A T ( A x b ) x x ¯ x n (6)

考虑 φ ( s ) = 2 c s s ( 0 , η 1 ) ,其中 η 1 ( 0 , + ] ,有

φ ( s ) = 1 c s φ ( s ) = 1 2 c s 3

φ ( s ) [ 0 , η 1 ) 上是凹函数。所以对于 δ 1 > 0 x B ( x ¯ , δ 1 ) [ f 1 ( x ¯ ) < f 1 ( x ) < f 1 ( x ¯ ) + η 1 ] ,有

φ ( f 1 ( x ) f 1 ( x ¯ ) ) = 1 c ( f 1 ( x ) f 1 ( x ¯ ) )

将(4)式与上式结合,可得到下式

φ ( f 1 ( x ) f 1 ( x ¯ ) ) 1 c λ A T ( A x b ) x x ¯ (7)

由距离函数的定义,有

dist ( 0 , f 1 ( x ) ) = λ A T ( A x b ) (8)

再将(7)式与(8)式结合,得

φ ( f 1 ( x ) f 1 ( x ¯ ) ) dist ( 0 , f 1 ( x ) ) 1 c λ A T ( A x b ) x x ¯ λ A T ( A x b ) (9)

将(6)式带入(9)式,就可以得到

φ ( f 1 ( x ) f 1 ( x ¯ ) ) dist ( 0 , f 1 ( x ) ) 1

f 1 ( x ) x ¯ 点处满足KL性质。由 x ¯ 的任意性,则 f 1 ( x ) 在集合 dom f 中的每点都具有KL性质,即 f 1 ( x ) 是KL函数。

由于函数 f 1 ( x ) x ¯ 处具有连续性,因此对 η 2 ( 0 , 1 ) δ 2 > 0 有下式成立

| f 1 ( x ) f 1 ( x ¯ ) | η 2 x B ( x ¯ , δ 2 )

f 1 ( x ) 是KL函数。 结合定义1,对 x B ( x ¯ , δ 1 ) [ f 1 ( x ¯ ) < f 1 ( x ) < f 1 ( x ¯ ) + η 1 ] 有下式成立:

dist ( 0 , f 1 ( x ) ) c f 1 ( x ) f 1 ( x ¯ ) (10)

δ = min { δ 1 , δ 2 } η = min { η 1 , η 2 } 。对 x B ( x ¯ , δ ) [ f ( x ¯ ) < f ( x ) < f ( x ¯ ) + η ] ,有

min z F ( x ) z = dist ( 0 , f ( x ) ) (11)

f ( x ) 是闭集,则 z f ( x ) ,满足 z = dist ( 0 , f ( x ) ) 。由于 f 1 x ¯ 处具有局部Lipschitz连续性,结合引理3,有 f ( x ) = f 1 ( x ) + 1 ( x ) 。结合引理1中的 1 ( x ) 的表达式以及式(11),有

z 2 = dist 2 ( 0 , f 1 ( x ) + 1 ( x ) ) λ 2 A T ( A x b ) 2 (12)

通过式(11)和式(12),有下式:

dist ( 0 , f ( x ) ) 2 = z 2 λ 2 A T ( A x b ) 2 = dist ( 0 , f 1 ( x ) ) 2 (13)

f 1 ( x ¯ ) < f 1 ( x ) f ( x ) < f ( x ¯ ) + η 结合,可以得到 x 1 x ¯ 1 。由式(10)和(13),有

dist ( 0 , f ( x ) ) 2 dist ( 0 , f 1 ( x ) ) 2 c ( f 1 ( x ) f 1 ( x ¯ ) ) c ( f ( x ) f ( x ¯ ) )

上述不等式说明了 f ( x ) x ¯ 处具有KL性质。由于 x ¯ dom f 中的任意性,故 f ( x ) 是KL函数。

4. 非单调线搜索收敛算法

本节将采用非单调线搜索算法求步长,采用BB [15] 步长作为每次迭代的初始步长,经过有限步的线搜索,就可以找到符合收敛条件的步长,该算法充分说明了该问题模型满足KL性质,基于KL性质可以分析算法的收敛性。

g k f ( x k ) ,Hessian矩阵 H k 2 f ( x k ) ,其中 k N ,则牛顿迭代为

x k + 1 = x k ( H k ) 1 g k (14)

用BB法可以得到Hessian矩阵 H k α k 1 I ,其中 α k > 0 ,I为 n × n 单位矩阵。将(14)式中的 H k 替换为 α k 1 I ,得到

x k + 1 = x k α k g k (15)

Δ g k = g k g k 1 Δ x k = x k x k 1 ,则

g k g k 1 = H k ( x k x k 1 )

Δ g k = H k Δ x k

将上式中的 H k 替换为 α k 1 I ,得到 Δ g k = α k 1 Δ x k ,所以 α k ( Δ x k ) T Δ g k = ( Δ x k ) T Δ x k ,即 α k = Δ x k 2 Δ x k , Δ g k 。令 B k = α k I ,迭代格式(15)也可以写为

x k + 1 = x k B k g k (16)

BB法选取的步长,满足下式:

α k = arg min α α Δ g k Δ x k 2 2 α k = arg min α Δ g k α 1 Δ x k 2 2

即步长为

α k B B 1 = Δ x k 2 Δ x k , Δ g k α k B B 2 = Δ x k , Δ g k Δ g k 2 (17)

α k 为第k次迭代的初始步长。设L是 A T A 的特征值中最大的一个数,初始步长 α 0 = 1 L ,为了保证全局收敛性,非单调步长准则,如下式

f ( x k α k g k ) C k c 1 α k g k 2 (18)

其中 C k 满足 C 0 = f ( x 0 ) C k + 1 = 1 Q k + 1 ( τ Q k C k + f ( x k + 1 ) ) Q 0 = 1 Q k + 1 = τ Q k + 1 c 1 , τ ( 0 , 1 )

算法的具体步骤如下:

输入: A , b , x 0 , τ ( 0 , 1 ) , c 1 ( 0 , 1 ) , ε ( 0 , 1 ) , Q 0 = 1 , λ = 100 , k = 0 , α 0

输出: x = x k + 1

步骤1:初始化 g 0 = λ A T ( A x 0 b ) + s i g n ( x 0 ) C 0 = f ( x 0 )

步骤2:当 g k > ε 时,向下进行;否则停止计算;

步骤3:当 f ( x k α k g k ) > C k c 1 α k g k 2 时, α k = 0.6 α k ;否则向下进行;

步骤4:迭代 x k + 1 = x k α k g k

步骤5:通过(17)式获取步长 α k + 1 ,且满足 α k + 1 [ 10 20 , 10 20 ] ,再计算 C k + 1 Q k + 1 ;令 k = k + 1 ,返回步骤2。

变量 C k 实际上为本次搜索准则的参照函数值,即充分下降性质的起始准则; C k + 1 是函数值 f ( x k + 1 ) C k 的凸组合,并不是仅仅依赖于 f ( x k + 1 ) ,而凸组合的两个系数由参数 τ 决定。这就说明 C k 包含之前算过的所有目标函数值。综上所述,本文提出的算法是可行有效的。

5. 实验结果分析及讨论

通过MATLAB实验来证明该算法在重构信号方面的优势,实验运行平台的配置如下:Intel(R)Core (TM)i5-10400F CPU @2.90 GHz;内存16.00 GB。在不同的MATLAB版本测试下,结果不会发生较大变化。本文实验在MATLAB R2020b中进行。

实验取 n = 2048 m = 512 a = 64 ,a表示原始信号中非零元的个数。其中非零元的位置均是随机选取的,A是与高斯矩阵独立同分布的随机矩阵 [16], b = A x 0 + ω ,其中 ω m ω 是分布为 N ( 0 , σ 2 I ) 的高斯噪声,取 σ 2 = 10 3 ε = 10 4 α k = min { max { α k , 10 20 } , 10 20 } 。在线搜索过程中,选取初始步长 α 0 = 0.8 c 1 = 10 6 τ = 0.85 。相对误差定义为

RelErr = x ^ x ¯ 2 x ¯ 2

其中 x ^ 表示重构信号, x ¯ 表示原始信号。相对误差用来测量重构信号的质量。当相邻两点的相对误差充分小,即

x k x k 1 2 x k 1 2 < ε

时,则停止迭代。原始信号,观测向量及重构信号如图1所示。

比较图1中的(a)和(c)可以分析出来,所有的蓝点均被红点包围,这说明原始信号几乎被精确重构。上述实验表明,该算法效果良好,相对误差较小,为恢复大的稀疏信号提供了一种有效的方法。

本文在接下来的实验中,采用了4种不同的信号,对不同的算法进行比较分析。取上述实验中的初始参数,对NBBL1n和NBBL1这2种算法进行了对比和分析。在初始和终止条件相同的情况下,从函数值下降的角度分析,无论针对那种信号,NBBL1n算法下降总是最快的,迭代次数也最少,为进一步说明结果,表1表2给出了详细的对比数据。由于设下随机种子,所以每次获取的信号是不变的。所以误差很小。

(a) 原始信号 (b) 观测向量(c) NBBL1n (RelErr = 1.38%)

Figure 1. The effectiveness of the NBBL1n algorithm

图1. NBBL1n算法的有效性

Table 1. Comparison of the number of iterations of the time required by the algorithm for different sparsity degrees

表1. 对于不同稀疏度,算法所需要的时间的迭代次数对比

表1中详细的列出了4种不同信号的运行时间和迭代次数。从表中可以分析出,随着信号规模的不断增大,运行算法所需要的时间也明显增加,但是通过对比,新版本算法的运行总用时较少,迭代次数也比较少。

Table 2. When c taken 0.4, each step iterative gradient norm vs c ( f ( x k ) f ( x ¯ ) )

表2. 取 c = 0.4 时,每步迭代梯度范数与 c ( f ( x k ) f ( x ¯ ) ) 的对比

表2取稀疏度为0.04,为了表示一般性,45次迭代做了对比分析,充分说明了新算法满足KL不等式。再对比表1,可以得到,满足KL性质的非单调线搜索算法,无论稀疏度如何,算法均可以达到最优。

6. 结束语

本文将结合非单调线搜索算法和BB步长,在满足KL性质的基础上,将算法应用于压缩感知问题模型,数值实验结果表明,当原始信号的稀疏度不同时,该方法是可行有效的。本文的方法无论是从运行时间还是迭代次数上均有不同程度的改进。但将本文的非单调线搜索应用到其他问题中是否有较好的效果,仍需要进一步研究。

基金项目

辽宁省教育厅科学研究一般项目(LJ2019005)。

NOTES

*通讯作者。

参考文献

[1] Donoho, D.L. (2006) Compressed Sensing. IEEE Transactions on Information Theory, 52, 1289-1306.
[2] Mallat, S.G. and Zhang, Z.F. (1993) Matching Pursuits with Time-Frequency Dictionaries. IEEE Transactions on Signal Processing, 41, 3397-3415.
[3] Tropp, J.A. and Gilbert, A.C. (2007) Signal Recovery from Random Measurements via Orthogonal Matching Pursuit. IEEE Transactions on Information Theory, 53, 4655-4666.
[4] Do, T.T., Gan, L., Nguyen, N., et al. (2008) Sparsity Adaptive Matching Pursuit Algorithm for Practical Compressed Sensing. Proceedings of the 42nd Asilomar Conference on Signals, Systems and Computers, Pacific Grove, 26-29 October 2008, 581-587.
[5] Chen, S.S., Donoho, D.L. and Saunders, M.A. (2001) Atomic Decomposition by Basis Pursuit. Society for Industrial and Applied Mathematics, 43, 29-159.
[6] 王文, 李永. 一种新非单调法求解压缩感知问题[J]. 电子科技, 2015, 28(2): 14-17.
[7] Lojasiewicz, S. (1963) Une propriété topologique des sous-ensembles analytiques réels. In: Les Équations aux Dérivées Partielles, Éditions du Centre National de la Recherche Scientifique, Paris, 87-89.
[8] Kurdyka, K. (1998) On Gradients of Functions Definable in O-Minimal Structures. Annales de l’Institut Fourier, 48, 769-783.
https://doi.org/10.5802/aif.1638
[9] 潘少华, 刘燕. 一类零模正则复合函数的指数1/2的KL性质[J]. 辽宁师范大学学报(自然科学版), 2019, 42(1): 1-4.
[10] 张立卫, 吴佳, 张艺. 变分分析与优化[M]. 北京: 科学出版社, 2013.
[11] 张立卫. 锥约束优化[M]. 北京: 科学出版社, 2010.
[12] Attouch, H., Bolte, J., Redont, P., et al. (2010) Proximal Alternating Minimization and Projection Methods for Nonconvex Problems: An Approach Based on the Kurdyka-Lojasiewicz Inequality. Mathematics of Operations Research, 35, 438-457. ‘
[13] Bolte, J., Sabach, S. and Teboulle, M. (2014) Proximal Alternating Linearized Minimization for Nonconvex and Nonsmooth Problems. Mathematical Programming, 146, 459-494.
https://doi.org/10.1007/s10107-013-0701-9
[14] Bonettini, S., Loris, I., Porta, F., et al. (2017) On the Convergence of a Linesearch Based Proximal-Gradient Method for Nonconvex Optimization. Inverse Problems, 33, Article ID: 055005.
https://doi.org/10.1088/1361-6420/aa5bfd
[15] Qiu, Y.Y., Yan, J.L. and Xu, F.Y. (2014) Nonmonotone Adaptive Barzilai-Borwein Gradient Algorithm for Compressed Sensing. Abstract and Applied Analysis, 2014, Article ID: 410104.
https://doi.org/10.1155/2014/410104
[16] 杨军, 孙世发, 马环宇, 等. 一种求解压缩感知问题的新方法[J]. 临沂大学学报, 2019, 41(6): 130-138.