线性规划之基追踪_python实现

非数学专业，本文不做理论推导及证明。
可进行指正说明探讨！！！
实时更新

最优化算法目录_python

---

一、线性规划

1.1 基追踪准则(BP)

解决重构压缩感知中的稀疏优化问题，实现对其$l_1$范数的最小化问题。

1.问题重述

线性方程组：
$$Ax=b$$
其中向量$x\in R^{n\times 1}，b\in R^{m\times 1}$，矩阵$A属于R^{m\times n}$，并且向量b的维数远远小于向量x的维数，即为$m\ll n$，因此方程(1)是欠正定的，存在无穷多个解，为“稀疏解”，即原始信号中有较多的零元素。(压缩感知)
$$x=[0,0,1,0,...,1,0,0]$$

2.构造$l_0$范数

构造一个128*256矩阵A，每个元素都服从高斯分布。

精确解u只有10%的元素非零，每个非零元素也服从高斯分布。以上理论是为了确保u是方程组(1)唯一的非零元素最少的解，即u是$l_0$范数问题的最优解：
$$\begin{gathered} min||x|{|}_0, \\ s.t.Ax=b. \end{gathered}$$
其中$||x|{|}_0$是指x中的非零元素的个数。由于$||x|{|}_0$是不连续的函数，并且只能取整数，故问题3实际上为非确定多项式(non-deterministic polynomial,NP)，求解困难。

因此定义$l_1$范数，$||x|{|}_1={\sum }_{i=1}^n|x{|}_i$，将上式替换到问题(3)中，得到等价于问题(1)的问题（称为$l_1$范数优化问题，基追踪问题）：
$$\begin{gathered} min||x|{|}_1, \\ s.t.Ax=b. \end{gathered}$$
通过证明，u也是$l_1$范数优化问题(4)的唯一最优解。

3.$l_1$范数最小化问题转化为线性规划问题

针对问题(4)，利用基追踪的定义，将问题(4)与线性规划(LP)联系起来。添加变量$\alpha$，将其定义两个非负变量的差：
$$\alpha =u-v(u,v\ge 0)$$
根据约束条件$Ax=b$，可以将公式(5)改写成$A(u-v)=b$，即为：
$$[A,-A]\left[\begin{array}{l}u\\ v\end{array}\right]=b$$
然后，根据范数的定义，目标函数可改写为：
$$||\alpha |{|}_1=\sum _{i=1}^n|\alpha {|}_i=\sum _{i=1}^n|u_i-v_i|$$
目标函数中包含绝对值，根据证明有，对于有限维的向量范数$l_1$范数最小化，即$min||\alpha |{|}_1={\sum }_{i=1}^n|\alpha {|}_i$，可等价为下属线性规划问题：
$$\begin{gathered} mine(u+v) \\ u-v=\alpha \\ u,v\ge 0 \end{gathered}$$
其中$e$为1行n列元素全为1的向量。定义$x\in R_m$，则可得到标准式的线性方程：
$$\begin{gathered} min{c}^{T}x, \\ s.t.Ax=b, \\ x\ge 0 \end{gathered}$$
公式(9)中，$c^{T}x$是目标函数，$Ax=b$是一组等式约束，$x\ge 0$定义边界。根据公式(5)、(6)，对上式中的变量做出如下映射，使其满足公式(8)：
$$m\iff 2n;x\iff \left[\begin{array}{l}u\\ v\end{array}\right];c\iff [1,1,,..,1{]}_{1\times 2n};A\iff [A,-A]$$
根据公式(9)和(10)，求解线性规划得到解$x_0=\left[\begin{array}{l}u\\ v\end{array}\right]$，则问题(5)的解为：
$${\alpha }_0=x_0(1:n)-x_0(n+1:2n)$$

4.基于linprog的基追踪Python代码

import numpy as np
from scipy import optimize as op # 线性规划库
from scipy.sparse import coo_matrix # 稀疏矩阵
import matplotlib.pyplot as plt
from matplotlib import rcParams # 绘图参数

config = {
    "font.family": 'serif',
    "font.size": 12,  #字体大小12，即小四
    "pgf.rcfonts": False,
    "text.usetex": True,
    "pgf.preamble": [
        r"\usepackage{unicode-math}",
        r"\setmainfont{Times New Roman}",
        r"\usepackage{xeCJK}",
        r"\setCJKmainfont{SimSun}",
    ],
}
rcParams.update(config)
# 坐标轴的刻度设置向内(in)或向外(out)
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

def sprase_rand(m, n, p):
    '''
    产生一个m行n列的服从均匀分布的随机稀疏矩阵
    :param m: 行数
    :param n: 列数
    :param p: 非零元素在此矩阵中所占的比例
    :return:  ss.coo_matric对象。使用toarray()方法或todense()方法转换为matrix对象或是array对象
    '''
    if type(m) != int or type(n) != int:
        raise TypeError('Rows(m) and Columns(n) must be an interger!')
    if p <= 0 or p > 1:
        raise ValueError('p must in (0, 1] !')
    # Counts of non-zero elements in sparse-matrix
    count = int(m * n * p) # 非零元素的个数
    # 随机生成服从高斯分布的行与列
    # Indexs of non-zero elements in sparse-matrix
    rows = np.random.randint(0, m, count)
    cols = np.random.randint(0, n, count)
    # 把非零值data按着行坐标rows和纵坐标cols写成两个列表
    # Values of non-zero elements in sparse-matrix
    data = np.random.randn(len(rows)) # 创建服从高斯分布(Gaussian distribution)的矩阵
    return coo_matrix((data, (rows, cols)), shape=(m, n))


def BP_linprog(Phi, s):
    '''
	L1范数求解最优解
    :param Phi: A矩阵
    :param s: b向量
    :return: L1范数下的最优解 alpha
    '''
    s,Phi = np.array(s),np.array(Phi) # 转化数据类型
    if np.size(s,0) < np.size(s,1): # 调整s向量，确保为列向量
        s = s.T
    p = np.size(Phi,1) # 提取A的列长
    # 根据公式10求解c，A_new，b_new
    # according to section 3.1 of the reference
    c = np.ones([2*p,1])
    A_new = np.hstack([Phi,-Phi])
    b_new = s
    bounds = [(0,None) for i in range(2*p)]
    # 求解线性规划方程
    x0 = op.linprog(c, A_eq=A_new, b_eq=b_new, bounds=bounds,method='revised simplex')['x']
    alpha = x0[:p] - x0[p:] # 带入公式求解
    return np.array([alpha])

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

alpha = BP_linprog(A, b) # 求解alpha
# 残差：实际值与预测值之间差的平方之和，残差平方和越小,其拟合程度越好。
print('\n残差：', np.linalg.norm(alpha.T - u.toarray(), 2))
# 绘图
plt.figure(figsize = (8,6))
plt.xlabel('length')
plt.ylabel('RSS')
plt.plot(u.toarray(), 'r', linewidth = 1, label = 'Original') # 绘制原信号
plt.plot(alpha.T, 'b--', linewidth = 1, label = 'Recovery') # 绘制恢复信号
plt.grid()
plt.legend()
plt.show()
 # 绘制原信号
plt.plot(alpha.T, 'b--', linewidth = 1, label = 'Recovery') # 绘制恢复信号
plt.grid()
plt.legend()
plt.show()

输出结果：

  	残差： 8.146147562440153e-15

参考文献

最优化：建模、算法与理论/最优化计算方法
python创建稀疏矩阵
在Python中创建、生成稀疏矩阵（均匀分布、高斯分布）
稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现