文章目录

一、前言
二、问题重述
三、构造 ℓ1\ell_1ℓ1 范数
四、ℓ1\ell_1ℓ1 范数最小化问题转换为线性规划问题
五、基于linprog的基追踪Python代码
六、运行测试
七、总结

一、前言

本文针对压缩重构感知中的稀疏优化问题，实现了对其 L1 范数最小化问题 的求解，文章内容较长，请耐心看完，代码部分在本文第五章。

二、问题重述

考虑线性方程组求解问题：

Ax=b(1)Ax = b \tag{1} Ax=b(1)

其中向量 x∈Rn×1,b∈Rm×1x \in R^{n\times 1},\ b \in R^{m\times 1}x∈Rn×1, b∈Rm×1，矩阵 A∈Rm×nA \in R^{m\times n}A∈Rm×n，且向量 bbb 的维数远小于向量 xxx 的维数，即 m≪nm \ll nm≪n。由于 m≪nm \ll nm≪n，方程组 (1) 是欠定的，因此存在无穷多个解，但是真正有用的解是所谓的“稀疏解”，即原始信号中有较多的零元素。

x=[0,0,1,0,...,1,0,0](2)x = [0, 0, 1, 0, ..., 1, 0, 0] \tag{2} x=[0,0,1,0,...,1,0,0](2)

如果加上稀疏性这一先验信息，且矩阵 AAA 以及原问题的解 uuu 满足某些条件，那么我们可以通过求解稀疏优化问题把 uuu 与方程组 (1) 的其他解区别开。这类技术广泛应用于压缩感知(compressive sensing)，即通过部分信息恢复全部信息的解决方案。

三、构造 ℓ1\ell_1ℓ1 范数

举一个具体的例子(在 Python 环境里构造 A,uA, uA,u 和 bbb)¹：

import numpy as np
m, n = 128, 256
# 128x256矩阵，每个元素服从Gauss随机分布
A = np.random.randn(m, n)
# 精确解 u 只有 10% 元素非零，每一个非零元素也服从高斯分布
# 可保证 u 是方程组唯一的非零元素最少的解
u = sprase_rand(n, 1, 0.1)
b = A * u

在这个例子中，我们构造了一个 128×256128 \times 256128×256 矩阵 AAA，它的每个元素都服从高斯 (Gauss) 随机分布(参照我的这篇博客：在Python中创建、生成稀疏矩阵（均匀分布、高斯分布）)。

精确解 uuu 只有 10% 的元素非零，每一个非零元素也服从高斯分布。这些特征可以在理论上保证 uuu 是方程组 (1) 唯一的非零元素最少的解，即 uuu 是如下 ℓ0\ell_0ℓ0 范数问题的最优解：

min⁡∥x∥0,s.t.Ax=b.(3)\min \left\|x\right\|_0,\ s.t.\ Ax = b. \tag{3} min∥x∥0, s.t. Ax=b.(3)

其中 ∥x∥0\left\|x\right\|_0∥x∥0是指 xxx 中非零元素的个数．由于 ∥x∥0\left\|x\right\|_0∥x∥0 是不连续的函数，且取值只可能是整数，问题 (3) 实际上是 NP(non-deterministic polynomial) 难的，求解起来非常困难。

若定义 ℓ1\ell_1ℓ1 范数：∥x∥1=∑i=1n∣xi∣\left\|x\right\|_1 = \sum_{i=1}^{n} \left|x_i\right|∥x∥1=∑i=1n∣xi∣，并替换到问题 (3) 中，我们得到了另一个形式上非常相似的问题(又称 ℓ1\ell_1ℓ1 范数优化问题，基追踪问题)：

min⁡∥x∥1,s.t.Ax=b.(4)\min \left\|x\right\|_1,\ s.t.\ Ax = b. \tag{4} min∥x∥1, s.t. Ax=b.(4)

可以从理论上证明：若 A,bA, bA,b 满足一定的条件²(例如使用前面随机产生的 AAA 和 bbb)，向量 uuu 也是 ℓ1\ell_1ℓ1 范数优化问题 (4) 的唯一最优解。

四、ℓ1\ell_1ℓ1 范数最小化问题转换为线性规划问题

在文献³中，给出了将问题 (4) 即 ℓ1\ell_1ℓ1 范数最小化问题转换为标准的线性规划问题的一种方法，首先对于如下问题：

min⁡∥α∥1,s.t.Φα=s.(P)\min \left\|\alpha\right\|_1,\ s.t.\ \Phi\alpha = s. \tag{P} min∥α∥1, s.t. Φα=s.(P)

利用基追踪 (Basis Pursuit) 的定义，我们尝试将问题 P 与线性规划 (LP) 问题连接起来。首先，式 P 中变量 α\alphaα 没有非负约束，在此我们将 α\alphaα 定义为两个非负变量的差：

α=u−v,u,v≥0(5)\alpha\ =\ u - v,\ u, v \ge 0 \tag{5} α = u−v, u,v≥0(5)

由于 uuu 可以大于也可以小于 vvv，所以 aaa 可以是正的也可以是负的⁴。接着约束条件 Φα=s\Phi\alpha = sΦα=s 重写为 Φ(u−v)=s\Phi(u - v) = sΦ(u−v)=s，也可改写为如下形式：

[Φ,−Φ][uv]=s(6)\begin{bmatrix} \Phi,-\Phi \end{bmatrix} \begin{bmatrix} u\\v \end{bmatrix}\ =\ s \tag{6} [Φ,−Φ][uv] = s(6)

然后，根据范数的定义，目标函数改写为：

∥α∥1=∑i=1n∣αi∣=∑i=1n∣ui−vi∣(7)\left\|\alpha\right\|_1\ =\ \sum_{i=1}^{n}\left|\alpha_i\right|\ =\ \sum_{i=1}^{n}\left|u_i-v_i\right| \tag{7} ∥α∥1 = i=1∑n∣αi∣ = i=1∑n∣ui−vi∣(7)

目标函数中包含绝对值，采用网页⁵中的证明，对于有限维度的向量 ℓ1\ell_1ℓ1 范数最小化，即 min⁡∥α∥1=∑i=1n∣αi∣\min\left\|\alpha\right\|_1 = \sum_{i=1}^{n} \left|\alpha_i\right|min∥α∥1=∑i=1n∣αi∣，等价于如下线性规划问题：

{mine(u+v)u−v=αu,v≥0(8)\left\{\begin{aligned} min\ \ \ \ &e(u\ +\ v)\\ &u\ -\ v\ =\ \alpha\\ &u, v \ge 0 \end{aligned}\right. \tag{8} ⎩⎪⎨⎪⎧min e(u + v)u − v = αu,v≥0(8)

其中 eee 为 1 行 nnn 列元素全为 1 的向量。接着若定义 x∈Rmx \in R^mx∈Rm，则可得到如下标准形式的线性方程：

min⁡cTx,s.t.Ax=b,x≥0.(9)\min\ c^Tx,\ s.t.\ Ax = b,\ x \ge 0. \tag{9} min cTx, s.t. Ax=b, x≥0.(9)

在式 (9) 中，cTxc^TxcTx 是目标函数，Ax=bAx = bAx=b 是一组等式约束，x≥0x \ge 0x≥0 定义了边界。结合式 (5)、(6)，我们对式 (9) 中的各变量做出如下映射使其符合式 (8)：

m⇔2n;x⇔[uv];c⇔[1,1,...,1]1×2n;A⇔[Φ,−Φ];b⇔s.(10)m\Leftrightarrow2n;\ x\Leftrightarrow\begin{bmatrix}u\\v\end{bmatrix};\ c\Leftrightarrow[1, 1, ..., 1]_{1\times 2n};\ A\Leftrightarrow[\Phi, -\Phi];\ b\Leftrightarrow s. \tag{10} m⇔2n; x⇔[uv]; c⇔[1,1,...,1]1×2n; A⇔[Φ,−Φ]; b⇔s.(10)

根据式 (9) 和式 (10) 求解线型规划得到解 x0=[uv]x_0 = \begin{bmatrix}u\\v\end{bmatrix}x0=[uv]，则问题 (P)(P)(P) 的解

α0=x0(1:n)−x0(n+1:2n)\alpha_0 = x_0(1:n) - x_0(n+1:2n) α0=x0(1:n)−x0(n+1:2n)

五、基于linprog的基追踪Python代码

同理可以用 MATLAB 做，原理在上面了，MATLAB 的可参考这篇博客：压缩感知重构算法之基追踪(Basis Pursuit, BP)

import numpy as np
from scipy import optimize as op def BP_linprog(Phi, s):'''s = Phi * alpha(Given Phi & s, try to derive alpha, alpha is a sparse vector)Parameters----------Phi : A matrix.s : A vector.Returns-------alpha : vectorOptimal solutions of equations under L1 norm.Reference---------Chen S S, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit[J]. SIAM review, 2001, 43(1): 129-159.(Available at: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf)Version: 1.0 writen by z.q.feng @2022.03.13'''s, Phi = np.array(s), np.array(Phi)if np.size(s, 0) < np.size(s, 1):s = s.Tp = np.size(Phi, 1)# according to section 3.1 of the referencec = np.ones([2 * p, 1])Aeq = np.hstack([Phi, -Phi])beq = sbounds = [(0, None) for i in range(2 * p)]x0 = op.linprog(c, A_eq = Aeq, b_eq = beq, bounds = bounds,method='revised simplex')['x']alpha = x0[:p] - x0[p:]return np.array([alpha])

六、运行测试

根据第二节的代码产生的矩阵来测试，编写运行绘图代码如下：

alpha = BP_linprog(A, b)
# 恢复残差
print('\n恢复残差：', np.linalg.norm(alpha.T - u.toarray(), 2))
# 绘图
plt.figure(figsize = (8, 6))
# 绘制原信号
plt.plot(u.toarray(), 'r', linewidth = 2, label = 'Original')
# 绘制恢复信号
plt.plot(alpha.T, 'b--', linewidth = 2, label = 'Recovery')
plt.grid()
plt.legend()
plt.show()

得到输出如下：

恢复残差： 4.0645022580106625e-15

绘制图像如下：

七、总结

Python 中的的科学计算库让其与 MATLAB 比较也有挺不错的实现价值。

刘浩洋，户将，李勇峰，文再文. 最优化：建模、算法与理论[M]. 北京: 高等教育出版社, 2020: 4-11. ↩︎
DONOHO D L. Compressedsensing[J/OL]. IEEE Transactions on information theory, 2006(4): 1289-1306. http://www.signallake.com/innovation/Compre
ssedSensing091604.pdf. ↩︎
DONOHO D L, CHEN S S, SAUNDERS M A. Atomicdecomposition by basis pursuit[J/OL]. SIAM review, 2001(1): 129-159. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.4272&rep=rep1&type=pdf. ↩︎
朱德通，孙文瑜，徐成贤. 最优化方法 (第二版)[M]. 北京: 高等教育出版社，2010: 49-51. ↩︎
L1 范数优化的线性化方法如何证明?[J/OL]. http://www.zhihu.com/question/21427075. ↩︎

稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现相关推荐

l1范数最小化快速算法【文献阅读】
1:解决的问题模型如下: 或者约束条件可以适当的松弛,即为如下模型: 当然约束条件取范数,数据获取的比较准确,结果逼近的效果更好,防止过拟合.如果取范数,则是获取的数据,受到污染比较严重. ...
l1范数最小化快速算法
1:解决的问题模型如下: 或者约束条件可以适当的松弛,即为如下模型: 当然约束条件取
超定线性方程组Ax=b极小L1范数求解——MATLAB/Python实现
文章目录一.前言二.问题重述二.极小模剩余向量的性质及求法三.基于基追踪准则的一种求解算法四.算法伪码五.超定线性方程组极小 ℓ1\ell_1ℓ1 范数求解Python代码六.检验与测 ...
l20范数最小化求解系数方程_贪婪组稀疏方法(Greedy group sparsity)
l20范数最小化求解系数方程_贪婪组稀疏方法(Greedy group sparsity) 本文章部分参考Fast group sparse classification l20范数最小化求解系数方程 ...
压缩感知的尽头: 原子范数最小化
文章目录前言问题建模 Toeplitz 矩阵的范德蒙德分解 DOA估计的一般框架 ℓ0\ell_0ℓ0-原子范数 ℓ0\ell_0ℓ0-原子范数与范德蒙德分解原子范数多维原子范数证明 ...
三维网格去噪算法（L0范数最小化，包含二维图像去噪）
参考文章(技术来源):Mesh denoising via L0 minimization 上面参考文章提出了一种基于L0范数最小化的三角网格去噪算法.该思想由二维图像平滑引申而来,所以先从基于L0范 ...
K-means的缺点（优化不仅仅是最小化误差）
K-means的缺点(优化不仅仅是最小化误差) #转载时,请注明英文原作David Robinson,译者Ding Chao.# 我最近遇到一个交叉验证的问题,我认为这个给我提供了一个很好的机会去用& ...
图像处理-基于图像梯度L0范数最小化(L0smooth)的保护边缘平滑滤波
算法程序备注: (1)下面是对一幅自然图像进行处理的结果: 可以看到图像有非常明显的变化,图像分成了一块一块,这是图像平滑后的结果,因为保护了边界,因此明显的边界仍然存在,但是不可避免的细节部分被磨平 ...
压缩感知重构算法之基追踪(Basis Pursuit, BP).基追踪并不能称为一个具体的算法，而是一种最优化准则,可以有很多实现方式,我认为指的是L0可以变为L1的准则
基追踪(basis pursuit)算法是一种用来求解未知参量L1范数最小化的等式约束问题的算法. 基追踪是通常在信号处理中使用的一种对已知系数稀疏化的手段.将优化问题中的L0范数转化为L1范数的求解 ...

稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现

文章目录

一、前言

二、问题重述

三、构造 ℓ1\ell_1ℓ1 范数

四、ℓ1\ell_1ℓ1 范数最小化问题转换为线性规划问题

五、基于linprog的基追踪Python代码

六、运行测试

七、总结

稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现相关推荐

最新文章

热门文章

稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现

文章目录

一、前言

二、问题重述

三、构造 ℓ1\ell_1ℓ1​ 范数

四、ℓ1\ell_1ℓ1​ 范数最小化问题转换为线性规划问题

五、基于linprog的基追踪Python代码

六、运行测试

七、总结

稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现相关推荐

最新文章

热门文章

三、构造 ℓ1\ell_1ℓ1 范数

四、ℓ1\ell_1ℓ1 范数最小化问题转换为线性规划问题