一、前言

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

二、问题重述

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

Ax = b \tag{1}

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

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

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

三、构造 $\ell_1$ 范数

举一个具体的例子(在 Python 环境里构造 $A, u$ 和 $b$ )^[刘浩洋，户将，李勇峰，文再文. 最优化：建模、算法与理论[M]. 北京: 高等教育出版社, 2020: 4-11.]：

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 \times 256$ 矩阵 $A$ ，它的每个元素都服从高斯 (Gauss) 随机分布(参照我的这篇博客：在Python中创建、生成稀疏矩阵（均匀分布、高斯分布）)。

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

\min \left\|x\right\|_0,\ s.t.\ Ax = b. \tag{3}

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

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

\min \left\|x\right\|_1,\ s.t.\ Ax = b. \tag{4}

可以从理论上证明：若 $A, b$ 满足一定的条件^[DONOHO D L. Compressedsensing[J/OL]. IEEE Transactions on information theory, 2006(4): 1289-1306. www.signallake.com/innovation/Compre ssedSensing091604.pdf.](例如使用前面随机产生的 $A$ 和 $b$ )，向量 $u$ 也是 $\ell_1$ 范数优化问题 (4) 的唯一最优解。

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

在文献^[DONOHO D L, CHEN S S, SAUNDERS M A. Atomicdecomposition by basis pursuit[J/OL]. SIAM review, 2001(1): 129-159. citeseerx.ist.psu.edu/viewdoc/dow….]中，给出了将问题 (4) 即 $\ell_1$ 范数最小化问题转换为标准的线性规划问题的一种方法，首先对于如下问题：

\min \left\|\alpha\right\|_1,\ s.t.\ \Phi\alpha = s. \tag{P}

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

\alpha\ =\ u - v,\ u, v \ge 0 \tag{5}

由于 $u$ 可以大于也可以小于 $v$ ，所以 $a$ 可以是正的也可以是负的^[朱德通，孙文瑜，徐成贤. 最优化方法 (第二版)[M]. 北京: 高等教育出版社，2010: 49-51.]。接着约束条件 $\Phi\alpha = s$ 重写为 $\Phi(u - v) = s$ ，也可改写为如下形式：

\begin{bmatrix} \Phi,-\Phi \end{bmatrix} \begin{bmatrix} u\\v \end{bmatrix}\ =\ s \tag{6}

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

\left\|\alpha\right\|_1\ =\ \sum_{i=1}^{n}\left|\alpha_i\right|\ =\ \sum_{i=1}^{n}\left|u_i-v_i\right| \tag{7}

目标函数中包含绝对值，采用网页^[L1 范数优化的线性化方法如何证明?[J/OL]. www.zhihu.com/question/21….]中的证明，对于有限维度的向量 $\ell_1$ 范数最小化，即 $\min\left\|\alpha\right\|_1 = \sum_{i=1}^{n} \left|\alpha_i\right|$ ，等价于如下线性规划问题：

\left\{\begin{aligned} min\ \ \ \ &e(u\ +\ v)\\ &u\ -\ v\ =\ \alpha\\ &u, v \ge 0 \end{aligned}\right. \tag{8}

其中 $e$ 为 1 行 $n$ 列元素全为 1 的向量。接着若定义 $x \in R^m$ ，则可得到如下标准形式的线性方程：

\min\ c^Tx,\ s.t.\ Ax = b,\ x \ge 0. \tag{9}

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

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}

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

\alpha_0 = x_0(1:n) - x_0(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 : vector
        Optimal 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.T
    p = np.size(Phi, 1)
    # according to section 3.1 of the reference
    c = np.ones([2 * p, 1])
    Aeq = np.hstack([Phi, -Phi])
    beq = s
    bounds = [(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 比较也有挺不错的实现价值。

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