最小二乘法，加权最小二乘法，迭代重加权最小二乘法

news2024/11/16 9:33:54

文章目录

    • 一：最小二乘法（OLS）
      • 1：概述
      • 2：代数式
      • 3：矩阵式（推荐）
        • 3.1：实现代码
    • 二：加权最小二乘法（WLS）
      • 1：增加对角矩阵 W
        • 1.1：实现代码
    • 三：迭代重加权最小二乘法（IRLS）
    • 四：应用
      • 1:拟合圆
      • 2: 曲线拟合
    • 五：总结

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

一：最小二乘法（OLS）

1：概述

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法是曲线拟合的常用方法，使用该方法对匹配函数的选取非常重要。。所谓匹配函数就是函数经过的路线在图中的点达到一个最佳匹配。不然就会出现过拟合和欠拟合的现象。
在这里插入图片描述在这里插入图片描述

在这里插入图片描述在这里插入图片描述

2：代数式

最小二乘法其思想主要是通过将理论值与预测值的距离的平方和达到最小。
在这里插入图片描述举例：曲线拟合中最基本和最常用的是直线拟合。设x和y之间的函数关系为：一元一次函数（F(a0,a1)＝a0+a1x）y＝a0+a1x 代数推导：
在这里插入图片描述

分别对a0和a1求偏导
在这里插入图片描述

整理为方程组
在这里插入图片描述

然后化简得：
在这里插入图片描述

3：矩阵式（推荐）

举例：曲线拟合中最基本和最常用的是直线拟合。设x和y之间的函数关系为：一元一次函数（F(a0,a1)＝a0+a1x）y＝a0+a1x 使用矩阵表达：Ax = B
在这里插入图片描述

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

在这里插入图片描述

如果是一元多项式函数：
在这里插入图片描述
其中m代表多项式的阶数，离散点与该多项式的平方和 为F(a0,a1,…, am ) 。其中n代表采样点数：

在这里插入图片描述

一元多项式矩阵表达和一元一次项矩阵表达是一样的：Ax = 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);

    //转置 AT
    AT = A.transpose();

    //x = (A^T * A)^-1 * A^T * B
    x = ((AT * A).inverse()) * (AT * B);
    return x;

}

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

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

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

1：增加对角矩阵 W

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

在这里插入图片描述
W^T * 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为空,默认构建对角线矩阵1
    if(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();
    //转置 AT
    AT = A.transpose();

    // x = (A^T * W^T * W * A)^-1 * A^T * W^T * W * B
    x = ((AT * WT * W * A).inverse()) * (AT * WT * W * B);
    return x;
}

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

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

1：迭代重新加权最小二乘法( IRLS ) 的方法用于解决具有p 范数形式的目标函数的某些优化问题：维基百科
在这里插入图片描述

通过迭代方法，其中每一步都涉及解决以下形式的加权最小二乘问题：
在这里插入图片描述在这里插入图片描述

2：下面引入一片论文中的一段迭代方法求解 论文地址：Burrus, C.S. (2014). 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/2012
function x = IRLS0(A,b,p,KK)
if nargin < 4, KK=10; end;
x = pinv(A)*b; 					% Initial L_2 solution
E = [];
for k = 1:KK 					% Iterate
	e = A*x - b; 				% Error vector
	w = abs(e).^((p-2)/2); 		% Error weights for IRLS
	W = diag(w/sum(w)); 		% Normalize weight matrix
	WA = W*A; 					% apply weights
	x = (WA'*WA)\(WA'*W)*b; 	% weighted L_2 sol.
	ee = norm(e,p); E = [E ee]; % Error at each iteration
end
plot(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/2012
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 value
x = pinv(A)*b; 						% Initial L_2 solution
E = [];
for k = 1:KK 						% Iterate
	if p >= 2, pk = min([p, K*pk]);	    % Homotopy change of p
	else pk = max([p, K*pk]); end
	e = A*x - b; 						% Error vector
	w = abs(e).^((pk-2)/2); 			% Error weights for IRLS
	W = diag(w/sum(w)); 				% Normalize weight matrix
	WA = W*A; 							% apply weights
	x1 = (WA'*WA)\(WA'*W)*b;		    % weighted L_2 sol.
	q = 1/(pk-1); 						% Newton's parameter
	if p > 2, x = q*x1 + (1-q)*x; nn=p; % partial update for p>2
	else x = x1; nn=2; end 				% no partial update for p<2
	ee = norm(e,nn); E = [E ee]; 		% Error at each iteration
end
plot(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
 * 参考论文地址：https://www.semanticscholar.org/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=1
    x = reweightedLeastSquares(A,B);

    //迭代  最大迭代次数kk
    for(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. (2014). Iterative Reweighted Least Squares ∗.

四：应用

以下针对超定方程组来求解的。数据个数大于未知数。

1:拟合圆

1：使用最小二乘效果，可以看出外部的噪点干扰还是比较大的。下面使用迭代重加权最小二乘算法优化。
在这里插入图片描述2：迭代重加权最小二乘
第1次迭代
在这里插入图片描述
第2次迭代
在这里插入图片描述第3次迭代
在这里插入图片描述第4次迭代
在这里插入图片描述…
第20次迭代
在这里插入图片描述

2: 曲线拟合

1：直线拟合 y = a0 + a1x
最小二乘：
在这里插入图片描述1.2：迭代重加权最小二乘
第1次迭代
在这里插入图片描述…
第100次迭代
在这里插入图片描述2：曲线拟合
一次函数拟合：
在这里插入图片描述
二次函数拟合：

在这里插入图片描述
三次函数拟合：
在这里插入图片描述
四次函数拟合：
在这里插入图片描述
五次函数拟合：

在这里插入图片描述
六次函数拟合：

在这里插入图片描述
七次函数拟合：
在这里插入图片描述
八次函数拟合：

在这里插入图片描述

九次函数拟合：
在这里插入图片描述
可以发现函数在第五次函数的时候拟合程度很好了，越往后越来过拟合了。

3：N点标定（包括9点标定）
9点标定是求视觉 中像素坐标和世界坐标建立的关系。
在这里插入图片描述可以看到 halcon算子块 vector_to_hom_mat2d 就是用最小二乘法来计算矩阵的。图中外部算法【2】其实是本文章中的最小二乘算法实现的。内部算法【1】实现是用求偏导计算的，在这篇文章有实现N点标定

五：总结

1：工具：主要Qt + Eigen库 + QCustomPlot类
Eigen库是一个用于矩阵计算，代数计算库
QCustomPlot类是一个用于绘图和数据可视化

2：上面完整代码已上传GitHub

3：其他链接
最小二乘代数推导
最小二乘矩阵推导
最小二乘法？为神马不是差的绝对值
正则化
鲁棒学习算法
最小二乘法的原理理解

最后：舞台再大，你不上台，永远是个观众。 平台再好，你不参与，永远是局外人。能力再大，你不行动，只能看别人成功！没有人会关心你付出过多少努力，撑得累不累，摔得痛不痛，他们只会看你最后站在什么置，然后美慕或鄙夷。

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/12228.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

oracle数据库控制语言—DCL

oracle数据库控制语言—DCL

文章目录1、授予系统权限1.1 授予创建其他对象权限2、撤销系统权限2.1 示例3、oracle 中的角色3.1 什么时角色3.2 创建角色并且授予权限给角色3.2.1 创建角色3.2.1.1 示例3.2.2 授予权限给一个角色3.2.2.1 示例3.2.3 授予一个角色给用户3.2.3.1 示例一3.2.3.2 示例二3.2.3.3 示…
阅读更多...
基于MySQL的事务管理

基于MySQL的事务管理

目录 概念&#xff1a;多条语句组成一个执行单位 事务的基本操作 MySQL中的事务必须满足A,C,I,D这四个基本特性 事务操作举例——&#xff08;转账&#xff09; 事务保存点——SAVEPOINT 事务隔离级别——多线程(并发同时访问) 总结 概念&#xff1a;多条语句组成一个执…
阅读更多...
Mongo的数据操作

Mongo的数据操作

文章目录一&#xff0c;创建数据库二&#xff0c;插入数据&#xff08;一&#xff09;插入单条数据1&#xff0c;insert2&#xff0c;save&#xff08;二&#xff09;插入多条数据三&#xff0c;修改数据四&#xff0c; 更新所有找到匹配的数据五&#xff0c;数据删除&#xff…
阅读更多...
极速Go语言入门(超全超详细)-基础篇

极速Go语言入门(超全超详细)-基础篇

文章目录 GoLang概述 Go语言三大牛谷歌创造Golang的原因Golang 的发展历程Golang 的语言的特点 Go语言开发工具Go开发环境配置(sdk下载及配置) 使用开发工具创建第一个Go项目 Go 程序开发的注意事项 官方参考文档 Go学习 Go变量 数据类型 标识符 运算符 键盘输入语句 程序流程…
阅读更多...
[附源码]计算机毕业设计JAVA高校贫困生认定系统

[附源码]计算机毕业设计JAVA高校贫困生认定系统

[附源码]计算机毕业设计JAVA高校贫困生认定系统 项目运行 环境配置&#xff1a; Jdk1.8 Tomcat7.0 Mysql HBuilderX&#xff08;Webstorm也行&#xff09; Eclispe&#xff08;IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持&#xff09;。 项目技术&#xff1a; SSM myba…
阅读更多...
列表的嵌套--Python

列表的嵌套--Python

#列表的嵌套&#xff1a;在每一个列表中都包含其他列表元素 #获取李四所在列表的值? #获取李四所在的子列表索引&#xff0c;并通过索引获取改子列表值 name_list [[小明,小红,小绿],[Tom,Lily,Rose],[张三,李四,王五]] print(name_list[2]) #在从子列表中通过李四所在的索引…
阅读更多...
[翻译] 使用FXGL创建一个非常基本的游戏

[翻译] 使用FXGL创建一个非常基本的游戏

游戏要求 首先&#xff0c;让我们为我们的简单游戏定义一些要求: 一个600x600的窗口。屏幕上的玩家&#xff0c;由蓝色矩形表示。可以通过按键盘上的W、S、A或D来移动玩家。UI由一行文本表示。当玩家移动时&#xff0c;UI文本会更新以显示玩家在其生命周期内移动了多少像素。 …
阅读更多...
今天给在家介绍一篇基于jsp的旅游网站设计与实现

今天给在家介绍一篇基于jsp的旅游网站设计与实现

项目描述 临近学期结束&#xff0c;还是毕业设计&#xff0c;你还在做java程序网络编程&#xff0c;期末作业&#xff0c;老师的作业要求觉得大了吗?不知道毕业设计该怎么办?网页功能的数量是否太多?没有合适的类型或系统?等等。这里根据疫情当下&#xff0c;你想解决的问…
阅读更多...
zookeeper报错length is greater than jute.maxbuffer=1048575

zookeeper报错length is greater than jute.maxbuffer=1048575

1、场景 最近在给上云项目部署系统&#xff0c;通过压测都已经正式上生产后发现kafka存在异常错误&#xff0c;经排查发现zookeeper也存在错误&#xff0c;怀疑kafka的问题可能是由于zk异常到的&#xff0c;报错如下 2022-11-17 06:26:43,052 [myid:] - WARN [NIOWorkerThr…
阅读更多...
HTML学生个人网站作业设计：游戏网站设计——原神首页 1页 带轮播图

HTML学生个人网站作业设计：游戏网站设计——原神首页 1页 带轮播图

⛵ 源码获取 文末联系 ✈ Web前端开发技术 描述 网页设计题材&#xff0c;DIVCSS 布局制作,HTMLCSS网页设计期末课程大作业 | 游戏官网 | 游戏网站 | 电竞游戏 | 游戏介绍 | 等网站的设计与制作 | HTML期末大学生网页设计作业&#xff0c;Web大学生网页 HTML&#xff1a;结构 …
阅读更多...
Burpsuite简介及MIME上传绕过实例

Burpsuite简介及MIME上传绕过实例

目录预备知识1.了解burpsuite2.了解服务端MIME类型检测实验目的实验环境实验步骤一使用Burpsuite的代理功能Target模块实验步骤二使用burpsuite上传绕过服务端MIME类型检测预备知识 1.了解burpsuite Burp Suite是用于攻击web应用程序的集成平台。包含了许多工具&#xff0c;并…
阅读更多...
[附源码]Python计算机毕业设计 学生宿舍管理系统

[附源码]Python计算机毕业设计 学生宿舍管理系统

项目运行 环境配置&#xff1a; Pychram社区版 python3.7.7 Mysql5.7 HBuilderXlist pipNavicat11Djangonodejs。 项目技术&#xff1a; django python Vue 等等组成&#xff0c;B/S模式 pychram管理等等。 环境需要 1.运行环境&#xff1a;最好是python3.7.7&#xff0c;…
阅读更多...
我们与元宇宙的距离

我们与元宇宙的距离

元宇宙的发展是一个循序渐进的过程&#xff0c;需要经过初始阶段、规划阶段、系统阶段和优化阶段。当前&#xff0c;虽然已经有许多元宇宙完成了搭建&#xff0c;但都是针对单一特定场景的模拟仿真&#xff0c;并且没有实现全面的推广应用&#xff0c;因此我们仍处于元宇宙的初…
阅读更多...
浅析多通道接收单元噪声系数的测试

浅析多通道接收单元噪声系数的测试

之前一个朋友要测试低噪声放大器(LNA)的噪声系数&#xff0c;但是声称遇到一些麻烦。LNA噪声系数的测试采用Y因子法非常简便&#xff0c;校准完成后直接连接待测件即可测试&#xff0c;可操作性非常强。麻烦在哪里呢&#xff1f; 原来待测件是一个含有四个通道的接收模块&…
阅读更多...
2023年三大网络安全威胁不容忽视

2023年三大网络安全威胁不容忽视

2022年已进入尾声&#xff0c;降低数字化风险、增强安全防御能力依然是众多企业组织数字化发展中的重要需求和目标。随着技术的不断进步&#xff0c;网络攻击者的攻击成本不断降低&#xff0c;同时攻击方式更加先进&#xff0c;美国《福布斯》网站在近日的报道中列出了2023年值…
阅读更多...
Zoho 如何使用低代码 #1 - 赋予人力资源以技术实力

Zoho 如何使用低代码 #1 - 赋予人力资源以技术实力

Zoho 为客户提供了一套跨功能产品&#xff0c;从运行简单的调查到简化复杂的企业组织职能&#xff0c;Zoho 几乎提供了企业的业务运行所需的一切。 组织在新的规范和挑战中不断进行扩展&#xff0c;这就不断需要构建可定制的解决方案。这就是为什么除了现成的应用程序之外&…
阅读更多...
第1关：ZooKeeper初体验

第1关：ZooKeeper初体验

ZooKeeper安装方法 由于本实验环境已经安装ZooKeeper并配置&#xff0c;下面主要讲述一般环境的安装方法。 可以从ZooKeeper的官方网站上下载稳定版&#xff0c;下载地址如下&#xff1a;Apache ZooKeeper 下载后&#xff0c;利用tar命令将压缩包解压到/opt/zookeeper-3.4.1…
阅读更多...
实验31：温湿度传感器实验

实验31：温湿度传感器实验

本实验也比较简单 用LCD1602显示温湿度传感器返回的温湿度值 本专栏就要结尾了 希望对大家有一定的帮助 学习就是一点一滴的 01 硬件电路设计 整体电路图 重点还是接口： 因为最后相放一个学生做的大创为结尾彩蛋，所以这里就不再添加报警灯什么的额外操作了 每一个实验…
阅读更多...
OVS 和 OVS-DPDK 对比

OVS 和 OVS-DPDK 对比

OVS 目前有两种比较突出的架构&#xff0c;一种是原生的 OVS 架构&#xff08;使用 kernel 作为 datapath&#xff09;&#xff0c;一种是基于 DPDK 的架构&#xff08;使用用户空间作为 datapath&#xff09;。 原生 OVS 原生 OVS 架构如下所示&#xff0c;主要包含两个组件…
阅读更多...
冷启动问题分析与解决办法

冷启动问题分析与解决办法

1、什么是冷启动问题&#xff1f; 在缺乏有价值数据的时候&#xff0c;如何有效地满足业务需求的问题&#xff0c;就是“冷启动问题”。为了沟通方便&#xff0c;下面统一从推荐系统的角度来讲“冷启动问题”&#xff0c;其他业务场景同理。 冷启动问题是机器学习系统中十分常…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023