精华内容
下载资源
问答
  • 最小二乘法 求解

    千次阅读 2019-07-12 12:25:50
    线性最小二乘法 ...非线性最小二乘法求解 https://zhuanlan.zhihu.com/p/42383070 很多问题最终归结为一个最小二乘问题,如SLAM算法中的Bundle Adjustment,位姿图优化等等。求解最小二乘的方法有很多,...

    线性最小二乘法

    可逆矩阵条件:

    1.A 为 方阵

    2. det(A) 不等于 0  即 A是满秩矩阵 即 A是非奇异矩阵

     

     

     

     

    非线性最小二乘法求解

    https://zhuanlan.zhihu.com/p/42383070

    很多问题最终归结为一个最小二乘问题,如SLAM算法中的Bundle Adjustment,位姿图优化等等。求解最小二乘的方法有很多,高斯-牛顿法就是其中之一。

     

    推导

    对于一个非线性最小二乘问题: 

    [公式]

    高斯牛顿的思想是把 [公式] 利用泰勒展开,取一阶线性项近似。

    [公式]

    带入到(1)式:

    [公式]

    对上式求导,令导数为0。

    [公式]

    令 [公式] , [公式] , 式(4)即为

    [公式]

    求解式(5),便可以获得调整增量 [公式] 。这要求 [公式]可逆(正定),但实际情况并不一定满足这个条件,因此可能发散,另外步长[公式]可能太大,也会导致发散。

    使用 H 代替了 二阶导数 海参矩阵,减少了计算量,,,(二阶展开就是海参矩阵)

    综上,高斯牛顿法(Gauss-Newton)的步骤为

    STEP1. 给定初值 [公式]
    STEP2. 对于第k次迭代,计算 雅克比 [公式] , 矩阵[公式] , [公式] ;根据(5)式计算增量 [公式] ;
    STEP3. 如果 [公式] 足够小,就停止迭代,否则,更新 [公式] .
    STEP4. 循环执行STEP2. SPTE3,直到达到最大循环次数,或者满足STEP3的终止条件。

    编程实现

    问题:非线性方程: [公式] ,给定n组观测数据 [公式] ,求系数 [公式] .

    分析:令 [公式] ,N组数据可以组成一个大的非线性方程组

    [公式]

    我们可以构建一个最小二乘问题:

    [公式] .

    要求解这个问题,根据推导部分可知,需要求解雅克比。

    [公式]

    使用推导部分所述的步骤就可以进行解算。代码实现:

    ydsf16/Gauss_Newton_solver​github.com

     


    /**
     * This file is part of Gauss-Newton Solver.
     *
     * Copyright (C) 2018-2020 Dongsheng Yang <ydsf16@buaa.edu.cn> (Beihang University)
     * For more information see <https://github.com/ydsf16/Gauss_Newton_solver>
     *
     * Gauss_Newton_solver is free software: you can redistribute it and/or modify
     * it under the terms of the GNU General Public License as published by
     * the Free Software Foundation, either version 3 of the License, or
     * (at your option) any later version.
     *
     * Gauss_Newton_solver is distributed in the hope that it will be useful,
     * but WITHOUT ANY WARRANTY; without even the implied warranty of
     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
     * GNU General Public License for more details.
     *
     * You should have received a copy of the GNU General Public License
     * along with Gauss_Newton_solver. If not, see <http://www.gnu.org/licenses/>.
     */
    
    #include <iostream>
    #include <eigen3/Eigen/Core>
    #include <vector>
    #include <opencv2/opencv.hpp>
    #include <eigen3/Eigen/Cholesky>
    #include <eigen3/Eigen/QR>
    #include <eigen3/Eigen/SVD>
    #include <chrono>
    
    /* 计时类 */
    class Runtimer{
    public:
        inline void start()
        {
            t_s_  = std::chrono::steady_clock::now();
        }
        
        inline void stop()
        {
            t_e_ = std::chrono::steady_clock::now();
        }
        
        inline double duration()
        {
            return std::chrono::duration_cast<std::chrono::duration<double>>(t_e_ - t_s_).count() * 1000.0;
        }
        
    private:
        std::chrono::steady_clock::time_point t_s_; //start time ponit
        std::chrono::steady_clock::time_point t_e_; //stop time point
    };
    
    /*  优化方程 */
    class CostFunction{
    public:
            CostFunction(double* a, double* b, double* c, int max_iter, double min_step, bool is_out):
            a_(a), b_(b), c_(c), max_iter_(max_iter), min_step_(min_step), is_out_(is_out)
            {}
            
            void addObservation(double x, double y)
            {
                std::vector<double> ob;
                ob.push_back(x);
                ob.push_back(y);
                obs_.push_back(ob);
            }
            
            void calcJ_fx()
            {
                J_ .resize(obs_.size(), 3);
                fx_.resize(obs_.size(), 1);
                
                for ( size_t i = 0; i < obs_.size(); i ++)
                {
                    std::vector<double>& ob = obs_.at(i);
                    double& x = ob.at(0);
                    double& y = ob.at(1);
                    double j1 = -x*x*exp(*a_ * x*x + *b_*x + *c_);
                    double j2 = -x*exp(*a_ * x*x + *b_*x + *c_);
                    double j3 = -exp(*a_ * x*x + *b_*x + *c_);
                    J_(i, 0 ) = j1;
                    J_(i, 1) = j2;
                    J_(i, 2) = j3;
                    fx_(i, 0) = y - exp( *a_ *x*x + *b_*x +*c_);
                }
            }
           
           void calcH_b()
            {
                H_ = J_.transpose() * J_;
                B_ = -J_.transpose() * fx_;
            }
            
            void calcDeltax()
            {
                deltax_ = H_.ldlt().solve(B_); 
            }
           
           void updateX()
            {
                *a_ += deltax_(0);
                *b_ += deltax_(1);
                *c_ += deltax_(2);
            }
            
            double getCost()
            {
                Eigen::MatrixXd cost= fx_.transpose() * fx_;
                return cost(0,0);
            }
            
            void solveByGaussNewton()
            {
                double sumt =0;
                bool is_conv = false;
                for( size_t i = 0; i < max_iter_; i ++)
                {
                    Runtimer t;
                    t.start();
                    calcJ_fx();
                    calcH_b();
                    calcDeltax();
                    double delta = deltax_.transpose() * deltax_;
                    t.stop();
                    if( is_out_ )
                    {
                        std::cout << "Iter: " << std::left <<std::setw(3) << i << " Result: "<< std::left <<std::setw(10)  << *a_ << " " << std::left <<std::setw(10)  << *b_ << " " << std::left <<std::setw(10) << *c_ << 
                        " step: " << std::left <<std::setw(14) << delta << " cost: "<< std::left <<std::setw(14)  << getCost() << " time: " << std::left <<std::setw(14) << t.duration()  <<
                        " total_time: "<< std::left <<std::setw(14) << (sumt += t.duration()) << std::endl;
                    }
                    if( delta < min_step_)
                    {
                        is_conv = true;
                        break;
                    }
                    updateX();
                }
               
               if( is_conv  == true)
                    std::cout << "\nConverged\n";
                else
                    std::cout << "\nDiverged\n\n";
            }
            
            Eigen::MatrixXd fx_;
            Eigen::MatrixXd J_; // 雅克比矩阵
            Eigen::Matrix3d H_; // H矩阵
            Eigen::Vector3d B_;
            Eigen::Vector3d deltax_;
            std::vector< std::vector<double>  > obs_; // 观测
            double* a_, *b_, *c_;
            
            int max_iter_;
            double min_step_;
            bool is_out_;
    };//class CostFunction
    
    
    
    
    int main(int argc, char **argv) {
        
        const double aa = 0.1, bb = 0.5, cc = 2; // 实际方程的参数
        double a =0.0, b=0.0, c=0.0; // 初值
        
        /* 构造问题 */
        CostFunction cost_func(&a, &b, &c, 50, 1e-10, true);
        
        /* 制造数据 */
        const size_t N = 100; //数据个数
        cv::RNG rng(cv::getTickCount());
        for( size_t i = 0; i < N; i ++)
        {
            /* 生产带有高斯噪声的数据 */
           double x = rng.uniform(0.0, 1.0) ;
           double y = exp(aa*x*x + bb*x + cc) + rng.gaussian(0.05);
           
           /* 添加到观测中 */
           cost_func.addObservation(x, y);
        }
        /* 用高斯牛顿法求解 */
        cost_func.solveByGaussNewton();
        return 0;
    }
    展开全文
  • 导航与定位问题,泰勒展开法与最小二乘法求解TDOA
  • 最小二乘法求解最优化问题 A solution is obtained by a Fletcher version of the Levenberg-Maquardt algoritm for minimization of a sum of squares of equation residuals.
  • 最小二乘法求解傅里叶级数系数 **自己推导的,如有错误,请大家批评指出 ** 题目 Solution

    最小二乘法求解傅里叶级数系数

    自己推导的,如有错误,请大家批评指出,谢谢!

    题目

    在这里插入图片描述

    Solution

    要求解 x ^ \hat{x} x^, 相当于找到一个合适的 x ^ \hat{x} x^使得估计的测量误差 H x ^ − z H\hat{x}-z Hx^z最小。
    可以采用欧几里得范数 ∣ H x ^ − z ∣ |H\hat{x}-z| Hx^z来表征估计误差,或者等价采用其平方形式 ε 2 ( x ^ ) = ∣ H x ^ − z ∣ 2 = Σ i = 1 m [ Σ j = 1 n h i j x j ^ − z i ] 2 \varepsilon^2(\hat{x})=|H\hat{x}-z|^2 = \Sigma^{m}_{i=1}[\Sigma^{n}_{j=1}h_{ij}\hat{x_j}-z_i]^2 ε2(x^)=Hx^z2=Σi=1m[Σj=1nhijxj^zi]2
    该微分函数在所有关于 x k ^ \hat{x_k} xk^的导数等于零的点取得最小值。对 x k x_k xk求偏导即可。

    在这里插入图片描述

    展开全文
  • 最小二乘法求解一元线性回归 介绍线性回归模型以及简单一元线性回归模型的解法。 通过代码实现最小二乘法求解一元线性回归实例,并对结果进行预测。 一、线性回归 二、回归问题的解决 三、最小二乘法介绍 四、...

    最小二乘法求解一元线性回归

    介绍线性回归模型以及简单一元线性回归模型的解法。

    通过代码实现最小二乘法求解一元线性回归实例,并对结果进行预测。

    一、线性回归

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

    在这里插入图片描述

    二、回归问题的解决

    在这里插入图片描述

    三、最小二乘法介绍

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

    四、最小二乘法求解线性回归

    在这里插入图片描述

    五、实例验证

    案例背景:数据中参数x为学习时间,y为得分。通过最小二乘法求解参数w,b,均方差。并预测x=80时的得分。

    数据链接:
    链接: https://pan.baidu.com/s/1KVw_9O5o9vqQnpgRNfLGVQ
    提取码:8u8e

    1.导入数据

    # 导入必要库
    import numpy as np
    import matplotlib.pyplot as plt
    
    points = np.genfromtxt('E:/PythonData/machine_learning/data.csv',delimiter=',')
    
    # 查看前5行数据
    points[:5]
    

    在这里插入图片描述

    2.绘制散点图

    # 分别提取points中的x和y数据
    x = points[:,0]
    y = points[:,1]
    
    # 绘制散点图
    plt.scatter(x,y)
    plt.show()
    

    在这里插入图片描述

    3.定义损失函数

    # 损失函数是系数w,b的函数,另外还要传入数据x,y
    def computer_cost(w,b,points):
        total_cost = 0
        M = len(points)
      
        # 逐点计算平方损失误差,然后求平均数
        for i in range(M):
            x = points[i,0]
            y = points[i,1]
            total_cost +=(y - w *x - b)**2
        # 取平均
        return total_cost/M
    

    4.定义算法拟合函数

    # 先定义求均值函数
    def average(data):
        sum = 0
        num = len(data)
        for i in range(num):
            sum+=data[i]
        return sum/num
        
    # 定义核心拟合函数
    def fit(points):
        M = len(points)
        x_bar = average(points[:,0])
        
        sum_yx = 0
        sum_x2 = 0
        sum_delta = 0
        
        for i in range(M):
            x = points[i,0]
            y = points[i,1]
            sum_yx += y*(x-x_bar)
            sum_x2 += x**2
            
        # 根据公式计算w
        w = sum_yx/(sum_x2 - M*(x_bar**2))
        
        # 再次创建for循环计算b
        for i in range(M):
            x = points[i,0]
            y = points[i,1]
            sum_delta += y-w*x
        b = sum_delta/M
        return w,b 
    

    5.测试(得到参数w,b,均方误差)

    # 将测试集传入拟合函数中
    w,b = fit(points)
    print('w is :',w)
    print('b is :',b)
    
    cost = computer_cost(w,b,points)
    print('cost is ',cost)
    

    在这里插入图片描述

    6.绘制拟合曲线

    plt.scatter(x,y)
    
    # 针对每一个x,绘制出预测的值
    pred_y = w*x+b
    
    # 画出拟合曲线,颜色设置为红色
    plt.plot(x,pred_y,c='r')
    

    在这里插入图片描述

    7.预测分数

    # 给出参数x,得出预测结果
    pred_y1 = w*80+b
    print(pred_y1)
    

    在这里插入图片描述

    展开全文
  • 最小二乘法求解线性方程组与伪逆

    千次阅读 2020-12-07 13:30:50
    最小二乘法求解线性方程组与伪逆 对于线性方程组Ax=bAx=bAx=b的求解。 如果A是可逆的,我们可以通过方程式左右两边乘A−1A^{-1}A−1求解: x=A−1bx=A^{-1}bx=A−1b 但是如果AAA是不可逆的方阵或非方阵呢? 这时可以...

    最小二乘法求解线性方程组与伪逆

    对于线性方程组 A x = b Ax=b Ax=b的求解。
    如果A是可逆的,我们可以通过方程式左右两边乘 A − 1 A^{-1} A1求解:
    x = A − 1 b x=A^{-1}b x=A1b
    但是如果 A A A是不可逆的方阵或非方阵呢?
    这时可以使用最小二乘法求解:

    首先我们来看看一些矩阵求导的相关性质:

    在这里插入图片描述

    在这里插入图片描述

    接下来我们用最小二乘法来求解方程组的解 x ∗ x^* x

    x ∗ = a r g max ⁡ x ∣ ∣ A x − b ∣ ∣ 2 x^*=arg\max_x||Ax-b||^2 x=argxmaxAxb2
    f ( x ) = ∣ ∣ A x − b ∣ ∣ 2 f(x)=||Ax-b||^2 f(x)=Axb2

    f ( x ) f(x) f(x)的极值,求其梯度: ∇ f ( x ) = 0 \nabla f(x)=0 f(x)=0

    ∇ f ( x ) = ∂ ∂ x ∣ ∣ A x − b ∣ ∣ 2 \nabla f(x)=\frac{\partial}{\partial x}||Ax-b||^2 f(x)=xAxb2 = ∂ ∂ x ( A x − b ) T ( A x − b ) =\frac{\partial}{\partial x}(Ax-b)^T(Ax-b) =x(Axb)T(Axb) = ∂ ∂ x ( x T A T − b T ) ( A x − b ) =\frac{\partial}{\partial x}(x^TA^T-b^T)(Ax-b) =x(xTATbT)(Axb) = ∂ ∂ x ( x T A T A x − x T A T b − b T A x + b T b ) =\frac{\partial}{\partial x}(x^TA^TAx-x^TA^Tb-b^TAx+b^Tb) =x(xTATAxxTATbbTAx+bTb) = 2 A T A x − 2 A T b = 0 =2A^TAx-2A^Tb=0 =2ATAx2ATb=0

    则解得: x ∗ = ( A T A ) - 1 A T b x^*=(A^TA)^{-1}A^Tb x=(ATA)1ATb

    ( A T A ) - 1 A T (A^TA)^{-1}A^T (ATA)1AT A A A的伪逆。

    参考:矩阵论简明教程, 徐仲 张凯院 陆全 冷国伟, 科学出版社

    展开全文
  • 一个最小的最小二乘法求解器,针对反复出现的小型密集问题的最小化,可提供极佳的性能。 为什么选择微型求解器? 以下是Tiny Solver与其他最小二乘最小化实现的不同之处 用户的成本函数直接内联到Levenberg-...
  • 采用递推最小二乘法求解超定线性方程组 Ax=b,其中 A 为 mxn 维的已知矩阵,b 为 m 维的已知向量,x 为 n 维的未知向量,其中 n=10,m=10000。A 与 b 中的元素服从独立同 分布的正态分布。绘出横坐标为迭代步数时的...
  • clearclcx=[20 25 30 35 40 20 25 30 35 40 20 25 30 35 40 20 25 30 35 40 20 25 30 35 40];y=[10 10 10 10 10 20 20 20 20 20 30 30 30 30 30 40 40 40 40 40 50 50 50 50 50];z=[11423 13000 12209 12826 13302 1...
  • 利用labview制作的简单最小二乘法求解线性方程工具 简单易操作 纯软件制作
  • 在图像变换中用最小二乘法求解仿射变换参数
  • 用于绘制和研究使用最小二乘法求解 ODE 的结果的脚本。 Rnorm 没有被微分以获得系数,而是我们搜索最小化 Rnorm 的 Cj。
  • 机器学习之线性回归的最小二乘法求解 假设现在一个普通的一阶线性方程,y=2*x+2*t。t是随机噪音,生成的散列点(x,y)会沿直线y=2*x上下摆动。利用最小二乘法做一次简单的一阶“曲线”拟合。用matlab做数据...
  • 最小二乘法求解步骤

    万次阅读 多人点赞 2019-09-03 23:02:38
    目标:线性回归问题,找到最佳参数使得损失函数最小 一、损失函数定义 线性方程:y=ax+by=ax+b y=ax+b 对于每个样本点 x(i)x(i) x^{(i)} ,其预测值为 y^(i)=ax(i)+by^(i)=ax(i)+b \hat y^{(i)}=ax^{(i)}+b ...
  • 基于jupyter notebook的python编程—–利用梯度下降算法求解多元线性回归方程,并与最小二乘法求解进行精度对比目录一、梯度下降算法的基本原理1、梯度下降算法的基本原理二、题目、表格数据、以及python环境搭建1、...
  • 使用Python编程,分别根据梯度法和最小二乘法求解多元函数问题分别使用梯度下降法和最小二乘法求解多元函数并进行比较,这里使用jupyter notebook平台进行Python编程一、题目描述二、使用梯度下降法求解多元函数(一...
  • 求解极值的第一步往往是求解一阶导数并让一阶导数等于0,最小二乘法也不能免俗。因此,首先在残差平方和RSS上对参数向量求导。这里的过程涉及少数矩阵求导的内容,需要查表来确定,感兴趣可以去维基百科去查看矩阵...
  • 可用于求非线性方程的最优解,设定初始值的范围,通过优化模型计算出满足范围条件下的最优解。
  • Python 最小二乘法求解线性回归模型 机器学习线性回归模型 线性回归(linear regression)是一种线性模型,它假设输入变量 x 和单个输出变量 y 之间存在线性关系 具体来说,利用线性回归模型,可以从一组输入变量...
  • 利用梯度下降算法求解多元线性回归方程,并与最小二乘法求解进行精度对比 基于jupyter notebook的python编程 这里写目录标题利用梯度下降算法求解多元线性回归方程,并与最小二乘法求解进行精度对比基于jupyter ...
  • 【机器学习】简单线性回归&最小二乘法十分钟学会!...seid=7923451902729373671&spm_id_from=333.337.0.0 ...最小二乘法求解线性回归 每个实际值和估计值的差的平方和,并希望值越小越好 平方和越大,则估计值和.
  • 以影像数据为例: 像素坐标(pixelX,pixelY)与实际地理坐标(geoX,geoY)之间的仿射变换函数为: geoX=a*pixelX+b*pixelY+c...如果知道一组像素坐标和实际地理坐标,可以通过最小二乘法,计算出abcdef六个参数。 最小..
  • 最小二乘法求解直线方程系数

    万次阅读 2017-06-03 09:03:19
    引言最小二乘法是经典的参数稳健估计方法。核心思想是使得估计出的模型与实际数据之间误差的平方和最小(趋于+∞\infty)。直线参数的估计又是最常用的曲线参数估计,稳健估计方法有很多:RANSAC、BaySAC、极大似然...
  • 提出了基于最小二乘优化系统的估值方法测量求解光器件物理偏振参量。对于输入偏振态可测或非可测的系统,在建立相应的估值解析方程和适当的测量条件基础上,此测量方案均是简单易行的。对以旋转波片为偏振发生器的...
  • 最小二乘法的简单lingo运算,仅供参考,请大家指证!
  • 最小二乘法求解多项式系数

    千次阅读 2020-06-05 11:47:22
    1原理 统计学习方法21页 2实现 #n次多项式的...array**(i+1)).reshape((m,1))]) from numpy.linalg import solve #求解系数得到最小的误差 X = solve(np.dot(A.T,A),np.dot(A.T,y_array.T)) print(X) 作图 自己实现看吧

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 15,228
精华内容 6,091
关键字:

最小二乘法求解