最小二乘法

最小二乘法

1. 最小二乘法的本质是什么
   以用不同的尺子量长度为例，得到几个测量值，如何确定真值？“让误差的平方和最小的y就是真值。”这是基于如果误差是随机的，测量值会在真值上下波动。尺子的例子，优化目标$min\sum (y-y_i{)}^2$，则目标函数对y求导，导数为零的y恰好为算数平均值。least square --> 最小+平方
2. 推广：用于多项式拟合
3. MATLAB实现
   3.1 解方程
   最小二乘法原理与MATLAB代码——线性拟合、多项式拟合、指定函数拟合
   3.1.1 直线方程
   $min\sum (y_i-a{x}_i-b{)}^2=0$，对a,b分别求偏导为0，得下面的方程组，即可解出a,b。另外可将方程组写成矩阵形式。
   $$\{\begin{array}{c}-\sum 2\left(y_i-a{x}_i-b\right)x_i=0\\ -\sum 2\left(y_i-a{x}_i-b\right)=0\end{array}$$ $$\left(\begin{array}{c}\begin{array}{cc}\sum 1& \sum x_i\end{array}\\ \begin{array}{cc}\sum x_i& \sum x_i^2\end{array}\end{array}\right)\left(\begin{array}{c}b\\ a\end{array}\right)=\left(\begin{array}{c}\sum y_i\\ \sum y_i{x}_i\end{array}\right)$$ 3.1.2 多项式，同样可写成类似上面的矩阵形式。XA=Y, 则A=inv(X)*Y
   3.2 对最优化问题，可以用梯度下降法求解。但下面的代码没有优化步长，费时间。

clear;
clf;
x0 = linspace(-1,1,100);
% 函数+噪声
y0 = sin(2*x0).*((x0.^2-1).^3+1/2)+(0.2*rand(1,100)-0.1);
% 原始数据图像
plot(x0,y0,'.');
% 用多项式函数来拟合f(x)
% f(x)=a0+a1*x+...+ak*x^k
% 目标函数：误差平方和 sum(y-f(x))^2
% 求最优值：对各个系数的偏导为0
% 法1：求解线性方程 XA=Y，则A=Y/X
% 法2：梯度下降法
% 指定多项式系数k
order = 10;
ini_A = zeros(1,order+1);  % 初始点
% 目标函数是最小化误差平方和，初始点给定，用梯度下降法求最优解
[A,k,err_ols] = GradientMethod_ols(x0, y0, order, ini_A);
y1 = polyfun(x0,order,A);
hold on;
plot(x0,y1);

function [A,k,err_ols] = GradientMethod_ols(x0, y0, order, ini_A)
% x,y：原始数据点
% order,ini_A: 确定了一个多项式函数
%% 第一步：给定初始迭代点A，收敛精度err，k=1
k = 1;
iterk = 1000;
A = ini_A;
err = 0.5;
%% 第二步:计算梯度和模并取搜索方向
err_ols = 100;
while(k<iterk)
    grad = gradofpolyols(x0,y0,order,A);
    dir = -grad;  % 搜索方向为负梯度方向
    %% 第三步:进行收敛判断
    if err_ols<=err  % 梯度等于0
        disp('误差足够小');
        break;
    end
%% 第四步求最优步长以及求新迭代点 
    lmd_opt = 0.1;  % 直接取固定值了
    A0 = A;
    while(1)   
    A = A0 + dir*lmd_opt;
    err_ols_1 = sum((y0-polyfun(x0,order,A)).^2);
    if err_ols_1<err_ols
        err_ols = err_ols_1;
        break
    end
    lmd_opt = lmd_opt*0.5;    
    end
    k = k+1;
end
if k >=iterk
disp('迭代次数达上限！');
end
end

function grad = gradofpolyols(x0,y0,order,A)
% 求梯度
grad = zeros(1,length(A));
for i = 1:length(A)
    grad(i) = -2*(y0-polyfun(x0,order,A))*(x0.^(i-1))';
end
end

function y = polyfun(x,order,A)
% 多项式函数的值
% y=a0+a1*x+...+ak*x^k
% x 自变量
% k 多项式阶数
% A 系数向量
y = zeros(1,length(x));
for i =0:order
    y = y + A(i+1)*x.^i;
end
end

4. 补充
   直线拟合（W=inv(A’*A)*A’*Y // W=inv( R )*Q’*Y）与lsqcurvefit函数
   MATLAB之最小二乘法
   b.

求多元函数的极值点

已经忘记多元函数的极值是对各项求偏导，所以复习一下。
1.一阶导数为0的点称为驻点，驻点不一定是极值点（$y=x^3$），极值点不一定是驻点（$y=|x|$）；对于多元函数，驻点是所有一阶导为零的点。
2.曲线凹凸性改变的点为拐点，如果函数二阶连续可导，拐点可推出二阶导为0；
3.极值点是函数图像的某段子区间内上极大值或者极小值点的横坐标，出现在驻点或者不可导点处。
多元函数的极值与其求法
（$f(x)=x^2,f^{\prime \prime }(x{)}^2=4>0,f^{\prime \prime }(x)=2>0有极小值$）
如果已知函数有且只有一个极值点，那“求偏导为0的点，可以得到极值点”，也对喽？