[含matlab程序]最小二乘法线性（以空间平面为例）与非线性（以二次函数为例）拟合，不看后悔_matlab最小二乘法拟合非线性方程-爱代码爱编程

最小二乘法解决的是什么问题？使用最小二乘法求参数的关键是什么？
○最小二乘法用于处理数据拟合问题，本质是实现向量向$\mathrm{X}$矩阵的列空间投影
○关键是找$\mathrm{X}$矩阵

一、以直线拟合作为切入点进行解释

做两个假设
1.我们手中有m个散点，要想用一条直线去拟合，散点的形式是：
X = [ x 0 x 1 ⋮ x m − 1 ] Y = [ y 0 y 1 ⋮ y m − 1 ] X = {\begin{bmatrix} {{x_0}}\\ {{x_1}}\\ \vdots \\ {{x_{m - 1}}} \end{bmatrix}} Y = {\begin{bmatrix} {{y_0}}\\ {{y_1}}\\ \vdots \\ {{y_{m - 1}}} \end{bmatrix}} X= ​x0​x1​⋮xm−1​​ ​Y= ​y0​y1​⋮ym−1​​ ​

2.假如这条直线的形式是已知的，如下：
h = θ 0 x + θ 1 h = {\theta _0}x+{\theta _1} h=θ0​x+θ1​

将我们已知的点的自变量x代入，可以得到以下形式，因此由这条拟合直线可以给我们输出m个结果组成向量H：

h 0 = θ 0 x 0 + θ 1 h 1 = θ 0 x 1 + θ 1 ⋮ h m − 1 = θ 0 x m − 1 + θ 1 \begin{array}{l} {h_0} = {\theta _0}{x_0} + {\theta _1}\\ {h_1} = {\theta _0}{x_1} + {\theta _1}\\ \vdots \\ {h_{m-1}} = {\theta _0}{x_{m-1}} + {\theta _1} \end{array} h0​=θ0​x0​+θ1​h1​=θ0​x1​+θ1​⋮hm−1​=θ0​xm−1​+θ1​​
H = [ h 0 h 1 ⋮ h m − 1 ] = [ x 0 1 x 1 1 ⋮ ⋮ x m − 1 1 ] [ θ 0 θ 1 ] = X θ (1) {H}= \begin{bmatrix} h_0 \\ h_1 \\ \vdots \\ h_{m-1} \\ \end{bmatrix}=\begin{bmatrix} x_0 &1\\ x_1 &1\\ \vdots &\vdots\\ x_{m-1} &1\\ \end{bmatrix}\begin{bmatrix} \theta_0 \\ \theta_1 \\ \end{bmatrix}={\rm X}{\theta}\tag{1} H= ​h0​h1​⋮hm−1​​ ​= ​x0​x1​⋮xm−1​​11⋮1​ ​[θ0​θ1​​]=Xθ(1)
注意这个$\mathrm{X}$从原来的$X$向量的基础上扩展成mx2的矩阵。
用下面这个公式描述直线的输出值和真实散点的y值的差距：

J ( θ ) = ∥ H − Y ∥ 2 = ∥ H − Y ∥ 2 = ∥ X θ − Y ∥ 2 J(\theta) = {\left\| {H - Y} \right\|^2} = {\left\| {H - Y} \right\|^2} = {\left\| {{\rm X}\theta - Y} \right\|^2} J(θ)=∥H−Y∥2=∥H−Y∥2=∥Xθ−Y∥2
最小二乘法的目的就在于实现解出$\theta$矩阵，使$J(\theta )$最小化，即：
arg ⁡ min ⁡ θ J ( θ ) \mathop {\arg \min }\limits_\theta J(\theta ) θargmin​J(θ)
具体是通过矩阵求导实现解出$\theta$矩阵的，参考这篇文章：一文让你彻底搞懂最小二乘法（超详细推导）
这里直接给出结果：
θ = ( X T X ) − 1 X T Y (2) \theta = {({{\rm X}^T}{\rm X})^{ - 1}}{{\rm X}^T}Y\tag{2} θ=(XTX)−1XTY(2)
这个形式让人联想起投影矩阵。
H = X θ = X ( X T X ) − 1 X T Y = P Y H = {\rm X}\theta ={\rm X}{({{\rm X}^T}{\rm X})^{ - 1}}{{\rm X}^T}Y={P}Y H=Xθ=X(XTX)−1XTY=PY
其中的P就是常说的投影矩阵，即：
P = X ( X T X ) − 1 X T P={\rm X}{({{\rm X}^T}{\rm X})^{ - 1}}{{\rm X}^T} P=X(XTX)−1XT
从形式上看，拟合的输出向量$H$实际上是真实散点的因变量（对于直线来说就是$Y$）向量向自变量$\mathrm{X}$矩阵的列空间投影的结果。
对于直线而言，$\mathrm{X}$的列空间就是由下面两个向量张成的空间：
[ x 0 x 1 ⋮ x m − 1 ] 和 [ 1 1 ⋮ 1 ] \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{m-1} \\ \end{bmatrix}和\begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \\ \end{bmatrix} ​x0​x1​⋮xm−1​​ ​和 ​11⋮1​ ​

向量的投影可以用来对向量进行正交化。
矩阵的投影则可以把向量投影到空间。
线性代数（十三）投影
对于矩阵
[ 3 0 0 4 0 0 ] \begin{bmatrix} 3&0 \\ 0&4 \\ 0&0 \\ \end{bmatrix} ​300​040​ ​
张成的空间，在下图中即是由向量$c$和向量$b$张成的黄色平面。图上标注的是一个三角，实际是一个平面，这个平面中的向量都可以用向量$c$和向量$b$线性组合的形式表达。

如果向量$a$要向这个平面投影，根据投影矩阵：

H = X ( X T X ) − 1 X T Y = P Y H ={\rm X}{({{\rm X}^T}{\rm X})^{ - 1}}{{\rm X}^T}Y={P}Y H=X(XTX)−1XTY=PY
在本题中：
Y = [ 0 3 4 ] X = [ 3 0 0 4 0 0 ] Y=\begin{bmatrix} 0 \\ 3 \\ 4\\ \end{bmatrix} {\rm X}=\begin{bmatrix} 3&0\\ 0&4 \\ 0&0\\ \end{bmatrix} Y= ​034​ ​X= ​300​040​ ​
因此可以求出投影到$\mathrm{X}$张成平面上的结果是:
H = [ 0 3 0 ] H=\begin{bmatrix} 0\\ 3\\ 0\\ \end{bmatrix} H= ​030​ ​
用这样一段程序（matlab）就能快速计算：

x = [3,0;0,4;0,0]
%x = [0,3;4,0;0,0]%交换次序因为张成空间不变，投影结果也不变
y = [0;3;4]
h = x/(x'*x)*x'*y

在这里顺势介绍一下向量的正交化，这块和最小二乘法关系不大但是和投影关系比较大。
如下图，向量$a$和向量$b$原本不是正交的（正交要求内积为0），$H$是$a$在$b$上的投影，上一步我们已经求得，而$e=a-H$且与$b$是正交的，$e$很好求得，由于向量是可以平移的，把红色的$e$向量的尾巴平移到原点处，就可以得到过原点的一个和$b$正交的向量，这是对$a$的正交化。

如果引入一个新的向量$d$，要想使$d$与正交化$a$得到的 e 1 e_1 e1​以及$b$都正交，要先作出$d$到 e 1 e_1 e1​与$b$张成平面的投影（下图中是一个蓝色的向量，没有字母标注），再用$d$减去这个投影得到的向量就可以解出向量 e 2 e_2 e2​。最终 e 1 e 2 e_1e_2 e1​e2​和$b$线性无关且两两正交。这三个向量分别除以各自的模长就得到两两正交的单位向量，这个操作也就是常说的单位化。

图片是本文原创，使用请标明出处。

二、引入平面的参数求解

注意平面环节中的符号和第一部分的直线环节的可能有相似之处，但意义不相同。
用不同的形式拟合散点，关键是构建$\mathrm{X}$矩阵，对于平面，$\mathrm{X}$矩阵和直线是不同的。
平面的形式如下：
a x + b y + c z + d = 0 或者用一个和直线环节介绍的类似的形式： h = θ 0 x 0 + θ 1 x 1 + θ 2 \begin{array}{l} ax + by +cz+d=0\\ \\ 或者用一个和直线环节介绍的类似的形式：\\ \\ h = {\theta _0}{x_0} + {\theta _1}{x_1} + {\theta _2} \end{array} ax+by+cz+d=0或者用一个和直线环节介绍的类似的形式：h=θ0​x0​+θ1​x1​+θ2​​
可以看出如果找到一个平面去拟合散点，这个平面的因变量变成两个自变量线性组合的形式，而对于直线只有一个因变量。两个自变量 x 0 {x_0} x0​， x 1 {x_1} x1​分别对应着平面一般式中的$x$和$y$，而因变量$h$则代表着一般式中的$z$。
如果这时候我们手中有很多散点，把他们带入到平面方程中：
h 0 = θ 0 x 0 , 0 + θ 1 x 0 , 1 + θ 2 h 1 = θ 0 x 1 , 0 + θ 1 x 1 , 1 + θ 2 ⋮ h m − 1 = θ 0 x m − 1 , 0 + θ 1 x m − 1 , 1 + θ 2 \begin{array}{l} {{h_0} = {\theta _0}{x_{0,0}} + {\theta _1}{x_{0,1}} + {\theta _2}}\\ {{h_1} = {\theta _0}{x_{1,0}} + {\theta _1}{x_{1,1}} + {\theta _2}}\\ \vdots \\ {{h_{m - 1}} = {\theta _0}{x_{m - 1,0}} + {\theta _1}{x_{m - 1,1}} + {\theta _2}} \end{array} h0​=θ0​x0,0​+θ1​x0,1​+θ2​h1​=θ0​x1,0​+θ1​x1,1​+θ2​⋮hm−1​=θ0​xm−1,0​+θ1​xm−1,1​+θ2​​

待求的参数是 θ 0 {\theta _0} θ0​， θ 1 {\theta _1} θ1​， θ 2 {\theta _2} θ2​，同样可以写成$H=\mathrm{X}\theta$的形式，对于求解平面的这三个参数来说，$\mathrm{X}$矩阵是长得这个样子:
[ x 0 , 0 x 0 , 1 1 x 1 , 0 x 1 , 1 1 ⋮ ⋮ ⋮ x m − 1 , 0 x m − 1 , 1 1 ] \begin{bmatrix} x_{0,0} &x_{0,1}&1\\ x_{1,0} &x_{1,1}&1\\ \vdots &\vdots &\vdots\\ x_{m-1,0} &x_{m-1,1}&1\\ \end{bmatrix} ​x0,0​x1,0​⋮xm−1,0​​x0,1​x1,1​⋮xm−1,1​​11⋮1​ ​
而$Y$矩阵就是我们手中散点的因变量，对于平面来说就是$z$值，注意区分符号。
现在$\mathrm{X}$和$Y$矩阵都是已知的，根据$(2)$式即可求解出$\theta$矩阵（其实是向量）:
[ θ 0 θ 1 θ 2 ] \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \end{bmatrix} ​θ0​θ1​θ2​​ ​

三、引入非线性表达式的参数求解

假如手中一些散点，要用二次函数进行拟合，一个二次函数的典型形式如下：
h = θ 0 x 2 + θ 1 x + θ 2 h = {\theta _0}x^2+{\theta _1} x+{\theta _2} h=θ0​x2+θ1​x+θ2​
对于非线性的场景，关键仍然是找$\mathrm{X}$矩阵，首先还是考虑怎么写成$(1)$式的形式：
H = [ h 0 h 1 ⋮ h m − 1 ] = [ x 0 2 x 0 1 x 1 2 x 1 1 ⋮ ⋮ ⋮ x m − 1 2 x m − 1 1 ] [ θ 0 θ 1 θ 2 ] = X θ {H}= \begin{bmatrix} h_0 \\ h_1 \\ \vdots \\ h_{m-1} \\ \end{bmatrix}=\begin{bmatrix} x_0^2 &x_0&1\\ x_1^2 &x_1&1\\ \vdots &\vdots&\vdots\\ x_{m-1}^2 &x_{m-1}&1\\ \end{bmatrix}\begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \\ \end{bmatrix}={\rm X}{\theta} H= ​h0​h1​⋮hm−1​​ ​= ​x02​x12​⋮xm−12​​x0​x1​⋮xm−1​​11⋮1​ ​ ​θ0​θ1​θ2​​ ​=Xθ
对于更高次幂的函数，也是如此找$\mathrm{X}$矩阵，因为对差距函数（损失函数）$J(\theta )$求导并使导函数等于零，这是一个对$\theta$求导的过程，不管是线性还是非线性，结果都是$(2)$式，因此只要找到$\mathrm{X}$矩阵，就可以求出$\theta$。

下面是matlab程序，首先介绍我自定义的LS_fit函数

LS_fit(自变量1,自变量2,因变量,是否是非线性,如果是非线性函数的话最高次项是几次)

function [theta] = LS_fit(Mat_x1,Mat_x2,Mat_x3,bool_nonlinear,exponent) 
    if bool_nonlinear == false
        %随机生成样本（点）的个数
        Mat_x4 = Mat_x1.^0;%常数项没有自变量，自变量的位置填充1
        %绘出点图
        plot3(Mat_x1,Mat_x2,Mat_x3,'redo')
        hold on

        %制作X矩阵
        x = [Mat_x1,Mat_x2,Mat_x4];
        %解出最优参数
        theta = (x'*x)\x'*Mat_x3%这里就是用投影矩阵计算投影向量

        %绘出平面图
        a = 0:1:30;
        b = 0:1:30;
        [X,Y] = meshgrid(a,b);
        Z = theta(1)*X+theta(2)*Y+theta(3);
        mesh(X,Y,Z)
        title("拟合图像")
        xlabel("x值")
        ylabel("y值")
        
    elseif (bool_nonlinear == true)
        if ischar(exponent)
            Mat_x4 = Mat_x1.^0;%常数项没有自变量，自变量的位置填充1
            x = [Mat_x1,Mat_x4];
            Mat_x2copy = Mat_x2;
            Mat_x2copy = log(Mat_x2copy);
            theta = (x'*x)\x'*Mat_x2copy;
            theta(2)=exp(theta(2));
            X = 0:0.1:2.5;
            Y = theta(2)*exp(theta(1)*X);
        else
            %制作X矩阵
            x = zeros(size(Mat_x1,1),exponent+1);
            for i = 1:(exponent+1)
                x(:,i) = Mat_x1.^(exponent-i+1);
            end
            %使用x矩阵解出参数
            theta = (x'*x)\x'*log(Mat_x2);
            X = -15:0.1:15;
            Y = 0;
            for i = 1:(exponent+1)
                Y = Y+theta(i)*X.^(exponent+1-i);
            end
        end 
        %绘出散点图
        plot(Mat_x1,Mat_x2,'o')
        hold on
        %绘出曲线图
        
        plot(X,Y,'red-')
        title("拟合图像")
        xlabel("x值")
        ylabel("y值")
        zlabel("z值")
        %text(x,y,'说明')
        legend('样本点','拟合曲线')
    else
        disp("错误的输入")
    end
end

新建一个脚本，调用刚才定义的LS_fit函数，首先展示平面拟合功能

num = 100;%随机生成100组点

Mat_x1 = 30*rand(num,1);%随机生成自变量
Mat_x2 = 20*rand(num,1);
%平面方程，z = (-3x-5y-3)/7
fuc_plane_z = (-3*Mat_x1-5*Mat_x2-3)/7;
Mat_z = fuc_plane_z-5+(5+5)*rand(num,1);%在z的基础上上下漂移
LS_fit(Mat_x1,Mat_x2,Mat_z,false,NaN) 

对于非线性的情况，以二次函数拟合为例，程序如下：

num = 100;%随机生成100组点

Mat_x1 = -15 + (30)*rand(num,1)
%用公式求x3，并随机加上正负50
offset = -50 + 100*rand(num,1)
Mat_y = 4*Mat_x1.^2+6*Mat_x1+1+offset%二次函数
%Mat_y = 5.5*Mat_x1.^3+4*Mat_x1.^2+6*Mat_x1+1+offset%三次函数
%Mat_y = 6*Mat_x1+1+offset%一次函数
LS_fit(Mat_x1,Mat_y,NaN,true,2) 

拟合的结果：

注意：即使是使用的函数的非线性模式，也可以解决一次函数拟合的问题，也就是一次函数是一种特殊的非线性函数。

num = 100;%随机生成100组点

Mat_x1 = -15 + (30)*rand(num,1)
%用公式求x3，并随机加上正负1
offset = -5 + 10*rand(num,1)
%Mat_y = 4*Mat_x1.^2+6*Mat_x1+1+offset%二次函数
%Mat_y = 5.5*Mat_x1.^3+4*Mat_x1.^2+6*Mat_x1+1+offset%三次函数
Mat_y = 6*Mat_x1+1+offset%一次函数
LS_fit(Mat_x1,Mat_y,NaN,true,1) 

针对 h = θ 0 e θ 1 x h = {\theta _0}{{\rm{e}}^{{\theta _1}x}} h=θ0​eθ1​x的形式，首先线性化，两边取对数可得 ln ⁡ h = ln ⁡ θ 0 + θ 1 x \ln h = \ln {\theta _0}{\rm{ + }}{\theta _{\rm{1}}}x lnh=lnθ0​+θ1​x，记 A = ln ⁡ θ 0 A=\ln {\theta _0} A=lnθ0​ ， h ˉ = ln ⁡ h \bar h=\ln h hˉ=lnh可得 h ˉ = θ 1 x + A \bar h = {\theta _{\rm{1}}}x{\rm{ + }}A hˉ=θ1​x+A，这是一个一次函数的形式。
当有多个样本点，需要提前把每个样本的 h ˉ \bar h hˉ解出来
h ˉ 0 = θ 1 x 0 + A h ˉ 1 = θ 1 x 1 + A ⋮ h ˉ m − 1 = θ 1 x m − 1 + A \begin{array}{l} {{\bar h_0} = {\theta _1}{x_0} +A}\\ {{\bar h_1} = {\theta _1}{x_1} +A}\\ \vdots \\ {{\bar h_{m - 1}} = {\theta _1}{x_{m - 1}} +A} \end{array} hˉ0​=θ1​x0​+Ahˉ1​=θ1​x1​+A⋮hˉm−1​=θ1​xm−1​+A​

显然$\mathrm{X}$矩阵和参数矩阵分别为：
[ x 0 1 x 1 1 ⋮ ⋮ x m − 1 1 ] 以及 [ θ 1 A ] \begin{bmatrix} x_0 &1\\ x_1 &1\\ \vdots &\vdots \\ x_{m-1} &1\\ \end{bmatrix} 以及 \begin{bmatrix} {\theta_1} \\ A\\ \end{bmatrix} ​x0​x1​⋮xm−1​​11⋮1​ ​以及[θ1​A​]
用解一次函数的方式解出两个参数， θ 0 = e A {\theta_0}=e^A θ0​=eA解出 θ 0 {\theta_0} θ0​，这样两个参数就都求解出来了。

Mat_x1 = [1,1.25,1.5,1.75,2]';
Mat_y = [5.1,5.79,6.53,7.45,8.46]';
LS_fit(Mat_x1,Mat_y,NaN,true,'EXP') 

求得参数分别是
[ θ 1 θ 0 ] = [ 0.5057 3.0725 ] \begin{bmatrix} {\theta_1} \\ {\theta_0}\\ \end{bmatrix}= \begin{bmatrix} 0.5057 \\ 3.0725\\ \end{bmatrix} [θ1​θ0​​]=[0.50573.0725​]
注意参数 θ 0 {\theta_0} θ0​和 θ 1 {\theta_1} θ1​的输出顺序，代入待拟合的函数中得到
h = 3.0725 e 0.5057 x h=3.0725{\rm{e}^{0.5057{x}}} h=3.0725e0.5057x

参考文章：
[1] 一文让你彻底搞懂最小二乘法（超详细推导）
[2] 线性代数（十三）投影
[3] QR分解

码字不易难免打错,
发现错误欢迎指出，
万分感谢！！!