前言

终于赶完期末大作业啦，抓不住我（不是。上次说到用LU分解的方式求解线性方程组，但是有时候当方程组不相容时常规方法就会失效，这时候我们可以使用最小二乘法，当然对于矩阵形式的拟合，我们可以利用后面介绍的M-P广义逆直接计算其最小二乘逼近解，开搞！

正文开始~ CSDN (゜-゜)つロ干杯

知识点

完全懵x的童鞋可以翻看相关的线性代数或者矩阵论知识~

对于常规意义下的矩阵逆我们有AB=BA=IAB=BA=IAB=BA=I形式，其中A为n阶可逆矩阵，B为其逆矩阵，记为A−1=BA^{-1}=BA−1=B。
若矩阵A满足条件AHA=AAHA^{H}A=AA^{H}AHA=AAH，则成A是一个正规矩阵，且有rank(A)=rank(AHA)=rank(AAH)rank(A)=rank(A^{H}A)=rank(AA^{H})rank(A)=rank(AHA)=rank(AAH)。正规矩阵的充分必要条件为A有n个线性无关的特征向量，且它们构成空间CnC^nCn的标准正交基。如果此时AHAA^{H}AAHA是正定矩阵，即rank(AHA)=nrank(A^{H}A)=nrank(AHA)=n，则AHAA^{H}AAHA可逆，对应于n个非零的特征值。
若存在A∈Cm×n,B∈Cn×mA \in C^{m \times n}, B \in C^{n \times m}A∈Cm×n,B∈Cn×m，使得BA=InBA=I_nBA=In。则称A是左可逆的，B为A的一个左逆矩阵，记为AL−1A_L^{-1}AL−1。在《矩阵论》（杨明，刘先忠）教材中有证明得到对于左可逆矩阵A，其有列满秩的性质，即rank(A)=nrank(A)=nrank(A)=n，则由知识点2可知，AHAA^{H}AAHA可逆，此时A存在一个左逆矩阵(AHA)−1AH(A^{H}A)^{-1}A^{H}(AHA)−1AH。
同样的，若存在C∈Cn×mC \in C^{n \times m}C∈Cn×m，使得AC=ImAC=I_mAC=Im，则成C是A的一个右逆矩阵，记为AR−1A_R^{-1}AR−1。此时可以证明得到A是行满秩的，A存在一个右逆矩阵AH(AAH)−1A^{H}(AA^{H})^{-1}AH(AAH)−1。
在左逆右逆的基础上，提出减号广义逆，有AGA=AAGA=AAGA=A，可以证明左逆右逆自动满足减号逆的定义。
在减号逆的基础上，Moore和Penrose提出广义加号逆，满足：
（1）AGA=AAGA=AAGA=A
（2）GAG=GGAG=GGAG=G
（3）(AG)H=AG(AG)^H=AG(AG)H=AG
（4）(GA)H=GA(GA)^H=GA(GA)H=GA
记广义加号逆G为A+A^+A+。

若矩阵A是列满秩的，则其加号逆为(AHA)−1AH(A^{H}A)^{-1}A^{H}(AHA)−1AH；
若矩阵A是行满秩的，则其加号逆为AH(AAH)−1A^{H}(AA^{H})^{-1}AH(AAH)−1。
结合前面我们提到的满秩分解，可证明得到以下结论，任意矩阵A∈Cm×nA \in C^{m \times n}A∈Cm×n都存在广义加号逆，设rank(A)=rrank(A)=rrank(A)=r，A的一个满秩分解为A=BC，B∈Cm×r,C∈Cr×mA=BC，B \in C^{m \times r}, C \in C^{r \times m}A=BC，B∈Cm×r,C∈Cr×m，则有
A+=CR−1BL−1=CH(CCH)−1(BHB)−1BHA^{+} = C_R^{-1}B_L^{-1}=C^{H}(CC^H)^{-1}(B^HB)^{-1}B^HA+=CR−1BL−1=CH(CCH)−1(BHB)−1BH。
借助正交投影变换（有兴趣的uu可查阅相关资料，这里不展开），我们可以证明得到对于不相容的方程组Ax=bAx=bAx=b，x0=A+bx_0=A^+bx0=A+b是其的最佳最小二乘解，最佳在这里的意义是该最小二乘解是在范数约束下的最佳逼近。

初等行变换求逆

由上面的内容可以知道为了得到不相容矩阵的最佳最小二乘解，我们需要求A的广义加号逆，那就要求A的满秩分解（其实也可以用奇异值分解，但这个坑我还没填上hhh），再对满秩分解的结果进行左右逆计算，所以最本质的操作是怎么计算矩阵的逆（这里我们可以保证所求矩阵AAHAA^{H}AAH或AHAA^{H}AAHA都是可逆的）。最最普遍的方法，从大一说到研一，那当然是初等行变换求逆啦。大致思路（还是啰嗦一句），将A矩阵与单位矩阵构成一个大矩阵[A∣I][A|I][A∣I]，通过初等行变换将A化为单位阵的同时，右边的单位阵在相同操作下得到的矩阵即为A的逆矩阵。

代码实现

@classmethod
def Primary_Row_Transformation_Inv(self, target, eps=1e-6, test=False):"""Primary_Row_Transformation_Inv Args:target ([np.darray]): [target matrix]eps ([float]): numerical threshold.test (bool, optional): [whether to show testing information]. Defaults to False.Returns:[np.darray]: [A-]Last edited: 22-01-09Author: Junno"""assert len(np.shape(target)) == 2# get row and col numbersM, N = np.shape(target)assert M == N, 'This matrix is not square-like, so it can\'t use PRT to get inverse and you can use MP-inverse to calculate its general inv'assert self.Get_row_rank(target, eps) == M, 'This matrix is not full-rank, so it can\'t use PRT to get inverse and you can use MP-inverse to calculate its general inv'E = np.eye(M)A = deepcopy(target)if test:print('original matrix:')print(target)for i in range(M):for j in range(i, N):if test:print('During process in row:{}, col:{}'.format(i, j))if sum(abs(A[i:, j])) > eps:# sort rows by non-zero values from small to largezero_index = []non_zero_index = []for k in range(i, M):if abs(A[k, j]) > eps:non_zero_index.append(k)else:zero_index.append(k)non_zero_index = sorted(non_zero_index, key=lambda x: abs(A[x, j]))sort_cols = non_zero_index+zero_indexif test:print("Sort_cols index: {}".format(sort_cols))# E should be processed with the same operations.A[i:, :] = A[sort_cols, :]E[i:, :] = E[sort_cols, :]if test:print('After resorting')print(A)print(E)# eliminate elements belowsprefix = -A[i+1:, j]/A[i, j]temp = (np.array(prefix).reshape(-1, 1))@A[i, :].reshape((1, -1))A[i+1:, :] += temptemp = (np.array(prefix).reshape(-1, 1))@E[i, :].reshape((1, -1))E[i+1:, :] += tempif test:print('After primary row transformation:')print(A)print(E)# eliminate elements abovefor k in range(i):if abs(A[k, j]) > eps:# be careful to using intermediate parametertemp = -(A[k, j]/A[i, j])A[k, :] += temp*A[i, :]E[k, :] += temp*E[i, :]if test:print('After eliminating elements above:')print(A)print(E)breakelse:continue# normalize the first element of each row to 1for m in range(M):if abs(A[m, m]) > eps:if A[m, m] != 1:temp = 1./A[m, m]A[m, :] *= tempE[m, :] *= tempreturn E

test example

>>> import numpy as np
>>> from Matrix_Solutions import Matrix_Solutions
>>> A=np.array([[ 1., -3.,  7.],
...        [ 2.,  4., -3.],
...        [-3.,  7.,  2.]], dtype=np.float32)>>> inv_A=Matrix_Solutions.Primary_Row_Transformation_Inv(A)
>>> inv_A
array([[ 0.148,  0.281, -0.097],[ 0.026,  0.117,  0.087],[ 0.133,  0.01 ,  0.051]])>>> A@inv_A
array([[ 1., -0., -0.],[-0.,  1., -0.],[ 0.,  0.,  1.]])
>>> inv_A@A
array([[ 1.,  0., -0.],[ 0.,  1.,  0.],[-0.,  0.,  1.]])

有了求逆工具后，我们的各种广义逆也就水到渠成，直接上代码

@classmethod
def Left_inv(self, target):"""Get_Left_Inv Args:target ([np.darray]): [target matrix]Returns:[np.darray]: [A_L]Last edited: 22-01-09Author: Junno"""assert len(np.shape(target)) == 2# get row ans col numbersM, N = np.shape(target)assert self.Get_col_rank(target) == Nreturn self.Primary_Row_Transformation_Inv(target.T@target)@target.T@classmethod
def Right_inv(self, target):"""Get_Right_Inv Args:target ([np.darray]): [target matrix]Returns:[np.darray]: [A_R]Last edited: 22-01-09Author: Junno"""assert len(np.shape(target)) == 2# get row ans col numbersM, N = np.shape(target)assert self.Get_row_rank(target) == Mreturn target.T@self.Primary_Row_Transformation_Inv(target@target.T)@classmethod
def Moore_Penrose_General_Inv(self, target):"""Moore_Penrose_General_Inv Args:target ([np.darray]): [target matrix]Returns:[np.darray]: [A+]Last edited: 22-01-09Author: Junno"""assert len(np.shape(target)) == 2# get row ans col numbersM, N = np.shape(target)B, C = self.Full_Rank_Fact(target)return self.Right_inv(C)@self.Left_inv(B)

test example

>>> test_m=np.array([
...     [1,0,-1,1],
...     [0,2,2,2],
...     [-1,4,5,3]
... ]).reshape((3,4)).astype(np.float32)>>> test_m
array([[ 1.,  0., -1.,  1.],[ 0.,  2.,  2.,  2.],[-1.,  4.,  5.,  3.]], dtype=float32)>>> MPI=Matrix_Solutions.Moore_Penrose_General_Inv(test_m)
>>> MPI
array([[ 0.278,  0.111, -0.056], [ 0.056,  0.056,  0.056], [-0.222, -0.056,  0.111], [ 0.333,  0.167,  0.   ]])# AGA=A
>>> test_m@MPI@test_m
array([[ 1.,  0., -1.,  1.],[-0.,  2.,  2.,  2.],[-1.,  4.,  5.,  3.]])# GAG=G
>>> MPI@test_m@MPI
array([[ 0.278,  0.111, -0.056],[ 0.056,  0.056,  0.056],[-0.222, -0.056,  0.111],[ 0.333,  0.167,  0.   ]])# (AG)^T=AG, 即AG是对称矩阵
>>> test_m@MPI
array([[ 0.833,  0.333, -0.167],[ 0.333,  0.333,  0.333],[-0.167,  0.333,  0.833]])# (GA)^T=GA, 即GA是对称矩阵
>>> (MPI@test_m)
array([[ 0.333,  0.   , -0.333,  0.333],[ 0.   ,  0.333,  0.333,  0.333],[-0.333,  0.333,  0.667, -0.   ],[ 0.333,  0.333,  0.   ,  0.667]])

应用：拟合数据点曲线

最小二乘法最常见的应用莫过于拟合曲线，使曲线与各个观测点的距离最小，得到x和y的对应关系。给出教材中的一个例子：
设有一组实验数据，（1，2），（2，3），（3，5），（4，7），我们希望使用直线
y=β0+β1xy=\beta_0+\beta_1 xy=β0+β1x来拟合数据点，求最佳拟合直线。列出矩阵形式：
[11121314][β0β1]=[2357]\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \end{bmatrix} \begin{bmatrix} \beta_0\\ \beta_1 \end{bmatrix}= \begin{bmatrix} 2 \\ 3 \\ 5 \\ 7 \end{bmatrix} ⎣⎢⎢⎡11111234⎦⎥⎥⎤[β0β1]=⎣⎢⎢⎡2357⎦⎥⎥⎤
设系数矩阵为A，则我们只要求得A的广义逆，那么β\betaβ就可以表示为A+bA^+bA+b。
最后拟合的误差可以通过∣∣δ∣∣2=∣∣Aβ−b∣∣2||\delta||_2=||A\beta-b||_2∣∣δ∣∣2=∣∣Aβ−b∣∣2来计算。当然我们也可以用高阶来拟合，这样的误差会更小，这在下面的实现中也会体现到。

代码实现

@classmethod
def Match_Points_LSM(self, ps, order):"""Match Points using LSM and show resultsArgs:ps ([type]): [arrays of (x,y) points]order ([type]): [fitting order]Last edited: 21-12-31Author: Junno"""M, N = ps.shape[0], order+1A = np.zeros((M, N))# generate coefficient matrix A with given orderfor c in range(N):A[:, c] = np.power(ps[:, 0], c)# get MP-inv of AL_Inv = self.Moore_Penrose_General_Inv(A)# solve betabeta = L_Inv@ps[:, 1]# show resultsprint('The matching line: y=', sep='', end='')for i in range(order+1):print('{:-.3f}*x^{}'.format(beta[i], i), sep='', end='')print('\n')# calculate fitting errorerror = np.linalg.norm(A@beta-ps[:, 1])print("Norm of Matching Error: {:.5f}".format(error))# prepare to show matching linex = np.arange(np.min(ps[:, 0])-1,np.max(ps[:, 0])+1, 0.01).reshape((-1,))B = np.zeros((x.shape[0], N))for c in range(N):B[:, c] = np.power(x, c)y = B@betaplt.figure(0)fig, ax = plt.subplots()raw = ax.scatter(ps[:, 0], ps[:, 1], marker='o', c='blue', s=10)match = ax.scatter(x, y, marker='v', c='red', s=1)plt.legend([raw, match], ['raw', 'match'], loc='upper left')plt.title('Fitting points with {} order'.format(order))plt.xlabel('x')plt.ylabel('y')plt.show()

test example

>>> a=np.array([
...     [1,2],
...     [2,3],
...     [3,5],
...     [4,7]
... ]).reshape((-1,2)).astype(np.float32)# 一阶拟合
>>> Matrix_Solutions.Match_Points_LSM(a,1)
The matching line: y=+0.000*x^0+1.700*x^1
Norm of Matching Error: 0.54772# 二阶拟合
>>> Matrix_Solutions.Match_Points_LSM(a,2)
The matching line: y=+1.250*x^0+0.450*x^1+0.250*x^2
Norm of Matching Error: 0.22361# 三阶拟合
>>> Matrix_Solutions.Match_Points_LSM(a,3)
The matching line: y=+3.000*x^0-2.333*x^1+1.500*x^2-0.167*x^3
Norm of Matching Error: 0.00000

总结

今天的内容先到这，对代码有任何问题或者指正欢迎留言，我会第一时间进行回复~
不求对原理的证明有多么深刻，更多的是能掌握这类问题的解决工具。下回继续填坑！

曲高未必人不识，自有知音与清词。

从左逆右逆广义逆到求解线性方程组的最小二乘解相关推荐

MATLAB中：左右除法、逆inv、广义逆pinv的区别
以下是通过实验得出的一些结论: 左除行相等,右除列相等.只要满足此条件便可运算,且左.右除意义不相同! 逆inv()仅针对"非奇异方阵|A|≠0",使得A^(-1) A=A ...
线性代数学习笔记10-4：左右逆、伪逆/M-P广义逆（从四个子空间和SVD角度理解）
下面讨论m×nm\times nm×n的秩为rrr的矩阵对于不同情况,讨论逆矩阵两侧逆矩阵 2-sided inverse 这也是一般所说的"逆矩阵"的含义方阵A\bolds ...
线性高斯反问题--广义逆
0 解与算符对于线性反问题Gm=d\mathbf{Gm=d}Gm=d的方法,长度方法基于分析解的两个属性:预测误差和解的简单程度.而广义逆的方法,则希望通过研究算符矩阵来获得更多关于反问题的属性. ...
【矩阵论】4. 矩阵运算——广义逆——减号逆
矩阵论 1. 准备知识--复数域上矩阵,Hermite变换) 1.准备知识--复数域上的内积域正交阵 1.准备知识--Hermite阵,二次型,矩阵合同,正定阵,幂0阵,幂等阵,矩阵的秩 2. 矩阵分 ...
【矩阵论】4. 矩阵运算——广义逆——加号逆应用
矩阵论 1. 准备知识--复数域上矩阵,Hermite变换) 1.准备知识--复数域上的内积域正交阵 1.准备知识--Hermite阵,二次型,矩阵合同,正定阵,幂0阵,幂等阵,矩阵的秩 2. 矩阵分 ...
【矩阵论】4. 矩阵运算——广义逆——加号逆的计算
矩阵论 1. 准备知识--复数域上矩阵,Hermite变换) 1.准备知识--复数域上的内积域正交阵 1.准备知识--Hermite阵,二次型,矩阵合同,正定阵,幂0阵,幂等阵,矩阵的秩 2. 矩阵分 ...
矩阵分析与应用（一）——矩阵基础知识、广义逆
文章目录前言部分符号约定关于矩阵理论的碎碎念一些基础知识与本门课知识串讲矩阵奇异与线性无关向量空间.内积.范数行列式.特征值.迹逆.广义逆矩阵方程.向量化.Kronecker积向量 ...
Moore-Penrose广义逆：可解决MATLAB报错“矩阵接近奇异值，或者缩放错误。结果可能不准确”
上一篇博文讲到:<方程AX=b的解的讨论(特解.通解.零空间向量等概念)及其MATLAB实现>,程序中用到的是mldivide或者A\b的方法(二者相同)来解方程. 但实际上运行过程中我们 ...
左神小和问题逆序对问题面试
归并排序的扩展问题 --先了解数据结构之归并排序,更容易理解以下的问题. 小和问题和逆序对问题. 在一个数组中,每一个数左边比当前数小的数累加起来,叫做这个数组的小和.在一个数组中,左边的数如果比右边 ...
求非线性方程组的最小二乘解的广义逆法C实现
求非线性方程组的最小二乘解的广义逆法 #include "math.h"#include "stdlib.h"#include "6gmiv.c&qu ...

从左逆右逆广义逆到求解线性方程组的最小二乘解

矩阵的广义逆及应用

前言

知识点

初等行变换求逆

代码实现

test example

test example

应用：拟合数据点曲线

代码实现

test example

总结

从左逆右逆广义逆到求解线性方程组的最小二乘解相关推荐

最新文章

热门文章