失眠网 > 最小二乘法加权最小二乘法迭代重加权最小二乘法 (含代码)【最小二乘线性求解】

最小二乘法加权最小二乘法迭代重加权最小二乘法 (含代码)【最小二乘线性求解】

时间：2023-06-19 20:16:59

文章目录

一：最小二乘法（OLS）1：概述2：代数式3：矩阵式（推荐）3.1：实现代码二：加权最小二乘法（WLS）1：增加对角矩阵 W1.1：实现代码三：迭代重加权最小二乘法（IRLS）四：应用1:拟合圆（算法：迭代重新加权最小二乘）2: 直线拟合（算法：迭代重新加权最小二乘）3:曲线拟合（算法：最小二乘）4：N点标定（包括9点标定）（算法：最小二乘）五：总结

原文链接：/Qvvjv

个人笔记：

最小二乘法，加权最小二乘法，迭代重加权最小二乘法。结合自己需要实现功能的目的，下面主要给出推导结果、代码实现和实际一些应用。推导过程最后会放一些个人参考的一些文章和资料。

这里推荐一个视频推导过程：利用矩阵求偏导得出 xxx=(ATA)−1ATB(A^TA)^{-1}A^TB(ATA)−1ATB矩阵乘法求导视频

一：最小二乘法（OLS）

1：概述

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

这里说个例子，比如已经确定了目标函数y=a0+a1x+a2x2y=a_0+a_1x+a_2x^2y=a0+a1x+a2x2，x和yx和yx和y是确定的实际数值，xxx是自变量，yyy是因变量，需要求a0,a1,a2a_0,a_1,a_2a0,a1,a2三个未知参数，这时候一般需要三个方程来组成一个方程组来求解得出这三个未知参数确定唯一的解。而实际中我们求三个未知参数一般都是在一个超定方程组（方程个数大于未知参数）中，这时候就需要用到最小二乘法来解决这个问题了，来求最优的解。下面给出代数式和矩阵式的解法。这里推荐使用矩阵式求解（非常方便）。

2：代数式

最小二乘法其思想主要是通过将理论值与预测值的距离的平方和达到最小。

举例：曲线拟合中最基本和最常用的是直线拟合。设xxx和yyy之间的函数关系为：一元一次函数y＝f(a0,a1)＝a0+a1xy＝f(a_0,a_1)＝a_0+a_1xy＝f(a0,a1)＝a0+a1x 代数推导：

分别对a0和a1a_0和a_1a0和a1求偏导，这里a0和a1a_0和a_1a0和a1就是未知参数

整理为方程组

然后化简得：

3：矩阵式（推荐）

举例：曲线拟合中最基本和最常用的是直线拟合。设xxx和yyy之间的函数关系为：一元一次函数y＝f(a0,a1)＝a0+a1xy＝f(a_0,a_1)＝a_0+a_1xy＝f(a0,a1)＝a0+a1x 使用矩阵表达：Ax=BAx = BAx=B，求xxx向量参数

推导过程我会在最后放上链接。那么：

如果是一元多项式函数：

其中m代表多项式的阶数，离散点与该多项式的平方和为F(a0,a1,...,am)F(a_0,a_1,..., a_m )F(a0,a1,...,am)。其中nnn代表采样点数：

一元多项式矩阵表达和一元一次项矩阵表达是一样的：Ax=BAx = BAx=B

3.1：实现代码

/* 普通最小二乘 Ax = B* (A^T * A) * x = A^T * B* x = (A^T * A)^-1 * A^T * B*/Array<double,Dynamic,1> GlobleFunction::leastSquares(Matrix<double,Dynamic,Dynamic> A, Matrix<double,Dynamic,1> B){//获取矩阵的行数和列数int rows = A.rows();int col = A.cols();//A的转置矩阵Matrix<double,Dynamic,Dynamic> AT;AT.resize(col,rows);//x矩阵Array<double,Dynamic,1> x;x.resize(col,1);//转置 ATAT = A.transpose();//x = (A^T * A)^-1 * A^T * Bx = ((AT * A).inverse()) * (AT * B);return x;}

Matrix是Eigen库中的一个矩阵类，这里引入Eigen库方便代数运算。

二：加权最小二乘法（WLS）

加权最小二乘法是对原模型进行加权，使之成为一个新的不存在异方差性的模型，然后采用普通最小二乘法估计其参数的一种数学优化技术。百度百科

1：增加对角矩阵 W

在最小二乘法的基础上增加一个对角矩阵WWW，对每组数据赋予不同的权重。

WT∗WW^T * WWT∗W对角矩阵每个数据的平方,消除负数。

1.1：实现代码

/* 加权最小二乘（WLS） W为对角线矩阵* W²(Ax - B) = 0* W²Ax = W²B* (A^T * W^T * W * A) * x = A^T * W^T * W * B* x = (A^T * W^T * W * A)^-1 * A^T * W^T * W * B*/Array<double,Dynamic,1> GlobleFunction::reweightedLeastSquares(Matrix<double,Dynamic,Dynamic> A, Matrix<double,Dynamic,1> B,Array<double,Dynamic,1> vectorW){//获取矩阵的行数和列数int rows = A.rows();int col = A.cols();//vectorW为空,默认构建对角线矩阵1if(vectorW.isZero()){vectorW.resize(rows,1);for(int i=0;i<rows;++i){vectorW(i,0) = 1;}}//A的转置矩阵Matrix<double,Dynamic,Dynamic> AT;AT.resize(col,rows);//x矩阵Array<double,Dynamic,1> x;x.resize(col,1);//W的转置矩阵Matrix<double,Dynamic,Dynamic> WT,W;W.resize(rows,rows);WT.resize(rows,rows);//生成对角线矩阵W = vectorW.matrix().asDiagonal();//转置WT = W.transpose();//转置 ATAT = A.transpose();// x = (A^T * W^T * W * A)^-1 * A^T * W^T * W * Bx = ((AT * WT * W * A).inverse()) * (AT * WT * W * B);return x;}

Matrix是Eigen库中的一个矩阵类，这里引入Eigen库方便代数运算。

三：迭代重加权最小二乘法（IRLS）

1：迭代重新加权最小二乘法（也叫迭代加权最小二乘法）( IRLS ) 的方法用于解决具有ppp范数形式的目标函数的某些优化问题：维基百科。迭代加权可以对目标函数和已知的数据来进行拟合，但是这时候有的数据对总体目标函数离得特别远，参与最小二乘法的时候就会对估计参数有很大的影响，这时候就要对参数进行优化，对数据较远的赋予较小的权重（不让它显得很重要，也就影响比较小），对数据较近的赋予较大的权重（影响大）。迭代加权最小二乘，是建立加权最小二乘上进行一个迭代来估值达到最优。

通过迭代方法，其中每一步都涉及解决以下形式的加权最小二乘问题：

2：下面引入一片论文中的一段迭代方法求解论文地址：Burrus, C.S. (). Iterative Reweighted Least Squares ∗.

MATLAB代码1：

% m-file IRLS0.m to find the optimal solution to Ax=b% minimizing the L_p norm ||Ax-b||_p, using basic IRLS.% csb 11/10/function x = IRLS0(A,b,p,KK)if nargin < 4, KK=10; end;x = pinv(A)*b; % Initial L_2 solutionE = [];for k = 1:KK % Iteratee = A*x - b; % Error vectorw = abs(e).^((p-2)/2); % Error weights for IRLSW = diag(w/sum(w)); % Normalize weight matrixWA = W*A; % apply weightsx = (WA'*WA)\(WA'*W)*b; % weighted L_2 sol.ee = norm(e,p); E = [E ee]; % Error at each iterationendplot(E)

MATLAB代码2：

% m-file IRLS1.m to find the optimal solution to Ax=b% minimizing the L_p norm ||Ax-b||_p, using IRLS.% Newton iterative update of solution, x, for M > N.% For 2<p<infty, use homotopy parameter K = 1.01 to 2% For 0<p<2, use K = approx 0.7 - 0.9% csb 10/20/function x = IRLS1(A,b,p,K,KK)if nargin < 5, KK=10; end;if nargin < 4, K = 1.5; end;if nargin < 3, p = 10; end;pk = 2; % Initial homotopy valuex = pinv(A)*b; % Initial L_2 solutionE = [];for k = 1:KK % Iterateif p >= 2, pk = min([p, K*pk]); % Homotopy change of pelse pk = max([p, K*pk]); ende = A*x - b; % Error vectorw = abs(e).^((pk-2)/2); % Error weights for IRLSW = diag(w/sum(w)); % Normalize weight matrixWA = W*A; % apply weightsx1 = (WA'*WA)\(WA'*W)*b; % weighted L_2 sol.q = 1/(pk-1); % Newton's parameterif p > 2, x = q*x1 + (1-q)*x; nn=p; % partial update for p>2else x = x1; nn=2; end % no partial update for p<2ee = norm(e,nn); E = [E ee]; % Error at each iterationendplot(E)

C++代码：

/* 迭代重加权最小二乘（IRLS） W为权重,p为范数* e = Ax - B* W = e^(p−2)/2* W²(Ax - B) = 0* W²Ax = W²B* (A^T * W^T * W * A) * x = A^T * W^T * W * B* x = (A^T * W^T * W * A)^-1 * A^T * W^T * W * B* 参考论文地址：/paper/Iterative-Reweighted-Least-Squares-%E2%88%97-Burrus/9b9218e7233f4d0b491e1582c893c9a099470a73*/Array<double,Dynamic,1> GlobleFunction::iterativeReweightedLeastSquares(Matrix<double,Dynamic,Dynamic> A, Matrix<double,Dynamic,1> B,double p,int kk){/* x(k) = q x1(k) + (1-q)x(k-1)* q = 1 / (p-1)*///获取矩阵的行数和列数int rows = A.rows();int col = A.cols();double pk = 2;//初始同伦值double K = 1.5;double epsilong = 10e-9; // εdouble delta = 10e-15; // δArray<double,Dynamic,1> x,_x,x1,e,w;x.resize(col,1);_x.resize(col,1);x1.resize(col,1);e.resize(rows,1);w.resize(rows,1);//初始x 对角矩阵w=1x = reweightedLeastSquares(A,B);//迭代最大迭代次数kkfor(int i=0;i<kk;++i){//保留前一个x值,用作最后比较确定收敛_x = x;if(p>=2){pk = qMin(p,K*pk);}else{pk = qMax(p,K*pk);}//偏差e = (A * x.matrix()) - B;//偏差的绝对值// 求矩阵绝对值：e = e.cwiseAbs(); 或 e.array().abs().matrix()e = e.abs();//对每个偏差值小于delta,用delta赋值给它for(int i=0;i<e.rows();++i){e(i,0) = qMax(delta,e(i,0));}//对每个偏差值进行幂操作w = e.pow(p/2.0-1);w = w / w.sum();x1 = reweightedLeastSquares(A,B,w);double q = 1 / (pk-1);if(p>2){x = x1*q + x*(1-q);}else{x = x1;}//达到精度,结束if((x-_x).abs().sum()<epsilong){return x;}}return x;}

C++实现代码和MATLAB基本一样，不过有稍微改进，其中参考了维基百科和Burrus, C.S. (). Iterative Reweighted Least Squares ∗.