超定线性方程组Ax=b极小L1范数求解——MATLABPython实现_python

文章目录

一、前言
二、问题重述
二、极小模剩余向量的性质及求法
三、基于基追踪准则的一种求解算法
四、超定线性方程组极小 ℓ 1 \ell_1 ℓ1 范数求解Python代码
五、检验与测试
- 1. 矩阵A满秩情况
- 2. 矩阵A列缺秩情况
六、总结

一、前言

本文前半部分算法实现参考的为姚健康老师的论文，这篇论文在万方学术平台可以免费下载：http://d.g.wanfangdata.com.cn/Periodical_jxkx200701002.aspx

后半部分算法涉及到稀疏优化，我在之前的博客中介绍了实现原理并给出了代码，链接如下：稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现

二、问题重述

给定超定线性方程组 A x = b Ax=b Ax=b：

A = [ a 11 a 12 . . . a 1 n a 21 a 22 . . . a 2 n . . . . . . . . . . . . a m 1 a m 2 . . . a m n ] ∈ R m × n , b = [ b 1 b 2 ⋮ b m ] ∈ R m × 1 A=\begin{bmatrix} a_{11}&a_{12}&...&a_{1n}\ a_{21}&a_{22}&...&a_{2n}\ ...&...&...&...\ a_{m1}&a_{m2}&...&a_{mn}\ \end{bmatrix}\in R^{m\times n},\ b=\begin{bmatrix} b_1\b_2\\vdots\b_m \end{bmatrix}\in R^{m\times 1} A=⎣⎢⎢⎡a11a21...am1a12a22...am2............a1na2n...amn⎦⎥⎥⎤∈Rm×n, b=⎣⎢⎢⎢⎡b1b2⋮bm⎦⎥⎥⎥⎤∈Rm×1

其中 m > n ≥ 2 m>n\ge 2 m>n≥2，求 x = [ x 1 , x 2 , . . . , x n ] T ∈ R n × 1 x={[x_1,x_2,...,x_n]}^T\in R^{n\times 1} x=[x1,x2,...,xn]T∈Rn×1，使得

min ⁡ ∥ A x − b ∥ 1 = min ⁡ ∑ i = 1 m ∣ ∑ j = 1 n a i j x j − b i ∣ (1) \min \left\|Ax-b\right\|_1=\min \sum_{i = 1}^{m}|\sum_{j = 1}^{n}a_{ij}x_j-b_i| \tag{1} min∥Ax−b∥1=mini=1∑m∣j=1∑naijxj−bi∣(1)

即求超定线性方程组 A x = b Ax=b Ax=b 在 ℓ 1 \ell_1 ℓ1 范数意义下的解，简称 ℓ 1 \ell_1 ℓ1 范数极小化问题¹。

定义：设 x ∗ x^* x∗ 为问题 ( 1 ) (1) (1) 的解，称 r ∗ = A x ∗ − b {r}^*={A}x^* - {b} r∗=Ax∗−b 为问题 ( 1 ) (1) (1) 的极小 ℓ 1 \ell_1 ℓ1 模剩余向量。

如果能求得极小 ℓ 1 \ell_1 ℓ1 模剩余向量 r ∗ {r}^* r∗，那么相容线性方程组 A x = b + r ∗ {A}x = {b} + {r}^* Ax=b+r∗ 的解即为问题 ( 1 ) (1) (1) 的解。

二、极小模剩余向量的性质及求法

依据文献²中的说明，我们得到如下定理：

定理：设有 2 个超定线性方程组 A x = b , B x = b {A}x={b}, {B}x={b} Ax=b,Bx=b，若存在可逆阵 P {P} P，使 B = A P {B} = {A}{P} B=AP，则这 2 个超定线性方程组的极小 ℓ 1 \ell_1 ℓ1 模剩余向量相同。

设 A = [ A 1 A 2 ] {A}=\begin{bmatrix} {A_1}\{A_2} \end{bmatrix} A=[A1A2]( A 1 {A_1} A1 为 n n n 阶可逆阵)，则方程 A x = b {A}x={b} Ax=b 可以转变为式 ( 2 ) (2) (2)：

[ I A 2 A 1 − 1 ] A 1 x = b (2) \begin{bmatrix} {I}\{A_2{A_1}^{-1}} \end{bmatrix}{A_1}x={b} \tag{2} [IA2A1−1]A1x=b(2)

设方程组 A x = b {A}x={b} Ax=b 的极小 ℓ 1 \ell_1 ℓ1 模剩余向量为 r = [ r 1 , r 2 , . . . , r m ] T {r}={[r_1, r_2, ..., r_m]}^T r=[r1,r2,...,rm]T，则由上述定理知，式 ( 2 ) (2) (2) 与方程组 A x = b {A}x={b} Ax=b 具有相同的极小 ℓ 1 \ell_1 ℓ1 模剩余向量，因而 r {r} r 也是式 ( 2 ) (2) (2) 的极小 ℓ 1 \ell_1 ℓ1 模剩余向量，记 A 1 x = y {A_1}x={y} A1x=y，故

[ I A 2 A 1 − 1 ] y = b + r (3) \begin{bmatrix} {I}\{A_2{A_1}^{-1}} \end{bmatrix}{y}={b} + {r} \tag{3} [IA2A1−1]y=b+r(3)

若我们令

A 2 A 1 − 1 = [ c 11 c 12 . . . c 1 n c 21 c 22 . . . c 2 n . . . . . . . . . . . . c m − n , 1 c m − n , 2 . . . c m − n , n ] ( m − n ) × n (4) {A_2{A_1}^{-1}}=\begin{bmatrix} c_{11}&c_{12}&...&c_{1n}\ c_{21}&c_{22}&...&c_{2n}\ ...&...&...&...\ c_{m-n,1}&c_{m-n,2}&...&c_{m-n,n}\ \end{bmatrix}_{(m-n)\times n} \tag{4} A2A1−1=⎣⎢⎢⎡c11c21...cm−n,1c12c22...cm−n,2............c1nc2n...cm−n,n⎦⎥⎥⎤(m−n)×n(4)

则极小 ℓ 1 \ell_1 ℓ1 模剩余向量 r r r 满足如下线性方程组：

{ r n + 1 − ( c 11 r 1 + c 12 r 2 + ⋯ + c 1 n r n ) = d 1 r n + 2 − ( c 21 r 1 + c 22 r 2 + ⋯ + c 2 n r n ) = d 2 ⋯ ⋯ r m − ( c m − n , 1 r 1 + c m − n , 2 r 2 + ⋯ + c m − n , n r n ) = d m − n (5) \left\{\begin{aligned} &r_{n+1}-(c_{11}r_1+c_{12}r_2+\cdots+c_{1n}r_n)=d_1&\ &r_{n+2}-(c_{21}r_1+c_{22}r_2+\cdots+c_{2n}r_n)=d_2&\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots\cdots&\ &r_{m}-(c_{m-n,1}r_1+c_{m-n,2}r_2+\cdots+c_{m-n,n}r_n)=d_{m-n}&\ \end{aligned}\right. \tag{5} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧rn+1−(c11r1+c12r2+⋯+c1nrn)=d1rn+2−(c21r1+c22r2+⋯+c2nrn)=d2 ⋯⋯rm−(cm−n,1r1+cm−n,2r2+⋯+cm−n,nrn)=dm−n(5)

其中

A 2 A 1 − 1 [ b 1 b 2 ⋮ b n ] − [ b n + 1 b n + 2 ⋮ b m ] = [ d 1 d 2 ⋮ d m − n ] (6) {A_2{A_1}^{-1}}\begin{bmatrix} b_1\b_2\\vdots\b_n \end{bmatrix} - \begin{bmatrix} b_{n+1}\b_{n+2}\\vdots\b_m \end{bmatrix}=\begin{bmatrix} d_1\d_2\\vdots\d_{m-n} \end{bmatrix} \tag{6} A2A1−1⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤−⎣⎢⎢⎢⎡bn+1bn+2⋮bm⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡d1d2⋮dm−n⎦⎥⎥⎥⎤(6)

具体转换过程请参考姚老师的论文。

于是要求极小 ℓ 1 \ell_1 ℓ1 模剩余向量 r {r} r，即求满足式 ( 5 ) (5) (5) 的 { r 1 , r 2 , . . . , r m } \{r_1, r_2, ..., r_m\} {r1,r2,...,rm}，使得 ∑ i = 1 m ∣ r i ∣ \sum_{i=1}^{m}|r_i| ∑i=1m∣ri∣ 最小。

三、基于基追踪准则的一种求解算法

由式 ( 5 ) (5) (5) 可得，原超定线性方程组 A x = b {A}x={b} Ax=b 的极小 ℓ 1 \ell_1 ℓ1 模剩余向量可由如下约束不可微优化问题求得：

{ m i n ∑ i = 1 m ∣ z i ∣ s . t . z n + 1 − ( c 11 z 1 + c 12 z 2 + ⋯ + c 1 n z n ) = d 1 , z n + 2 − ( c 21 z 1 + c 22 z 2 + ⋯ + c 2 n z n ) = d 2 , ⋯ ⋯ z m − ( c m − n , 1 z 1 + c m − n , 2 z 2 + ⋯ + c m − n , n z n ) = d m − n (P) \left\{\begin{aligned} &min\sum_{i=1}^{m}|z_i|\ s.t.\ \ \ \ &z_{n+1}-(c_{11}z_1+c_{12}z_2+\cdots+c_{1n}z_n)=d_1,\ &z_{n+2}-(c_{21}z_1+c_{22}z_2+\cdots+c_{2n}z_n)=d_2,\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \cdots\cdots\ &z_{m}-(c_{m-n,1}z_1+c_{m-n,2}z_2+\cdots+c_{m-n,n}z_n)=d_{m-n} \end{aligned}\right. \tag{P} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧s.t. mini=1∑m∣zi∣zn+1−(c11z1+c12z2+⋯+c1nzn)=d1,zn+2−(c21z1+c22z2+⋯+c2nzn)=d2, ⋯⋯zm−(cm−n,1z1+cm−n,2z2+⋯+cm−n,nzn)=dm−n(P)

记问题 ( P ) (P) (P) 的最优解 z ∗ = [ z 1 ∗ , z 2 ∗ , . . . , z m ∗ ] T {z}^*={[z_1^*, z_2^*, ..., z_m^*]}^T z∗=[z1∗,z2∗,...,zm∗]T，方程组 A x = b + z ∗ {A}x={b}+{z}^* Ax=b+z∗ 为原方程的一个相容线性方程组，其解 x ∗ {x}^* x∗ 即为问题 ( 1 ) (1) (1) 的解：

x ∗ = A 1 − 1 [ b 1 b 2 ⋮ b n ] + A 1 − 1 [ z 1 ∗ z 2 ∗ ⋮ z n ∗ ] (7) {x}^*={{A_1}^{-1}}\begin{bmatrix} b_1\\b_2\\\vdots\\b_n \end{bmatrix} + {{A_1}^{-1}}\begin{bmatrix} z_1^*\\z_2^*\\\vdots\\z_n^* \end{bmatrix} \tag{7} x∗=A1−1⎣⎢⎢⎢⎡b1b2⋮bn⎦⎥⎥⎥⎤+A1−1⎣⎢⎢⎢⎡z1∗z2∗⋮zn∗⎦⎥⎥⎥⎤(7)

在问题 ( P ) (P) (P) 中，若我们定义 C ∈ R ( m − n ) × m {C} \in {R}^{(m-n)\times m} C∈R(m−n)×m：

C = [ − c 11 − c 12 ⋯ − c 1 n 1 0 ⋯ 0 − c 21 − c 22 ⋯ − c 2 n 0 1 ⋯ 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ − c m − n , 1 − c m − n , 2 ⋯ − c m − n , n 0 0 ⋯ 1 ] (8) {C}=\begin{bmatrix} -c_{11}&-c_{12}&\cdots&-c_{1n}&1&0&\cdots&0\ -c_{21}&-c_{22}&\cdots&-c_{2n}&0&1&\cdots&0\ \cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\ -c_{m-n,1}&-c_{m-n,2}&\cdots&-c_{m-n,n}&0&0&\cdots&1\ \end{bmatrix} \tag{8} C=⎣⎢⎢⎡−c11−c21⋯−cm−n,1−c12−c22⋯−cm−n,2⋯⋯⋯⋯−c1n−c2n⋯−cm−n,n10⋯001⋯0⋯⋯⋯⋯00⋯1⎦⎥⎥⎤(8)

则问题 ( P ) (P) (P) 变成了如下问题：

m i n ∑ i = 1 m ∣ z i ∣ , s . t . C z = d . (Q) min\sum_{i=1}^{m}|z_i|,\ s.t.\ {C}{z}={d}. \tag{Q} mini=1∑m∣zi∣, s.t. Cz=d.(Q)

对于极小 ℓ 1 \ell_1 ℓ1 模剩余向量 z {z} z，考虑到我们的求解应当尽可能逼近精确解，即极小 ℓ 1 \ell_1 ℓ1 模剩余向量 z {z} z 中的元素应当尽可能为0。

这样问题 ( Q ) (Q) (Q) 的形式又回到了 ℓ 1 \ell_1 ℓ1 范数最小化问题中的稀疏优化问题。

参考我的这篇博客求解上述稀疏优化问题：稀疏优化L1范数最小化问题求解之基追踪准则(Basis Pursuit)——原理及其Python实现

四、超定线性方程组极小 ℓ 1 \ell_1 ℓ1 范数求解Python代码

采用 Python 中 scipy 模块内置的 linprog 函数求解线性方程组，编写代码如下：

import numpy as np 
from scipy import optimize as op

def OLE_linprog(A, b):
    '''
    Ax = b (Given A & b, try to derive x)
    
    Parameters
    ----------
    A : matrix like.
    b : array like.

    Returns
    -------
    x : array
        Minimal L1 norm solution of the system of equations.
        
    Reference
    ---------
    YAO Jian-kang. An Algorithm for Minimizing l1-Norm to Overdetermined Linear Eguations[J]. 
    JIANGXI SCE7NICE, 2007, 25(1): 1-4.(Available at: 
    http://d.g.wanfangdata.com.cn/Periodical_jxkx200701002.aspx)
        
    Version: 1.0 writen by z.q.feng @2022.03.13
    '''
    A, b = np.array(A), np.array(b)
    if np.size(A, 0) < np.size(A, 1):
        raise ValueError('Matrix A rows must greater than columns!')
    m, n = A.shape
    # Trans A into two matrix(n x n and (m - n) x n)
    A1, A2 = A[:n, :], A[n:, :]
    if np.linalg.matrix_rank(A) >= n:
        # inverse of A1 
        A1_ = np.linalg.inv(A1)
    else:
        # Generalized inverse of A1 
        A1_ = np.linalg.pinv(A1)
    # c_ij = A2 * A1_
    c = np.dot(A2, A1_)
    # r(n+1:m) = A2*inv(A1)*r(1:n) + d
    d = np.dot(c, b[:n]) - b[n:]
    # Basic-Pursuit, target function = sum(u, v)
    t = np.ones([2 * m, 1])
    # Aeq_ = [c I(m-n)]
    Aeq_ = np.hstack([-c, np.eye(m - n, m - n)])
    # Aeq[u v] = Aeq_ * (u - v)
    Aeq = np.hstack([Aeq_, -Aeq_])
    # u, v > 0
    bounds = [(0, None) for i in range(2 * m)]
    # r0 = [u; v]
    r0 = op.linprog(t, A_eq = Aeq, b_eq = d, bounds = bounds,
                    method='revised simplex')['x']
    # Minimal L1-norm residual vector, r = u - v
    r = np.array([r0[:m] - r0[m:]])
    # Solving compatible linear equations Ax = b + r
    # Generalized inverse solution
    x = np.linalg.pinv(A).dot(b + r.T)
    return x

采用 MATLAB 求解的由于篇幅问题我就不放了，私信我。

五、检验与测试 1. 矩阵A满秩情况

对于线性方程 A x = b {A}x = {b} Ax=b，我们通过产生给定的 A {A} A 与 u {u} u，并令 b = A ∗ u {b} = {A} * {u} b=A∗u，则可以认为 u {u} u 即为问题 min ⁡ ∥ A x − b ∥ 1 \min\ \left\|{A}x-{b}\right\|_1 min ∥Ax−b∥1 的最优解。

接着我们利用上述函数求解得到其经验解 u ∗ {u}^* u∗，利用 ∥ u ∗ − u ∥ 2 \left\|{u}^*-{u}\right\|_2 ∥u∗−u∥2 定义其恢复残差，并绘图确定求解准确性。

if __name__ == '__main__':
    m, n = 256, 128
    # 256x128矩阵，每个元素服从Gauss随机分布
    A = np.random.randn(m, n)
    # 128x1矩阵，每个元素服从Gauss随机分布
    # b = A * u, 则 Au - b = 0，u 即为我们要求的精确解
    u = np.random.randn(n, 1)
    b = np.dot(A, u)
    alpha = OLE_linprog(A, b)
    # 恢复残差 
    print('\n恢复残差：', np.linalg.norm(alpha - u, 2))
    # 绘图
    plt.figure(figsize = (8, 6))
    # 绘制原信号
    plt.plot(u, 'r', linewidth = 2, label = 'Original')
    # 绘制恢复信号
    plt.plot(alpha, 'b--', linewidth = 2, label = 'Recovery')
    plt.grid() 
    plt.legend()
    plt.show()

得到恢复残差为6.0407e-14，恢复信号图如图所示(信号随机生成，每次结果均不一样)。

2. 矩阵A列缺秩情况

根据式 ( 2 ) (2) (2) 可以看出，若矩阵 A {A} A 行不满秩，即 n < r a n k ( A 1 ) < m n < rank({A_1}) < m n<rank(A1)<m，我们的方法是不需要做对应改变的，而若矩阵 A {A} A 列缺秩，则需要使用矩阵 A 1 {A_1} A1 的Moore-Penrose逆矩阵 A 1 † {A_1^\dagger} A1† 替代 A 1 − 1 {{A_1}^{-1}} A1−1 在式 ( 2 ) (2) (2) 中构造 b {b} b。

下面我们通过代码让矩阵 A {A} A 列秩亏缺1：

A(:, n) = zeros(m, 1);

运行代码，得到恢复残差为0.2551，恢复信号图如图所示(信号随机生成，每次结果均不一样)。

可以看到，虽然误差增大了，但是还在允许范围内。

六、总结

目前，求解问题 ( 1 ) (1) (1) 的方法已有不少，主要有：

线性规划方法；
投影下降法；
有效集法；
区间方法等。

上述方法各有优缺点，如线性规划方法使规模扩大，投影下降算法过于复杂，不易被工程技术人员掌握，有效集法较其它方法直观，算法简单，区间方法过于复杂，无最优解区间的检验比较麻烦。

本文利用极小 ℓ 1 \ell_1 ℓ1 模剩余向量将问题 ( 1 ) (1) (1) 转化为先求一个约束不可微最优化问题，再解一个相容线性方程组，算法具有不扩大规模、简单、易于 *** 作的优点。

姚健康. 超定线性方程组极小 ℓ1 范数解的一个算法[J/OL]. 江西科学, 2007
(1): 1-3. http://d.g.wanfangdata.com.cn/Periodical_jxkx200701002.aspx. ↩︎
王嘉松，权日光. 不相容线性方程组极小极大解的一种新算法[J]. 1996: 04. ↩︎

欢迎分享，转载请注明来源：内存溢出

原文地址: http://outofmemory.cn/langs/570967.html

超定线性方程组Ax=b极小L1范数求解——MATLABPython实现

发表评论

评论列表（0条）