自学脚手架——“Data-Driven Science and Engineering” by steven L. brunton（Chapter 4.0 - 4.7）

4.0 Overview of Machine Learning

• 相关资料：

“神经网络与深度学习” by 邱锡鹏

4.1 Classic Curve Fitting

• Over\Under-fitting

上方的图代表拟合的结果，看似拟合的恰到好处，但是放宽横坐标：

更宽的区域并没有很好拟合，这就是过拟合的一种形式。

• 为什么拟合曲线会被称为回归呢？

“回归”这个词源于弗朗西斯·高尔顿爵士（Sir Francis Galton，1822年2月16日-1911年1月17日）：

高尔顿（右）研究父子身高的方法，是线性回归的源头。高尔顿是统计学的先驱，他的衣钵传给了他的学生现代统计学奠基人卡尔·皮尔逊（Karl Pearson, 1857~1936）（左），皮尔逊相关系数的提出者，再到他的儿子埃贡·皮尔逊（Egon Pearson, 1895~1980）一脉相承，开创并极大发展了统计学。

他发现
1）高个子父母的儿子身高往往也比普通人高，但是没有他父辈高；
2）矮个子父母的儿子身高往往也比普通人矮，但是没有他父辈矮
他称这种现象为复归或回复（reversion），后来改为"回归平庸"（Regression towards mediocrity)。也就是说人类的身高都会回到平均值附近，他将这种现象称为均值回归。

高尔顿回归线，横轴显示子辈的身高，纵轴显示父辈（mid-parents父母身高均值，母亲身高乘以1.08）的身高，如果知道父辈身高，直线OM可以预测儿子身高，反之依然。

假设统计父亲（横坐标）和儿子（纵坐标）的身高数据：

由上可见在线性拟合出的曲线中，交点$(1.77,1.77)$左侧，父亲矮儿子高，右侧为父亲高儿子矮，这从另一个角度表述了高尔顿的观察结果。当然这个对应关系并不是确定关系，而是一个概率关系。

所以这条拟合出来的直线，其实就表示了均值回归现象，因此拟合直线的过程被称为线性回归（Linear Regression）。

• 但是这种回归的规律其实并不符合现实的情况。当高的父辈生了比自己矮的儿子，儿子长大变成父辈，又生比自己矮的儿子，最终子孙后代的身高会逼近普通人身高的平均值，即我们现在的人之间的身高应该都一样才对，但是事实并非如此。或者现在人的身高标准差比以前低？

我们尝试用Cobweb plot或Verhulst diagram（可以翻译为蛛网图，但是蛛网图包含了其他经常的形式）来演示这种现象：

• 下面是差分方程（映射）$x_{n+1}=\cos x_n$用Cobwebs表示的情形：

• Cobwebs画法
  给定$x_{n+1}=f(x_n)$和初始条件$x_0$，从$x_0$出发画一条垂线直到与$f$的图像相交，交点的高度就是输出$x_1$。然后可以返回到横轴，从$x_1$开始重复上面过程可以得到$x_2$。但是，一个更简单的方法是，画一条对角线$x_{n+1}=x_n$，并从交点开始做平行于横轴的线与之相交，然后从这个交点画垂线与$f$的图像相交的就是$x_2$。重复这个过程$n$次可得到在线上的前$n$个点。
• 逻辑斯蒂映射$x_{n+1}=r{x}_n(1-x_n)$的Cobwebs

现在我们可以画出身高映射$x_{n+1}=0.516x_n+0.8567$的蛛网图：

可以很容易看出其不动点为(1.77,1.77)，即均值回归。

• 拟合曲线与回归之间的关系

• 在优化问题的方法中，有一个常用的Levenberg–Marquardt algorithm

这种方法在matlab的Curve Fitting Tool中的Custom Equation使用：

matlab的fit函数的算法中也包含了这种算法：

Algorithm to use for the fitting procedure, specified as the comma-separated pair consisting of ‘Algorithm’ and either ‘Levenberg-Marquardt’ or ‘Trust-Region’.

Available when the Method is NonlinearLeastSquares.

• Code 4.1中的问题

源代码在matlab R2021b运行中有些错误，修改为下方加粗的地方

% The data
x=[1 2 3 4 5 6 7 8 9 10];
y=[0.2 0.5 0.3 3.5 1.0 1.5 1.8 2.0 2.3 2.2];
p1=fminsearch(@fit1,[1 1],[],x,y);

p2=fminsearch(@fit2,[1 1],[],x,y);
p3=fminsearch(@fit3,[1 1],[],x,y);
xf=0:0.1:11;
y1=polyval(p1,xf); y2=polyval(p2,xf); >y3=polyval(p3,xf);
subplot(2,1,2)
plot(xf,y1,‘k’), hold on
plot(xf,y2,‘k–’,‘Linewidth’,[2])
plot(xf,y3,‘k’,‘Linewidth’,[2])
plot(x,y,‘ro’,‘Linewidth’,[2]), hold on

function E=fit1(x0,x,y)
E=max(abs( x0(1)*x+x0(2)-y ));
end

function E=fit2(x0,x,y)
E=sum(abs( x0(1)*x+x0(2)-y ));
end

function E=fit3(x0,x,y)
E=sum(abs( x0(1)*x+x0(2)-y ).^2 );
end

注意在fminsearch函数的输入参数中有一个是[]，如果不加这个，会出现bug：

>> untitled
错误使用 fminsearch (第 106 行)}
第 3 个参数必须为 options 结构体。

出错 untitled (第 4 行)
p1=fminsearch(@fit1,[1 1],x,y);

详见fminsearch函数官方文档：

• fminsearch函数使用无导数法计算无约束的多变量函数的最小值，搜索由以下公式指定的问题的最小值：
  $$\underset{x}{\min }f(x)$$
  f(x)为返回标量的函数，x是向量或矩阵。

• x = fminsearch(fun,x0) 在点 x0 处开始并尝试求 fun 中描述的函数的局部最小值 x。

• x = fminsearch(fun,x0,options) 使用结构体 options 中指定的优化选项求最小值。使用 optimset 可设置这些选项。

• polyval函数
  
  y = polyval(p,x) 计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）：
  $$p(x)=p_1{x}^n+p_2{x}^{n-1}+\cdots +p_{n}x+p_{n+1}$$

• polyfit函数
  
  p = polyfit(x,y,n) 返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1：
  $$p(x)=p_1{x}^n+p_2{x}^{n-1}+\cdots +p_{n}x+p_{n+1}$$

• matlab拟合的通常方法

参见matlab数据拟合

• Code 4.7中interp2函数的使用：

Vq = interp2(X,Y,V,Xq,Yq) 使用线性插值返回双变量函数在特定查询点的插入值。结果始终穿过函数的原始采样。X 和 Y 包含样本点的坐标。V 包含各样本点处的对应函数值。Xq 和 Yq 包含查询点的坐标。

另外运行代码之后有bug：

>> untitled
函数或变量 ‘jj’ 无法识别。

出错 untitled (第 1 行)
x(1)=x0(jj); y(1)=y0(jj);

猜测可能是在画Figure 4.3 (b)时候的程序参数。

• 为何需要插值，以及如何使用数据求梯度？P127

• 相关资料：

“Nonlinear dynamics and chaos” by Steven H. Strogatz

4.2 Nonlinear Regression and Gradient Descent

• gradient $\nabla f(\mathbf{x})$

平常所谓“梯度”是指一个空间位置函数的变化率，在 数学上就是它的微商。对于多元函数，它对每个空间坐标变量都有一个偏微商，如$\frac{\partial \Phi }{x}，\frac{\partial \Phi }{y}，\frac{\partial \Phi }{z}$等。这些偏微商表示标量场$\Phi (x,y,z)$沿三个坐标方向的变化率。如果要问$\Phi (x,y,z)$沿任意方向$\Delta l$的变化率是多少呢？

如上图所示，$P$是标量场中的某个点，设此点标量场的数值是$\Phi (P)$，由$P$点引一个位移矢量$\Delta l$，到达附近的另一点$Q$，设$Q$点标量场的数值为$\Phi (Q)=\Phi (P)+\Delta \Phi$，令$Q$点向$P$点趋近，$\Delta l\to 0$，则标量场沿$\Delta l$方向的变化率为：

$$\frac{\partial \Phi }{\partial l}=\underset{\Delta l\to 0}{\lim }\frac{\Delta \Phi }{\Delta l}$$

$\frac{\partial \Phi }{\partial l}$叫做标量场$\Phi$在$P$点沿$\Delta l$方向的方向微商。

显然，在同一地点$P$，$\Phi$沿不同方向的方向微商一般说来是不同的。那么沿哪个方向的方向微商最大呢？前图中做通过$P$、$Q$两点$\Phi$的等值面，在两等值面上标量场的数值分别是$\Phi (P)$和$\Phi (P)+\Delta \Phi$。在局部范围看来，两等值面近似平行。通过$P$点引等值面的法线与另一等值面交于$Q^{\prime }$点。法线方向的位移矢量$\Delta \mathbf{n}=\overrightarrow{P{Q}^{\prime }}$是两等值面间最短的位移矢量，其它方向的位移矢量都比$\Delta \mathbf{n}$长。例如对于上述位移矢量$\Delta \mathbf{l}$，设它与$\Delta \mathbf{n}$的夹角为$\theta$，则由图中不难看出：

$$\Delta \mathbf{n}=\Delta \mathbf{l}\cos \theta ⩽\Delta \mathbf{l},or,\Delta \mathbf{l}=\frac{\Delta \mathbf{n}}{\cos \theta }⩾\Delta \mathbf{n}$$

沿$\Delta \mathbf{n}$方向的方向微商为：

$$\frac{\partial \Phi }{\partial n}=\underset{\Delta n\to 0}{\lim }\frac{\Delta \Phi }{\Delta n}=\underset{\Delta l\to 0}{\lim }\frac{\Delta \Phi }{\Delta l}\frac{1}{\cos \theta }=\frac{\partial \Phi }{\partial l}\frac{1}{\cos \theta }⩾\frac{\partial \Phi }{\partial l}$$

由此可见，沿$\Delta \mathbf{n}$方向的方向微商比任何其它方向的方向微商都大。

标量场的梯度定义为这样一个矢量：它沿方向微商最大的方向（即$\Delta \mathbf{n}$方向），数值上等于这个最大的方向微商（即$\frac{\partial \Phi }{\partial \mathbf{n}}$）。标量场$\Phi$的梯度通常记为$\mathrm{grad}\Phi$或$\nabla \Phi$。根据上面的分析可知，$\Phi$的梯度方向总是与$\Phi$的等值面垂直的，且指向变大的方向。

注意梯度指向的方向与朝向最小值的方向相反。

标量场的梯度是个矢量场。其在各种坐标系中具体表示式如下：

1）直角坐标系

$$\nabla \Phi =\frac{\partial \Phi }{\partial x}\mathbf{i}+\frac{\partial \Phi }{\partial y}\mathbf{j}+\frac{\partial \Phi }{\partial z}\mathbf{k}$$

2）柱坐标系

$$\nabla \Phi =\frac{\partial \Phi }{\partial \rho }{\mathbf{e}}_{\rho }+\frac{1}{\rho }\frac{\partial \Phi }{\partial \phi }{\mathbf{e}}_{\phi }+\frac{\partial \Phi }{\partial z}{\mathbf{e}}_z$$

3）球坐标系

$$\nabla \Phi =\frac{\partial \Phi }{\partial r}{\mathbf{e}}_r+\frac{1}{r}\frac{\partial \Phi }{\partial \theta }{\mathbf{e}}_{\theta }+\frac{1}{r\sin \theta }\frac{\partial \Phi }{\partial \phi }{\mathbf{e}}_{\phi }$$

值得注意的是前图也可以画成下面形状：

在这个图中我们可以发现有$\Delta \mathbf{l}$比$\Delta \mathbf{n}$短，即$\frac{\partial \Phi }{\partial n}<\frac{\partial \Phi }{\partial l}$。其实应当强调梯度是一个局域的概念，当我们把$P$点邻域放大，把$\Delta \Phi$减小，最终我们会得到如前图比较平坦的形状，进而保持梯度为最大方向微商的性质。这个与后面讨论步长$\delta$有关，因为离散值象征着讨论内容不再是局域的，使用的时候就需要小心。

• 为何梯度不与等高线垂直？

按照前面讨论的结论：$\Phi$的梯度方向总是与$\Phi$的等值面垂直的。并且（4.31）式可以转换为：

$${\mathbf{x}}_{k+1}-{\mathbf{x}}_k=-\delta \nabla f({\mathbf{x}}_k)$$

即连线方向与梯度方向平行。

但是在文章Figure 4.6右图中，其线并不与等高线垂直，这是怎么回事？注意这里说的等高线是计算梯度时的某点处的等高线，即迭代过程中的起始点。

另外注意到（4.33）式：

$$\frac{\partial F}{\partial \delta }=-\nabla f({\mathbf{x}}_{k+1})\nabla f({\mathbf{x}}_k)=0\text{\hspace{1em}(4.33)}$$

这个式子说明了上一点的梯度与下一点的梯度垂直。而两点之间的连线又是与梯度平行，所以上下两步的连线也应该互相垂直。所以文中图画的可能有点问题，应该画的更精确一些。

其中（4.33）式是将（4.31）代入到（4.32），并做微分得到：

$${\mathbf{x}}_{k+1}(\delta )={\mathbf{x}}_k-\delta \nabla f({\mathbf{x}}_k)\text{\hspace{1em}(4.31)}$$
$$F(\delta )=f({\mathbf{x}}_{k+1}(\delta ))\text{\hspace{1em}(4.32)}$$

$$\begin{array}{rl}\frac{\partial F}{\partial \delta }& =\frac{\partial f({\mathbf{x}}_{k+1}(\delta ))}{\partial \delta }\\ & =\frac{\partial f({\mathbf{x}}_{k+1})}{\partial {\mathbf{x}}_{k+1}}\frac{\partial ({\mathbf{x}}_{k+1}(\delta ))}{\partial \delta }\\ & =\frac{\partial f({\mathbf{x}}_{k+1})}{\partial {\mathbf{x}}_{k+1}}\frac{\partial ({\mathbf{x}}_k-\delta \nabla f({\mathbf{x}}_k))}{\partial \delta }\\ & =-\nabla f({\mathbf{x}}_{k+1})\nabla f({\mathbf{x}}_k)\\ & =0\end{array}$$

上面的推导虽然差不多，但是不严格，主要是因为文中使用的$\frac{\partial F}{\partial \delta }$会误导读者。其实上述是用了全微分的公式：

$$dz=\frac{\partial z}{\partial x}dx+\frac{\partial z}{\partial y}dy$$

或者：

$$\frac{dz}{dt}=\frac{\partial z}{\partial x}\frac{dx}{dt}+\frac{\partial z}{\partial y}\frac{dy}{dt}$$

假设${\mathbf{x}}_{k+1}=(x_{k+1}^1,x_{k+1}^2,x_{k+1}^3)$应用全微分公式可以将推导过程重新写为：

$$\begin{array}{rl}\frac{dF}{d\delta }& =\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^1}\frac{d(x_{k+1}^1(\delta ))}{d\delta }+\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^2}\frac{d(x_{k+1}^2(\delta ))}{d\delta }+\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^3}\frac{d(x_{k+1}^3(\delta ))}{d\delta }\\ & =\left(\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^1},\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^2},\frac{\partial f({\mathbf{x}}_{k+1})}{\partial x_{k+1}^3}\right)\cdot \left(\frac{d(x_{k+1}^1(\delta ))}{d\delta },\frac{d(x_{k+1}^2(\delta ))}{d\delta },\frac{d(x_{k+1}^3(\delta ))}{d\delta }\right)\\ & =\nabla f({\mathbf{x}}_{k+1})\frac{d({\mathbf{x}}_{k+1})}{d\delta }\\ & =\nabla f({\mathbf{x}}_{k+1})\frac{d({\mathbf{x}}_k-\delta \nabla f({\mathbf{x}}_k))}{d\delta }\\ & =-\nabla f({\mathbf{x}}_{k+1})\nabla f({\mathbf{x}}_k)\\ & =0\end{array}$$

值得进一步强调的是，上述这种确定步长和轨迹方向的方法极大提高了运行效率。

• 如果等高线为正圆形状，则轨迹会如何迭代？

设待优化函数：

$$f(x,y)=x^2+y^2$$

则：

$${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\delta \nabla f({\mathbf{x}}_k)=(1-2\delta )x\hat{\mathbf{x}}+(1-2\delta )y\hat{\mathbf{y}}$$

则：

$$F(\delta )=f({\mathbf{x}}_{k+1}(\delta ))=(1-2\delta {)}^2{x}^2+(1-2\delta {)}^2{y}^2$$

则其导数：

$$F^{\prime }(\delta )=-4(1-2\delta )x^2-4(1-2\delta )y^2$$

令$F^{\prime }(\delta )=0$：

$$-4(1-2\delta )x^2-4(1-2\delta )y^2=0$$

无解，可见文中步长确定的方法有它的适用范围。在这种情况下可以取变步长（与梯度下降速度有关），但是不用要求$F^{\prime }(\delta )=0$，虽然这样计算速度不如前者，但是也可以得到最终解。

实际上，只要初始点与最优化点的连线与所有等值线或者等值面垂直，就会有可能发生上述的情况。

• $F^{\prime }(\delta )=0$如何保证求的是最小值呢？是否有其他局域值？可见需要提前根据目标函数进行确定。

• 为何（4.31）中要用到梯度？

$${\mathbf{x}}_{k+1}(\delta )={\mathbf{x}}_k-\delta \nabla f({\mathbf{x}}_k)\text{\hspace{1em}(4.31)}$$

梯度的首要作用的为下一步的迭代确定方向，因为其指向值变大的方向。另外注意到梯度是一个向量，这个向量的模并不是单位1的，那么其值的大小是否对迭代有影响？如果$\delta$固定，则肯定会有影响，而$\delta$的存在就是调节这种影响到最优的一个参数。

但是如果只是上述这种看法是狭隘的，应该说对于下一步$f({\mathbf{x}}_{k+1})$的值同时受制于$\delta$,${\mathbf{x}}_k=(x_k^1,x_k^2,x_k^3)$，$\nabla f({\mathbf{x}}_k)=(\frac{\partial f}{\partial x_k^1},\frac{\partial f}{\partial x_k^2},\frac{\partial f}{\partial x_k^3})=(f_1^{\prime },f_2^{\prime },f_3^{\prime })$，也就是说下一步的值可以重写为：

$$f({\mathbf{x}}_{k+1})=f(x_k^1,x_k^2,x_k^3,f_1^{\prime },f_2^{\prime },f_3^{\prime },\delta )$$

广义来看，其实这就是一个运用上一步预测下一步的过程，但是这个预测必定是有限的，因为它一并没有具备所有导数的完备信息，二是不可能实现利用局域的信息来完全预测出下一步即为极值点的效果。

那么现在的问题是，公式（4.31）是否是最佳的解决方案呢？还是说这个公式也是经验性的一种形式？

• Alternating Descent

文中最后说：

Moreover, the method is derivative free, which
is attractive in many applications.
此外，该方法是无导数的，在许多应用中很有吸引力。

为何是没有导数呢？重新审视一下算法过程：

The basic strategy is simple: optimize along one variable at a time, seeking the minimum while holding all other variables fixed. After passing through each variable once, the process is repeated until a desired convergence is reached.

注意上述加粗的部分，这个单变量寻找最小值的过程其实文中没有具体给出，而且并没有说明这个单变量的最小值是否是全局的还是局域的。

但是由于此处取最小值对应的只是变化某个变量的最小值，而其他变量都是固定的，所以可以直接使用变化的变量的偏微分为0求出极值即可。比如对于待优化函数（4.29）式而言，假设初始点为(3,2)，首先做$x$的偏导数使其为0，并固定$y=2$，可以得到：

$$\begin{gathered} \frac{\partial f}{\partial x}=2x=0 \\ \Downarrow \\ x=0 \end{gathered}$$

所以经过一次迭代之后的坐标为(0,2)，由于此处处理的是凸函数，所以确定偏微分就是取极小值，而没有再做讨论。然后可以固定$x$，变化$y$，继续把坐标迭代下去。

• 注意Figure 4.5(b)的标签$\frac{\partial f}{\partial y}$是与图例中的说法是不一致的，如果图例是正确的，则(b)的标签应该为$\frac{\partial f}{\partial x}$。

• 全局最优算法

遗传算法，PSO（粒子群\鸟群算法），蚁群算法，MCMC，模拟退火算法，禁忌算法，人工鱼群算法，混合蛙跳算法，量子遗传算法，混沌优化算法，人工免疫算法，细菌觅食算法，人工萤火虫算法，蝙蝠算法，果蝇优化算法，生物地理优化算法，入侵野草优化算法，引力搜索算法，竞选算法，人工植物优化算法，文化算法，和声搜索算法，灰狼优化算法，布谷鸟搜索算法，化学反应优化算法等等

• 优化算法工具箱

PlatEMO-matlab

PlatEMO-github

• 不同error metrics的计算量

通过上述的了解才可以知道对于不同errror metrics（即（4.8））计算量来源于优化，确切的说来自于在找到最优解之前的迭代轨迹。在Chapter 3中P89中也有类似的讨论，其中提到：

$$\hat{\mathrm{s}}=\underset{\mathrm{s}}{\mathrm{argmin}}||\mathrm{s}|{|}_0,\mathrm{subject}\mathrm{to}\mathrm{y}=\mathrm{C}\mathrm{\Psi s}\text{\hspace{1em}(3.4)}$$
The optimization in (3.4) is non-convex, and in general the solution can only be found with a brute-force search that is combinatorial in n and K. In particular, all possible K-sparse vectors in Rn must be checked; if the exact level of sparsity K is unknown, the search is even broader. Because this search is combinatorial, solving (3.4) is intractable for even moderately large n and K, and the prospect of solving larger problems does not improve with Moore’s law of exponentially increasing computational power.

上面加粗的地方即表示类似于公式（4.8a）：

$$E_{\infty }(f)=\underset{1<k<n}{\max }|f(x_k)-y_k|,MaximumError({\ell }_{\infty })\text{\hspace{1em}(4.8a)}$$

这两者相同的地方就是要搜遍全部的可能组合，才能最终确定最优值，而在用最小二乘的error matrics下（凸函数）下，可以使用其函数本身的性质，比如使用梯度下降法来提高搜索效率，即不会遍历所有可能而只有一条搜索轨迹，其他最优算法也是这样。这也正体现了凸函数在优化中的重要性。而这也是作者一直强调的：

The choice of norm has a profound impact on the optimal solution achieved.

• 特殊情况

$$E_2(f)={\left(\frac{1}{n}\sum _{k=1}^n|f(x_k)-y_k{|}^2\right)}^{1/2}\text{\hspace{1em}(4.8c)}$$

当（4.8c）公式中平方下的函数$f(x_k)-y_k$长下面这个样子的时候情况会变得复杂：

在将其取模平方之后，负值会被折到正值上，此时极值点会形成极值环，此时会有无数个极值解，将会把问题变得复杂：

4.3 Regression and Ax = b: Over- and Under-Determined Systems

• 如何将回归问题对应到$\mathbf{Ax}=\mathbf{b}$

对于回归，我们可以用最一般的形式$\mathbf{Y}=f(\mathbf{X},\mathbf{\beta })$将其表示出来，其中$\mathbf{Y}$表示一组数据，$\mathbf{X}$表示另一组数据，$\mathbf{\beta }$为待定参数。那么当我们使用矩阵表示三者关系的时候是如何实现的？

文中的$\mathbf{Ax}=\mathbf{b}$，其中Loadings$\mathbf{x}$对应公式（4.4）$\mathbf{Y}=f(\mathbf{X},\mathbf{\beta })$中的待定参数$\mathbf{\beta }$。Outcomes$\text{b}$对应到$\text{Y}$，表示一组待拟合的数据。那么另一组待拟合的数据$\text{X}$对应矩阵中哪个元素呢？

注意到矩阵$\text{A}$的名称为Model terms，由于就剩下这么一个矩阵没有被定义，而公式（4.4）$\mathbf{Y}=f(\mathbf{X},\mathbf{\beta })$中尚有函数关系和数据$\text{X}$并没有对应到矩阵上去，所以可以想到的是这个Model terms矩阵$\text{A}$包含了两方面的信息：模型和数据，这一点在公式（4.44）中最明显，显然这个矩阵已经根据模型对数据进行了处理。

而这一点与传统的对线性方程组的表示略有不同，传统的表示中$\text{A}$为系数矩阵，$\mathbf{\beta }$为待解变量矩阵。而此处$\text{A}$为数据变量矩阵，$\mathbf{\beta }$为待解系数矩阵。

$$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_m\end{array}\right]$$

$$\left[\begin{array}{ccccc}|& |& |& \cdots & |\\ 1& x_j& x_j^2& \cdots & x_j^{n-1}\\ |& |& |& \cdots & |\end{array}\right]\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ ⋮\\ {\beta }_n\end{array}\right]=\left[\begin{array}{c}|\\ f(x_j)\\ |\end{array}\right]$$

这样在Figure 4.8中，$\mathbf{x}$对应传统线性方程组中的待求解变量，在有唯一值的时候，其行数等于$\mathbf{A}$的行数和列数，且等于$\mathbf{b}$的行数。但是此处并不只限定如此。

值得注意的是，拟合二次多项式函数，在最小二乘意义下的解，其可以表示为类似形式（4.13）：

$$\left(\begin{array}{cc}{\sum }_{k=1}^n{x}_k^2& {\sum }_{k=1}^n{x}_k\\ {\sum }_{k=1}^n{x}_k& n\end{array}\right)\left(\begin{array}{c}{\beta }_1\\ {\beta }_2\end{array}\right)=\left(\begin{array}{c}{\sum }_{k=1}^n{x}_k{y}_k\\ {\sum }_{k=1}^n{y}_k\end{array}\right)$$

能否推广到多次多项式？

在这种多项式模型的特殊情况下我们可以看到，其实一个欠定系统的最小二乘解，可以表示为一个具有唯一解的系统（线性方程组）。

• $||\text{Ax}-\text{b}|{|}_2$中随机矩阵$\mathbf{A}$有何意义？

其中的$\mathbf{A}$其实并没有模型附加在上面，单纯的就是个矩阵而已，那么列举这样一个矩阵有什么意义呢？如果我们反过来看，当我们定义一个模型矩阵$\text{A}$，其即使形式多样，代入数据之后，最后得到其实是一个每个元素均为数值的矩阵，看到这个全是数值的矩阵我们并不会联想到某个具体的模型概念。所以这里的随机矩阵代表着随意某一类模型代入数据之后的最终表示。而文中着重研究的并不是模型，而是${\ell }_1$和${\ell }_2$范数对最终解的影响。所以随机撒随机数并不影响最终结果，反而能够一定程度上体现一般性质（随机的无限可能性）。

• 伪逆与（4.40）比较

• 注意（4.41）式中没有了$||\text{Ax}-\text{b}|{|}_2$，但实质上式默认了其值为0，因为每个解都必须是其矩阵的解，即$||\text{Ax}-\text{b}|{|}_2=0$。

• 欠定与超定解的区别

• 一般情况下给定的数据的数量要远远大于系统的参数数量，也就是说，一般的模型都是超定的。

4.4 Optimization as the Cornerstone of Regression

• Figure 4.14中的box plots（箱线图\箱型图）

Figure 4.14中的Box-plot（箱形图）中的蓝色长方框（箱）内的数据是将所有数据从小到大排列，取其中数据序号（注意不是数值）在25%~75% （分别叫做下四分位数和上四分位数）上的数据，即将数据四等分，只取其中的两部分，上下方的长方框黑线顶端的黑横线是表示数据的最大值最小值。由于是按数据的个数进行分别，所以最终箱子上下的黑线长度是不同的，箱中的红色减号为中位数（横线），当然其也不一定在中间。

注意到红色加号和红色减号，红色加号应该为异常值，且可以有多个。

其中判断是否为异常值的标准为：

$$\mathbf{data}>Q3+1.5\times IQR,or,\mathbf{data}>Q1-1.5\times IQR$$

其中IQR为四分位距（interquartile range），又称为四分差，与方差、标准差一样，表示统计资料中各变量分散情形，但四分差更多为一种稳健统计（robust statistic）。计算公式为：

$$IQR=Q3-Q1$$

其中$Q3$，$Q1$分别为上四分位数（第三四分位数）和下四分位数（第一四分位数）。

有时还会在箱内加上其他符号，如下图：

• Figure 4.14 中的error

由于优化的过程需要一个停止的信号，这个信号一般就是error函数的值，所以不管哪些算法，只要参数够多，最终的error都可以是一样的（当然也有压不下去的情况，即有欠拟合的情况或者error不可能无穷小，而且这里也不是一个欠定的问题），那么这里的error是指什么呢？

如果一个算法对100组随机数据分别进行了拟合，得到了100组$\mathbf{\beta }$，然后为了展示算法的抗噪声干扰能力，可以将拟合曲线和$f(x)=x^2$作差，然后把残差通过箱型图展示出来。但是值得注意的是对比Figure 4.13 可以看出纵坐标轴其实是在一个量级上下的，再考虑到（c）为作者认为最好的结果，其二次项的系数为1，可以推测纵坐标就是$\mathbf{\beta }$值。所以这里纵坐标标注的有问题，更准确应该是$\mathbf{\beta }$或者是$Loadings$。当然后面可以通过Code 4.13判定。

• Figure 4.15中的Error

这里的error也比较奇怪，因为其数量级与Figure4.14中各方法和使用least-square regression的差很多。按照文章的意思这里是（a）比较不同方法的error；（b）比较不同p下的error（注意这里横坐标的虽然与Figure 4.14相同，但是这里的p是指最大的p，即分几次用不同长度的级数去拟合得到的总error）。这里有三种可能：

1. error代表的是另一个拟合对象的$\mathbf{\beta }$的误差的平均值，之所以是平均值，因为$\mathbf{\beta }$有好几个；
2. error代表的是长时间迭代优化后，最后再也压不下去的error函数值，考虑这一点主要是因为其值比较小；
3. 用一个算法对100组随机数据分别进行了拟合，得到了100组$\mathbf{\beta }$，然后为了展示算法的抗噪声干扰能力，可以将拟合曲线和$f(x)=x^2$作差，然后把残差通过箱型图展示出来。

再注意到Code 4.12中的代码：

n=100; L=4;
x=linspace(0,L,n);
f=(x.^2).’; % parabola with 100 data points

M=20; % polynomial degree
for j=1:M
phi(:,j)=(x.’).^(j-1); % build matrix A
end

for j=1:4
fn=(x.^2+0.1*randn(1,n)).’;
an=pinv(phi)fn; fna=phian; % least-square fit
En=norm(f-fna)/norm(f);
subplot(4,2,4+j),bar(an)
end

其中第一个加粗的部分为含噪声的函数值，其被用作训练的数据；
第二个加粗的部分先使用求伪逆函数pinv求出$\text{Ax}=\text{b}$中$\text{A}$的伪逆，然后根据公式（1.20）：
$${\mathbf{A}}^{†}\mathbf{A}\tilde{\mathbf{x}}={\mathbf{A}}^{†}\mathbf{b}⟹\tilde{\mathbf{x}}=\tilde{\mathbf{V}}{\tilde{\mathbf{\Sigma }}}^{-1}{\tilde{\mathbf{U}}}^{\ast }\mathbf{b}$$

最终利用伪逆得到等价意义上的最小二乘解（详情见Chapter 1.4）。

En值为通过利用含噪声的数据得出的参数反演得到拟合曲线，然后与$f(x)=x^2$作差。所以确定是上面第三种error。

• Figure 4.15（c）

注意Figure 4.15（c）中的误差为同区间误差，是拟合曲线与目标曲线之间的差。在下一节中，作者说在高阶的时候也可以体现过拟合的现象。而在下一节，作者又强调过拟合可以使是拟合曲线在test set（测试区间）表现的特征，而非training set。所以过拟合的标准很受人为定义的影响。下一节的Fgiure 4.17总结了上述两种不同的误差大小（被分别称为Interpolation Error和Extrapolation Error）比如我们用带噪声的数据进行拟合，但是我们不觉得$f(x)=x^2$是我们的期待的函数，而取其他函数形式，则最终误差可能会很大。所以确定哪个是我们期待的函数这个问题就有很大的操作空间。

4.5 The Pareto Front and Lex Parsimoniae

由于一般拟合的数据量都很大，超定形式的系统是最常见的，参数少的情况，容易出现而欠定的系统容易出现过拟合，当然两者都会出现。

一定要注意拟合曲线，并不代表最好的效果是穿过过数据中的每个点，而不穿过每个点并不意味着欠拟合，欠拟合和过拟合并不是严格对立的概念，需要有阈值进行区别。当然由于一般的超定系统也不会有穿过每个点的解。

• Figure 4.17

图例中：

In the interpolation regime x ∈[0, 4], the model error stays constrained, with increasing error due to overfitting for polynomials of degree greater than 2.

加粗部分应该是指各个模型最终截止的拟合误差标准。但是文中并没有讨论这个标准，因为这个标准直接影响到两种Error的大小。

而这个training set如果是在数据上均匀随机取点呢？效果如何？

• 过拟合

• 注意欠定和超定都是的可以是欠拟合、过拟合或者正好合适的，但是欠拟合更容易过拟合，因为参数比较多，但是并不等价。由于判断过拟合还是欠拟合需要一个最终的心中期待的目标函数，所以判断过拟合和欠拟合的标准会随需要与主观观念的不同而不同。

• 另外注意拟合的好并不意味着非得要穿过所有数据点。判断过拟合还是欠拟合比较微妙，需要谨慎对待。

4.6 Model Selection: Cross-Validation

• Code 4.15

n=200; L=8;
x=linspace(0,L,n);
x1=x(1:100); % train
x2=x(101:200); % test
n1=length(x1);
n2=length(x2);
ftrain=(x1.^2).’; % train parabola x=[0,4]
ftest=(x2.^2).’; % test parbola x=[4,5]
figure(1), subplot(3,1,1),
plot(x1,ftrain,’r’,x2,ftest,’b’,’Linewidth’,[2])

其中加粗的部分中转置的符号为.'，这与一般的不同。matlab中矩阵后面加单引号是共轭转置，加点和单引号是转置：

• interpolation error

为何这里要用插值误差这个术语呢？百科中对插值的解释：

• 在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。
• 插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。
• 插值：用来填充图像变换时像素之间的空隙。

由此可见插值概念与拟合曲线的概念很接近。文中与interpolation相对的是extrapolation，所以这里的插值应该是针对拟合数据区间内的时候使用的点，得到拟合参数之后再利用这些离散点还原（呈现）已拟合曲线的过程被称为插值，而由此造成的误差叫做interpolation error。

另外在wiki中的词条Regression analysis：

• Interpolation and extrapolation
  Regression models predict a value of the Y variable given known values of the X variables. Prediction within the range of values in the dataset used for model-fitting is known informally as interpolation. Prediction outside this range of the data is known as extrapolation. Performing extrapolation relies strongly on the regression assumptions. The further the extrapolation goes outside the data, the more room there is for the model to fail due to differences between the assumptions and the sample data or the true values.

上述说法与前面分析的差不多，还是指对训练或者拟合的数据点所在的区间内的曲线预测。

• Figure 4.19

在除掉小于threshold的部分之后，是否应当对留下来的做一个值上的补偿？

• withheld data

As for partitioning the data, a common strategy is
to break the data into 70% training data, 20% validation data, and 10% withheld data.

指的是出training data和validation data之外的无用数据？

其中data set，training set，test set与data，training data，validation data之间不等价？

4.7 Model Selection: Information Criteria

• Code 4.19

其中函数trapz的官方使用手册说明：

• Q = trapz(Y) 通过梯形法计算 Y 的近似积分（采用单位间距）。Y 的大小确定求积分所沿用的维度：
  
  如果 Y 为向量，则 trapz(Y) 是 Y 的近似积分；
  
  如果 Y 为矩阵，则 trapz(Y) 对每列求积分并返回积分值的行向量；
  
  如果 Y 为多维数组，则 trapz(Y) 对其大小不等于 1 的第一个维度求积分。该维度的大小变为 1，而其他维度的大小保持不变。

• Q = trapz(X,Y) 根据 X 指定的坐标或标量间距对 Y 进行积分。
  
  如果 X 是坐标向量，则 length(X) 必须等于 Y 的大小不等于 1 的第一个维度的大小。
  
  如果 X 是标量间距，则 trapz(X,Y) 等于 X*trapz(Y)。

• 注意在没有truth model的情况下，$f(x)$由truth data代替。所以增加了KL divergence的使用难度。

• maximum likelihood

最大似然估计（极大似然估计，似然指可能性，为文言文）是建立在这样的思想上：已知某个参数能使这个样本出现的概率最大，我们当然不会再去选择其他小概率的样本，所以干脆就把这个参数作为估计的真实值。

设X为离散型随机变量，$\theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$为多维参数向量，如果随机变量相互独立，则可得概率函数$P(X_1=x_1,\dots ,X_n=x_n)={\prod }_{i=1}^{n}p(x_1;{\theta }_1,{\theta }_2,\dots ,{\theta }_k)$，在$\theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$固定时，上式表示$X_1=x_1,\dots ,X_n=x_n$的概率；当$X_1=x_1,\dots ,X_n=x_n$已知的时候，它又变成$\theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$的函数，可以把它记为$L({\theta }_1,{\theta }_2,\dots ,{\theta }_k)={\prod }_{i=1}^{n}p(x_1;{\theta }_1,{\theta }_2,\dots ,{\theta }_k)$，称此函数为似然函数。似然函数值的大小意味着该样本值出现的可能性的大小，既然已经得到了样本值$X_1=x_1,\dots ,X_n=x_n$，那么它出现的可能性应该是较大的，即似然函数的值也应该是比较大的，因而最大似然估计就是选择使$L({\theta }_1,{\theta }_2,\dots ,{\theta }_k)$达到最大值的那个$\theta$作为真实的估计。

• log-likelihood

• ground truth model

wiki对于ground truth的解释

• Ground truth is information that is known to be real or true, provided by direct observation and measurement (i.e. empirical evidence) as opposed to information provided by inference.
• Statistics and machine learning
  “Ground truth” may be seen as a conceptual term relative to the knowledge of the truth concerning a specific question. It is the ideal expected result.[2] This is used in statistical models to prove or disprove research hypotheses. The term “ground truthing” refers to the process of gathering the proper objective (provable) data for this test. Compare with gold standard. For example, suppose we are testing a stereo vision system to see how well it can estimate 3D positions. The “ground truth” might be the positions given by a laser rangefinder which is known to be much more accurate than the camera system.
  
  Bayesian spam filtering is a common example of supervised learning. In this system, the algorithm is manually taught the differences between spam and non-spam. This depends on the ground truth of the messages used to train the algorithm – inaccuracies in the ground truth will correlate to inaccuracies in the resulting spam/non-spam verdicts.

知乎中的解释：

• 机器学习包括有监督学习(supervised learning)，无监督学习(unsupervised learning)，和半监督学习（semi-supervised learning）.在*有监督学习中，数据是有标注的，以(x, t)的形式出现，其中x是输入数据，t是标注.正确的t标注是ground truth， 错误的标记则不是。（也有人将所有标注数据都叫做ground truth）
• Ground truth是摄影、测量与遥感学领域常用词汇，其解释就是字面意思：地面真值，地面实况；延伸到图像处理、机器学习等其他领域一般表示真实值，正确答案（或正确测量数据）。它是一个正确的基准值，一般用来进行误差估算和效果评价。例子可参考高赞答案。
• 不能单从字面意思理解，这个词是从遥感、气象领域borrow过来的，在这些领域里按照字面意思就很好理解了，就是指卫星在对地=地面进行观测时，可能会有误差，所以这时候就要用Ground Truth 来修正观测的数据。理解了它在原领域中的意思，那么它在ML中的含义自然也就明了了。就是指监督学习中的label，而算法的预测label就相当于卫星的观测数据。

• Autoregressive model (AR)

在统计学、计量经济学和信号处理学中，自回归(AR)模型是一种随机过程的表示;因此，它被用来描述自然界、经济学等中的某些时变过程。自回归模型指定输出变量线性地依赖于它自己以前的值和一个随机项（一个不完全可预测的项);因此，该模型是一个随机差分方程(或递归关系，不应与微分方程混淆）的形式。它与移动平均（moving-average，MA）模型一起，是随机结构更为复杂的时间序列的一般自回归移动平均（autoregressive–moving-average，ARMA）模型和自回归综合移动平均（autoregressive integrated moving average，ARIMA）模型的特例和关键组成部分;它也是向量自回归模型（vector autoregressive，VAR）的一个特例，它由一个以上的互锁随机差分方程（interlocking stochastic difference equation ）在多个演化随机变量中组成的系统。
p阶的自回归模型$AR(p)$定义为：
$$X_t=c+\sum _{i=1}^p{\phi }_i{X}_{t-i}+{\epsilon }_t$$
其中${\phi }_1,\dots ,{\phi }_p$是模型中的系数，$c$是常数，${\epsilon }_t$是白噪声。如果使用backshift operator$B$来表示：
$$X_t=c+\sum _{i=1}^p{\phi }_i{B}^i{X}_t+{\epsilon }_t$$
其中操作符backshift operator$B$在时间序列$X=X_1,X_2,\dots$中，定义为：
$$B{X}_t=X_{t-1}forallt>1$$
如果使用lag opoerator$L$，则同样可以表示为：
$$L{X}_t=X_{t-1}forallt>1$$
如果回退多次则：
$$L^k{X}_t=X_{t-k}forallt>k$$

• Moving-average model (MA)

在时间序列分析中，移动平均模型(MA模型)，又称移动平均过程，是一种常用的单变量时间序列建模方法。移动平均模型规定，输出变量线性依赖于一个随机(不完全可预测)项的当前值和各种过去值。移动平均模型与自回归(AR)模型一样，是随机结构更为复杂的时间序列一般ARMA和ARIMA模型的特殊情况和关键组成部分。移动平均模型（Moving-average model）不应该与移动平均（moving average）混淆，这是一个不同的概念，尽管有一些相似之处。与AR模型相反，有限MA模型总是平稳的。
$MA(q)$用来表示q阶移动平均模型，定义为：
$$X_t=\mu +{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q{\epsilon }_{t-q}=\mu +\sum _{i=1}^q{\theta }_i{\epsilon }_{t-i}+{\epsilon }_t$$
其中$\mu$为序列的平均值，${\theta }_1,\dots ,{\theta }_q$是模型中的参数，${\epsilon }_t,{\epsilon }_{t-1},\dots ,{\epsilon }_{t-q}$为白噪声误差项。上式也可以重新写为：
$$X_t=\mu +(1+{\theta }_{1}B+\cdots +{\theta }_q{B}^q){\epsilon }_t$$
因此，移动平均模型在概念上是该系列的当前值与当前和以前（观察到的）白噪声误差项或随机冲击的线性回归。假设每个点的随机冲击（random shocks）是相互独立的，来自相同的分布，通常是一个正态分布，位置在零和常数尺度。

• Autoregressive integrated moving average (ARIMA)

在统计学和计量经济学中，特别是在时间序列分析中，自回归综合移动平均（ARIMA）模型是自回归移动平均（ARMA）模型的推广。这两种模型都适合于时间序列数据，以更好地理解数据或预测序列中的未来点（预测）。ARIMA模型应用在某些情况下,数据表明non-stationarity的意思(但不是方差/自协方差)，在一个初始差分步骤（对应于模型的“集成”部分）可以应用一个或多个时间消除的non-stationarity意味着函数（例如，趋势）。
给定实数时间序列$X_t$，则ARMA(p’,q)模型定义为：
$$X_t-{\alpha }_1{X}_{t-1}-\cdots -{\alpha }_{p^{\prime }}{X}_{t-p^{\prime }}={\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q{\epsilon }_{t-q}$$
还可以表示为：
$$(1-\sum _{i=1}^{p^{\prime }}{\alpha }_i{L}^i)X_t=(1+\sum _{i=1}^q{\theta }_i{L}^i){\epsilon }_t$$
其中${\alpha }_i$为自回归部分的参数，${\theta }_i$为移动平均部分的参数，${\epsilon }_t$为误差项，通常被假定为独立的、从均值为零的正态分布抽样的同分布变量。
如果$(1-{\sum }_{i=1}^{p^{\prime }}{\alpha }_i{L}^i)$有d个单位重根（unit root，因子为(1-L)）上式可转化为：
$$(1-\sum _{i=1}^{p^{\prime }}{\alpha }_i{L}^i)X_t=(1-\sum _{i=1}^{p^{\prime }-d}{\phi }_i{L}^i)(1-L{)}^d$$
一个ARIMA(p,d,q)过程可以用$p=p^{\prime }-d$表示多项式分解的性质，定义为：
$$(1-\sum _{i=1}^p{\phi }_i{L}^i)(1-L{)}^d{X}_t=(1+\sum _{i=1}^q{\theta }_i{L}^i){\epsilon }_t$$
因此可以认为是ARMA(p+d,q)过程具有具有d个单位根的自回归多项式的特殊情况。(因此，由 d > 0 的 ARIMA 模型准确描述的过程都不是广义平稳的)

上述情况可以概括如下。

$$(1-\sum _{i=1}^p{\phi }_i{L}^i)(1-L{)}^d{X}_t=\delta +(1+\sum _{i=1}^q{\theta }_i{L}^i){\epsilon }_t$$
上式用平移项（drift）$\frac{\delta }{1-\sum {\phi }_i}$表示了ARIMA(p,d,q)过程。

• Computing AIC and BIC Scores

注意此处的训练数据是autoregressive model产生的，而拟合用的模型是三个不同时间延迟（对应$L$操作符）的ARIMA模型。

• Code 4.20

T = 100; % Sample size
DGP = arima(’Constant’,-4,’AR’,[0.2,0.5],’Variance’,2);
y = simulate(DGP,T);

EstMdl1 = arima(’ARLags’,1);
EstMdl2 = arima(’ARLags’,1:2);
EstMdl3 = arima(’ARLags’,1:3);

logL = zeros(3,1); % Preallocate loglikelihood vector
[_,,logL(1)] = estimate(EstMdl1,y,’print’,false);
[_,,logL(2)] = estimate(EstMdl2,y,’print’,false);
[_,,logL(3)] = estimate(EstMdl3,y,’print’,false);

[aic,bic] = aicbic(logL, [3; 4; 5], T*ones(3,1))

• 相关资料：

Regression Towards Mediocrity in Hereditary Stature

如何解释$\lceil$线性回归$\rfloor$的含义？

“最优化方法及其MATLAB实现” by许国根

机器学习里经常出现ground truth这个词，能否准确解释一下？

宋亚磊 21:23
是的是的，拟合是这样的，这体现在优化时无法选择或不必选择0伪范数