文章目录

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

一、前言

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

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

整篇文章目前已总结作为论文发表，可下载：基于残差向量l1l_1l1范数最小化与基追踪的多元线性模型参数估计方法

二、问题重述

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

A=[a11a12...a1na21a22...a2n............am1am2...amn]∈Rm×n,b=[b1b2⋮bm]∈Rm×1A=\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≥2m>n\ge 2m>n≥2，求 x=[x1,x2,...,xn]T∈Rn×1x={[x_1,x_2,...,x_n]}^T\in R^{n\times 1}x=[x1,x2,...,xn]T∈Rn×1，使得

min⁡∥Ax−b∥1=min⁡∑i=1m∣∑j=1naijxj−bi∣(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)

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

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

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

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

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

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

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

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

[IA2A1−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)

若我们令

A2A1−1=[c11c12...c1nc21c22...c2n............cm−n,1cm−n,2...cm−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 模剩余向量 rrr 满足如下线性方程组：

{rn+1−(c11r1+c12r2+⋯+c1nrn)=d1rn+2−(c21r1+c22r2+⋯+c2nrn)=d2⋯⋯rm−(cm−n,1r1+cm−n,2r2+⋯+cm−n,nrn)=dm−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)

其中

A2A1−1[b1b2⋮bn]−[bn+1bn+2⋮bm]=[d1d2⋮dm−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−1b1b2⋮bn−bn+1bn+2⋮bm=d1d2⋮dm−n(6)

具体转换过程请参考姚老师的论文。于是要求极小 ℓ1\ell_1ℓ1 模剩余向量 r{r}r，即求满足式 (5)(5)(5) 的 {r1,r2,...,rm}\{r_1, r_2, ..., r_m\}{r1,r2,...,rm}，使得 ∑i=1m∣ri∣\sum_{i=1}^{m}|r_i|∑i=1m∣ri∣ 最小。

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

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

{min∑i=1m∣zi∣s.t.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)\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∗=[z1∗,z2∗,...,zm∗]T{z}^*={[z_1^*, z_2^*, ..., z_m^*]}^Tz∗=[z1∗,z2∗,...,zm∗]T，方程组 Ax=b+z∗{A}x={b}+{z}^*Ax=b+z∗ 为原方程的一个相容线性方程组，其解 x∗{x}^*x∗ 即为问题 (1)(1)(1) 的解：

x∗=A1−1[b1b2⋮bn]+A1−1[z1∗z2∗⋮zn∗](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−1b1b2⋮bn+A1−1z1∗z2∗⋮zn∗(7)

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

C=[−c11−c12⋯−c1n10⋯0−c21−c22⋯−c2n01⋯0⋯⋯⋯⋯⋯⋯⋯⋯−cm−n,1−cm−n,2⋯−cm−n,n00⋯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) 变成了如下问题：

min∑i=1m∣zi∣,s.t.Cz=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 opdef OLE_linprog(A: list, b: list):'''Ax = b (Given A & b, try to derive x)Parameters----------A : [[a11,a12,...,a1n],[a21,a22,...,a1n],...,[am1,am2,...,amn]].b : [[b1],[b2],...,[bm]].Returns-------x : arrayMinimal L1 norm solution of the system of equations.Reference---------冯志强,张鸿燕.基于残差向量l_1范数最小化与基追踪的多元线性模型参数估计方法[J].海南师范大学学报(自然科学版),2022,35(03):250-259+267. (Available at:http://hsxblk.hainnu.edu.cn/ch/reader/view_abstract.aspx?file_no=20220303)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 A1A1_ = np.linalg.inv(A1)else:# Generalized inverse of A1A1_ = np.linalg.pinv(A1)# c_ij = A2 * A1_c = np.dot(A2, A1_)# r(n+1:m) = A2*inv(A1)*r(1:n) + dd = 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 > 0bounds = [(0, None) for i in range(2 * m)]# r0 = [u; v]r0 = op.linprog(t, A_eq=Aeq, b_eq=d, bounds=bounds,method='highs')['x']# Minimal L1-norm residual vector, r = u - vr = np.array([r0[:m] - r0[m:]])# Solving compatible linear equations Ax = b + r# Generalized inverse solutionx = np.linalg.pinv(A).dot(b + r.T)return x

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

六、检验与测试

1. 矩阵A满秩情况

对于线性方程 Ax=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⁡∥Ax−b∥1\min\ \left\|{A}x-{b}\right\|_1min ∥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<rank(A1)<mn < rank({A_1}) < mn<rank(A1)<m，我们的方法是不需要做对应改变的，而若矩阵 A{A}A 列缺秩，则需要使用矩阵 A1{A_1}A1 的Moore-Penrose逆矩阵 A1†{A_1^\dagger}A1† 替代 A1−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. ↩︎
冯志强,张鸿燕.基于残差向量l_1范数最小化与基追踪的多元线性模型参数估计方法[J].海南师范大学学报(自然科学版),2022,35(03):250-259+267. ↩︎

超定线性方程组Ax=b极小L1范数求解——MATLAB/Python实现相关推荐

【代数之美】线性方程组Ax=0的求解方法
在3D视觉中,我们常常会遇到这样一个问题:求解线性方程组Ax=0Ax=0Ax=0,从矩阵映射的角度来说,所有解组成了矩阵AAA的零空间.一个典型的场景比如用八点法求解本质矩阵EEE,参见我前面的博文: ...
MATLAB当中线性方程组、不定方程组、奇异方程组、超定方程组的介绍
系列文章目录 MATLAB绘图函数的相关介绍--海底测量.二维与三维图形绘制 MATLAB求函数极限的简单介绍文章目录一.线性方程组 1.1.线性方程组简介 1.2.矩阵的初等变换 1.3.MAT ...
超定方程的求解、最小二乘解、Ax=0、Ax=b的解，求解齐次方程组，求解非齐次方程组（推导十分详细）
本篇主要介绍的是超定方程组的求解,如果你不想看繁琐的推导过程,你可以直接看红字部分的结论! 1. 齐次线性方程组 Ax = 0 对于方程Ax=0\bm A \bm x = 0Ax=0,在我们实际的使用 ...
求超定方程组最小二乘解的三种方法
目录 1.超定线性方程组与最小二乘解 2.求解超定方程组的三种方法 3.参考链接 1.超定线性方程组与最小二乘解超定线性方程组:方程的个数大于解个数,方程组是无解的,但是我们可以求得其最小二 ...
机器学习中的正则化——L1范数和L2范数
机器学习中的正则化--L1范数和L2范数正则化是什么?为什么要正则化? LP范数 L0范数(补充了解) L1范数 L2范数 L1范数和L2范数的区别以深度学习的角度看待L1范数和L2范数正则化是 ...
【图像重建】基于matlab L1范数自适应双边总变分超分辨率图像序列重建【含Matlab源码 2209期】
一.正则化图像超分辨重建简介 1 超分辨率重建数学模型设有N帧低分辨率观测图像yk(k=1,2,-,N),图像大小为M×M,将每帧低分辨率(LR)图像yk按列方向排成向量的形式,记作Yk,大小为[M ...
matlab中欠定方程组超定方程组_一篇文章入门大规模线性方程组求解
前面介绍过主要的线性方程组求解库,参考附录.求解大规模线性方程组是仿真软件求解器的底层技术,求解器时间基本都消耗在方程组求解上.线性方程组的解法比较成熟,方法也有很多,而且不同的方法对应不同类型方程组 ...
适定、超定和欠定方程及压缩传感技术
原文"适定.超定和欠定方程",链接http://blog.sina.com.cn/s/blog_4b700c4c0102e5v3.html. 不定方程 http://wenku.b ...
线性方程组AX=b，AX=0以及非线性方程组的最小二乘解（解方程组-＞优化问题）
一.非齐次线性方程组AX=b的最小二乘解超定方程组无解是因为方程组包含了过多的约束条件,无法满足所有的约束条件,在这种情况下,方程组的某些方程必然是矛盾的,也就是说,他们描述的条件是不兼容的,无法同 ...

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

文章目录

一、前言

二、问题重述

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

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

四、算法伪码

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

六、检验与测试

1. 矩阵A满秩情况

2. 矩阵A列缺秩情况

七、总结

超定线性方程组Ax=b极小L1范数求解——MATLAB/Python实现相关推荐

最新文章

热门文章

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

文章目录

一、前言

二、问题重述

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

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

四、算法伪码

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

六、检验与测试

1. 矩阵A满秩情况

2. 矩阵A列缺秩情况

七、总结

超定线性方程组Ax=b极小L1范数求解——MATLAB/Python实现相关推荐

最新文章

热门文章

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