AYu

参考

一、实验目的

掌握最小二乘法求解（无惩罚项的损失函数）、掌握加惩罚项（2范数）的损失函数优化、梯度下降法、共轭梯度法、理解过拟合、克服过拟合的方法(如加惩罚项、增加样本)。

二、实验要求及实验环境

要求：

生成数据，加入噪声；
用高阶多项式函数拟合曲线；
用解析解求解两种loss的最优解（无正则项和有正则项）
优化方法求解最优解（梯度下降，共轭梯度）；
用你得到的实验数据，解释过拟合。
用不同数据量，不同超参数，不同的多项式阶数，比较实验效果。

实验环境：

Google colab：python

三、设计思想（本程序中的用到的主要算法及数据结构）

数据结构：数组/向量/矩阵

算法：

解析法求最优解

无正则项时，通过对目标函数

求导可以得到

加入正则项，则目标函数变为

解为

梯度下降法

1）确定当前位置的损失函数的梯度，对于θi,其梯度表达式如下：

2）用步长乘以损失函数的梯度，得到当前位置下降的距离，即：

3）确定是否所有的θi，梯度下降的距离都小于ε，如果小于ε则算法终止（或者指定迭代次数），当前所有的θi(i=0,1,…n)即为最终结果。否则进入步骤4.

4）更新所有的θ，对于θi，其更新表达式如下。更新完毕后继续转入步骤1.

矩阵表达方式：

共轭梯度法

在回归模型中，回归系数β正是线性方程组 Ax=b 的解，其中 A=X′X，b=X′y。共轭梯度法，就是想像上面这个式子一样，把解x表达成共轭向量基的线性组合：只要依次算出所有的共轭向量pi和对应的系数αi，就可以得出x。而在实际情况中，有可能用更少数目的pi就能得到对x的良好近似。

四、实验结果与分析

当数据量小（=10）时：

不带正则项的解析解得出的参数预测结果在M=9时出现了严重的过拟合现象：

随着M的增大，在训练集上的误差越来越小，但在测试集上的误差在M>6后开始增大：

当加上正则项后，M=9时也没有出现过拟合现象：

超参数λ的选择会影响正则化的效果，当λ太小和接近1时训练集的误差都很大，在训练集上的误差由于控制了参数的大小也随着增大：

使用梯度下降法，学习率取太大会使M=9的模型参数变得太大而无法计算，而太小会影响计算速度，最终取α= 0.01，迭代次数为1000000，结果如下：

使用共轭梯度法时，设定误差阈值为1e-10，在M=9时也会出现过拟合现象：

当数据量由10增加至30时，过拟合现象明显减缓.

解析解不带正则项方法：

对应Erms图：

在测试集上的误差明显减小。

在解析解带正则项方法中表现也有所改进：

共轭梯度法M=9的模型过拟合现象也减轻了：

五、结论

解析解方法能够得到回归的精确解，在特征量合适时可以使用。
梯度下降法能够通过控制迭代次数控制解的优化程度，适合于不太需要精确解的场合。
共轭梯度法既可以求出精确解，也可以在中途停下，同时因为减少了曲折的迭代路线比梯度下降法效率更高，是较优的选择。
当出现过拟合现象时，原因可能是数据量不够或者模型设定得过于复杂，可以通过增加正则项、增大数据量或者减少模型复杂度解决。

六、源代码（带注释）

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def generate_data_1(size):  #随机生成数据点
  x = np.random.rand(size)
  y = np.sin(2 * np.pi * x) + np.random.normal(0,0.3,size)
  return x, y

def generate_data_2(size):  #等距生成数据点
  x = np.linspace(0,1,size)
  y = np.sin(2 * np.pi * x) + np.random.normal(0,0.3,size)
  return x, y

size = 30
x_train, y_train = generate_data_2(size)
x_func = np.linspace(0,1,100)
y_func = np.sin(2 * np.pi * x_func)
plt.plot(x_train,y_train,'.b')
plt.plot(x_func,y_func,'-g')
plt.show()


def polynomial_X(x, exp):   #构造输入特征矩阵
  if(exp == 0):
    return np.ones(x.size)
  X = np.vstack((np.ones(x.size),x))   #将X的第一列设为全1
  for i in range(2, exp+1):
    X = np.vstack((X, np.power(x, i)))
  return X.T


def train_analytic(x_train, exp):  #不带正则项解析法求参
  X = polynomial_X(x_train, exp)
  if(exp == 0):      # M=0的情况特殊处理，不能使用pinv函数
    theta = 1/(X.T.dot(X)) * X.T.dot(y_train)
  else:
    theta = np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y_train)
  return theta

def predict_analytic(x_train, x_test, exp):  #不带正则项解析法预测
  theta = train_analytic(x_train, exp)
  return polynomial_X(x_test, exp).dot(theta)

j = 1
for i in [0, 1, 3, 9]:
  y_predict = predict_analytic(x_train, x_func, i)
  plt.subplot(2,2,j)
  j += 1
  plt.plot(x_func, y_predict, '-r')
  plt.plot(x_train,y_train,'.b')
  plt.plot(x_func,y_func,'-g')
plt.show()


def calculate_erms(y, t):   #计算Erms
  return np.sqrt(np.mean(np.square(y-t)))

x = np.array([0,1,2,3,4,5,6,7,8,9])  #Erms图的x轴，代表degree
training_erms = np.zeros(10)
test_erms = np.zeros(10)

theta = []
for i in range(10):
  theta.append(train_analytic(x_train, i))
  y_train_predict = polynomial_X(x_train, i).dot(theta[i])
  training_erms[i] = calculate_erms(y_train_predict, y_train)
  y_test_predict = predict_analytic(x_train, x_func, i)
  test_erms[i] = calculate_erms(y_test_predict, 
y_func+np.random.normal(0,0.3,len(y_func)))  #真实数据有高斯分布的误差

plt.plot(x, training_erms, 'o-b', label="Training")
plt.plot(x, test_erms, 'o-r', label="Test")
plt.legend()
plt.xlabel("degree")
plt.ylabel("Erms")
plt.show()


def train_analytic_with_regularization(x_train, exp, lamb):   #带正则项解析法求参
  X = polynomial_X(x_train, exp)
  if(exp == 0):
    theta = 1/(X.T.dot(X)+lamb) * X.T.dot(y_train)
  else:
    theta = np.linalg.pinv(X.T.dot(X)+lamb*np.eye(X.shape[1])).dot(X.T).dot(y_train)
  return theta

def predict_analytic_with_regularization(x_train, x_test, exp, lamb):  #带正则项解析法预测
  theta = train_analytic_with_regularization(x_train, exp, lamb)
  return polynomial_X(x_test, exp).dot(theta)

j = 1
lamb = 0.001
for i in [0, 1, 3, 9]:
  y_predict = predict_analytic_with_regularization(x_train,x_func,i,lamb)
  plt.subplot(2,2,j)
  j += 1
  plt.plot(x_func,y_predict, '-r')
  plt.plot(x_train,y_train,'.b')
  plt.plot(x_func,y_func,'-g')
plt.show()


x = np.linspace(-40, 0)   #Erms图x轴，代表lnλ的值
training_erms = np.zeros_like(x)
test_erms = np.zeros_like(x)

X = polynomial_X(x_train, 9)  #选择M=9观察正则项效果
theta = []
for i,index in zip(x, range(x.size)):
  theta.append(train_analytic_with_regularization(x_train, 9, np.exp(i)))
  y_train_predict = X.dot(theta[index])
  training_erms[index] = calculate_erms(y_train_predict, y_train)
  y_test_predict = predict_analytic_with_regularization(x_train, x_func, 9, np.exp(i))
  test_erms[index] = calculate_erms(y_test_predict, 
y_func+np.random.normal(0,0.3,len(y_func)))

plt.plot(x, training_erms, 'o-b', label="Training")
plt.plot(x, test_erms, 'o-r', label="Test")
plt.legend()
plt.xlabel("lnλ")
plt.ylabel("Erms")
plt.show()


def gradient_descent(X, y, alpha, iterations):  #梯度下降法求参
  if(np.size(X.shape) == 1):  #若X只有一列，代表只有一个参数
    theta = 0
  else:
    theta = np.zeros(np.size(X, 1))
  for i in range(iterations):
    theta = theta - alpha * X.T.dot(X.dot(theta) - y)
  return theta

def gradient_descent_predict(X, y_train, X_test, exp, alpha, iterations):
  theta = gradient_descent(X, y_train, alpha, iterations)
  if(exp == 9):
    print("M=9时，theta=", theta)  # 查看M=9时得到的theta值
  return X_test.dot(theta)

j = 1
iterations = 1000000
# X = polynomial_X(x_train, 5)
# X_test = polynomial_X(x_func, 5)
# for alpha in [0.01, 0.03, 0.1, 0.3]:
#   plt.subplot(2,2,j)
#   j += 1
#   y_predict = gradient_descent_predict(X, y_train, X_test, 5, alpha, iterations)
#   plt.plot(x_func,y_predict,'-r')
#   plt.plot(x_train,y_train,'.b')
#   plt.plot(x_func,y_func,'-g')
# plt.show()

for i in [0, 1, 3, 9]:
  plt.subplot(2,2,j)
  j += 1
  X = polynomial_X(x_train, i)
  X_test = polynomial_X(x_func, i)
  y_predict = gradient_descent_predict(X, y_train, X_test, i, 0.01, iterations)
  plt.plot(x_func,y_predict,'-r')
  plt.plot(x_train,y_train,'.b')
  plt.plot(x_func,y_func,'-g')
plt.show()

# X = polynomial_X(x_train, 9)
# X_test = polynomial_X(x_func, 9)
# y_predict = gradient_descent_predict(X, y_train, X_test, 9, 0.0001, 1000000000)
# plt.plot(x_func,y_predict,'-r')
# plt.plot(x_train,y_train,'.b')
# plt.plot(x_func,y_func,'-g')
# plt.show()


def conjugate_gradient(X, y):   #共轭梯度法求参
  esp = 1e-10          #当误差小于该值时停止迭代
  A = X.T.dot(X)
  b = X.T.dot(y)
  if (np.size(X.shape) == 1):
    theta = 0
    r = b - A*theta
    p = r
    r2 = r*r
    alpha = r2 / (p * A * p)
    theta = theta + alpha * p
  else:
    theta = np.zeros(np.size(X, 1))
    r = b - A.dot(theta)
    p = r
    r2 = r.T.dot(r)
    err = 1
    # for i in range(np.size(X, 1)):
    while err > esp:
      alpha = r2 / (p.T.dot(A).dot(p))
      theta = theta + alpha * p
      r = r - alpha * A.dot(p)
      r2_new = r.T.dot(r)
      err = np.sqrt(r2_new)
      beta = r2_new / r2
      p = r + beta * p
      r2 = r2_new
  return theta

def conjugate_gradient_predict(X_train, y_train, X_test):
  theta = conjugate_gradient(X_train, y_train)
  # if(exp == 9):
  #   print("M=9时，theta=", theta)  # 查看M=9时得到的theta值
  return X_test.dot(theta)

j = 1
for i in [0, 1, 3, 9]:
  plt.subplot(2,2,j)
  j += 1
  y_predict = conjugate_gradient_predict(polynomial_X(x_train, i), 
y_train, polynomial_X(x_func, i))
  plt.plot(x_func,y_predict,'-r')
  plt.plot(x_train,y_train,'.b')
  plt.plot(x_func,y_func,'-g')
plt.show()

参考