地址：https://github.com/Raymond-Yang-2001/AndrewNg-Machine-Learing-Homework
单变量线性回归找到一维方程，拟合一条直线。
h
w
,
b
(
x
)
=
b
+
w
x
h_{w,b}(x)=b+wx
hw,b(x)=b+wx
w
w
w和
b
b
b是参数，为了方便运算，可以给
x
x
x加上一个
x
0
=
1
x_0=1
x0=1
h
w
,
b
(
x
)
=
b
x
0
+
w
x
1
h_{w,b}(x)=bx_{0}+wx_{1}
hw,b(x)=bx0+wx1
J
(
w
,
b
)
=
1
2
m
∑
i
=
1
m
(
h
w
,
b
(
x
(
i
)
)
−
y
(
i
)
)
2
J(w,b)=\frac{1}{2m}\sum_{i=1}^{m}(h_{w,b}(x^{(i)})-y^{(i)})^{2}
J(w,b)=2m1i=1∑m(hw,b(x(i))−y(i))2
为了避免不恰当的数据范围带来损失过大或过小，（例如若数据数值过大，损失可能会在
1
0
5
10^5
105或者
1
0
6
10^6
106这个数量级，不适合直观分析）在评估损失的时候，可以对
h
w
,
b
(
x
(
i
)
)
h_{w,b}(x^{(i)})
hw,b(x(i))和
y
(
i
)
y^{(i)}
y(i)先进行标准化，使得损失数值在可评估的范围内。但在进行梯度下降时，不进行此操作
w
j
=
w
j
−
α
∂
∂
w
j
J
(
w
,
b
)
=
w
j
−
α
1
m
∑
i
=
1
m
(
h
w
,
b
(
x
(
i
)
)
−
y
(
i
)
)
x
(
i
)
w_j=w_{j}-\alpha\frac{\partial}{\partial{w_j}}{J(w,b)}=w_{j}-\alpha \frac{1}{m}\sum_{i=1}^{m}{(h_{w,b}(x^{(i)})-y^{(i)})x^{(i)}}
wj=wj−α∂wj∂J(w,b)=wj−αm1i=1∑m(hw,b(x(i))−y(i))x(i)
b
j
=
b
j
−
α
∂
∂
b
j
J
(
w
,
b
)
=
w
j
−
α
1
m
∑
i
=
1
m
(
h
w
,
b
(
x
(
i
)
)
−
y
(
i
)
)
b_j=b_{j}-\alpha\frac{\partial}{\partial{b_j}}{J(w,b)}=w_{j}-\alpha \frac{1}{m}\sum_{i=1}^{m}{(h_{w,b}(x^{(i)})-y^{(i)})}
bj=bj−α∂bj∂J(w,b)=wj−αm1i=1∑m(hw,b(x(i))−y(i))
这里，我们可以使用
θ
\theta
θ统一标识参数，包括
w
w
w和
b
b
b。
即，第
j
j
j个参数
θ
j
\theta_j
θj的更新可以写为：
θ
j
=
θ
j
−
α
∂
∂
w
j
J
(
θ
;
x
)
=
w
j
−
α
1
m
∑
i
=
1
m
(
h
θ
(
x
(
i
)
)
−
y
(
i
)
)
x
(
i
)
\theta_{j}=\theta_{j}-\alpha\frac{\partial}{\partial{w_j}}{J(\theta;\mathbf{x})}=w_{j}-\alpha \frac{1}{m}\sum_{i=1}^{m}{(h_{\theta}(x^{(i)})-y^{(i)})x^{(i)}}
θj=θj−α∂wj∂J(θ;x)=wj−αm1i=1∑m(hθ(x(i))−y(i))x(i)
其中
α
\alpha
α是学习率。
多变量线性回归试图找出多个变量和预测值之间的关系。例如，房子大小、房子卧室数量和房价之间的关系。
在样本的不同特征数值差异过大的时候，基于梯度的优化方法会出现一些问题。例如存在如下回归方程：
h
θ
(
x
)
=
θ
0
+
θ
1
x
1
+
θ
2
x
2
h_{\theta}(x)=\theta_{0}+\theta_{1}x_{1}+\theta_{2}x_{2}
hθ(x)=θ0+θ1x1+θ2x2
假设
x
2
x_{2}
x2的范围是
0
∼
1
0\sim1
0∼1，
x
1
x_1
x1的范围是
1
0
3
∼
1
0
4
10^3\sim10^4
103∼104。我们根据梯度同时优化
θ
0
∼
θ
2
\theta_0\sim\theta_2
θ0∼θ2，使得其均改变了相同的大小，那么显然在输入样本相同的情况下，
θ
1
\theta_1
θ1的变动会比
θ
2
\theta_2
θ2导致更大的输出的变化。这也可以理解为模型对
θ
1
\theta_1
θ1更敏感。如下损失等线图所示，
θ
1
\theta_1
θ1的微小变动会带来损失的剧烈变化。这种情况下，参数的优化会更加困难。
解决这个问题的方法之一就是特征缩放，将两个特征缩放到相同的范围内。例如，可以进行z-score标准化工作：
x
n
e
w
=
x
−
μ
σ
x_{new} = \frac{x-\mu}{\sigma}
xnew=σx−μ
其中，
μ
\mu
μ是数据集的均值，
σ
\sigma
σ是标准差，新数据的分布是均值为0，标准差为1的分布。
数据标准化后的参数等损失图如下所示：
由于对数据进行了缩放，所以最后得到的参数也会出现相应的缩放。其具体关系如下：
θ
0
+
θ
1
∼
d
+
1
x
1
∼
d
+
1
−
μ
x
σ
x
=
y
−
μ
y
σ
y
\theta_{0}+\theta_{1\sim d+1}\frac{x_{1\sim d+1}-\mu_{x}}{\sigma_{x}}=\frac{y-\mu_{y}}{\sigma_{y}}
θ0+θ1∼d+1σxx1∼d+1−μx=σyy−μy
这里我们对
y
y
y也进行了标准化，事实上也可以不这么做，对性能没有任何影响。但是对y的标准化使得参数变得更小，对于初始化为0的参数能更快达到收敛。
在标准化
y
y
y这种情况下，参数的逆缩放公式为：
θ
1
∼
d
+
1
n
e
w
=
θ
1
∼
d
+
1
σ
x
σ
y
\theta_{1\sim d+1}^{new}=\frac{\theta_{1\sim d+1}}{\sigma_{x}}\sigma_{y}
θ1∼d+1new=σxθ1∼d+1σy
得到：
θ
0
σ
y
+
θ
1
∼
d
+
1
n
e
w
(
x
1
∼
d
+
1
−
μ
x
)
=
y
−
μ
y
\theta_{0}\sigma_{y}+\theta_{1\sim d+1}^{new}(x_{1\sim d+1}-\mu_{x})=y-\mu_{y}
θ0σy+θ1∼d+1new(x1∼d+1−μx)=y−μy
θ
0
n
e
w
=
θ
0
σ
y
+
μ
y
−
θ
1
∼
d
+
1
n
e
w
μ
x
\theta_{0}^{new}=\theta_{0}\sigma_{y}+\mu_{y}-\theta_{1\sim d+1}^{new}\mu_{x}
θ0new=θ0σy+μy−θ1∼d+1newμx
其中，在向量化运算时，
θ
1
∼
d
+
1
n
e
w
\theta_{1\sim d+1}^{new}
θ1∼d+1new和
μ
x
\mu_{x}
μx均为(1,d)的向量，乘法应该采用向量内积。
设数据 x \boldsymbol{x} x的维度是 ( n , d ) (n,d) (n,d)，其中n是样本数量，d是样本特征的维度。为了计算方便，我们在样本上添加一个额外数值全为1的特征维度，使其维度变为 ( n , d + 1 ) (n,d+1) (n,d+1)
import numpy as np
def square_loss(pred, target):
"""
计算平方误差
:param pred: 预测
:param target: ground truth
:return: 损失序列
"""
return np.sum(np.power((pred - target), 2))
def compute_loss(pred, target):
"""
计算归一化平均损失
:param pred: 预测
:param target: ground truth
:return: 损失
"""
pred = (pred - pred.mean(axis=0)) / pred.std(axis=0)
target = (pred - target.mean(axis=0)) / target.std(axis=0)
loss = square_loss(pred, target)
return np.sum(loss) / (2 * pred.shape[0])
class LinearRegression:
"""
线性回归类
"""
def __init__(self, x, y, val_x, val_y, epoch=100, lr=0.1):
"""
初始化
:param x: 样本, (sample_number, dimension)
:param y: 标签, (sample_numer, 1)
:param epoch: 训练迭代次数
:param lr: 学习率
"""
self.theta = None
self.loss = []
self.val_loss = []
self.n = x.shape[0]
self.d = x.shape[1]
self.epoch = epoch
self.lr = lr
t = np.ones(shape=(self.n, 1))
self.x_std = x.std(axis=0)
self.x_mean = x.mean(axis=0)
self.y_mean = y.mean(axis=0)
self.y_std = y.std(axis=0)
x_norm = (x - self.x_mean) / self.x_std
y_norm = (y - self.y_mean) / self.y_std
self.y = y_norm
self.x = np.concatenate((t, x_norm), axis=1)
self.val_x = val_x
self.val_y = val_y
def init_theta(self):
"""
初始化参数
:return: theta (1, d+1)
"""
self.theta = np.zeros(shape=(1, self.d + 1))
def validation(self, x, y):
x = (x - x.mean(axis=0)) / x.std(axis=0)
y = (y - y.mean(axis=0)) / y.std(axis=0)
outputs = self.predict(x)
curr_loss = square_loss(outputs, y) / (2 * y.shape[0])
return curr_loss
def gradient_decent(self, pred):
"""
实现梯度下降求解
"""
# error (n,1)
error = pred - self.y
# gradient (d+1, 1)
gradient = np.matmul(self.x.T, error)
# gradient (1,d+1)
gradient = gradient.T / pred.shape[0]
# update parameters
self.theta = self.theta - (self.lr / self.n) * gradient
def train(self):
"""
训练线性回归
:return: 参数矩阵theta (1,d+1); 损失序列 loss
"""
self.init_theta()
for i in range(self.epoch):
# pred (1,n); theta (1,d+1); self.x.T (d+1, n)
pred = np.matmul(self.theta, self.x.T)
# pred (n,1)
pred = pred.T
curr_loss = square_loss(pred, self.y) / (2 * self.n)
val_loss = self.validation(self.val_x, self.val_y)
self.gradient_decent(pred)
self.val_loss.append(val_loss)
self.loss.append(curr_loss)
print("Epoch: {}/{}\tTrain Loss: {:.4f}\tVal loss: {:.4f}".format(i + 1, self.epoch, curr_loss, val_loss))
# un_scaling parameters
self.theta[0, 1:] = self.theta[0, 1:] / self.x_std.T * self.y_std[0]
self.theta[0, 0] = self.theta[0, 0] * self.y_std[0] + self.y_mean[0] - np.dot(self.theta[0, 1:], self.x_mean)
return self.theta, self.loss, self.val_loss
def predict(self, x):
"""
回归预测
:param x: 输入样本 (n,d)
:return: 预测结果 (n,1)
"""
# (d,1)
t = np.ones(shape=(x.shape[0], 1))
x = np.concatenate((t, x), axis=1)
pred = np.matmul(self.theta, x.T)
return pred.T
数据集可视化
训练集与测试集划分
from LinearRegression import LinearRegression
epochs = 200
alpha = 1
linear_reg = LinearRegression(x=train_x_ex,y=train_y_ex,val_x=val_x_ex, val_y=val_y_ex, lr=alpha,epoch=epochs)
start_time = time.time()
theta,loss, val_loss = linear_reg.train()
end_time = time.time()
Train Time: 0.0309s
Val Loss: 6.7951
训练过程可视化
与sk-learn比较预测曲线
数据可视化与训练集、验证集
from LinearRegression import LinearRegression
alpha = 0.1
epochs = 1000
multi_lr = LinearRegression(train_x,train_y_ex,val_x=val_x,val_y=val_y_ex, epoch=epochs,lr=alpha)
start_time = time.time()
theta, loss, val_loss = multi_lr.train()
end_time = time.time()
Train Time: 0.1209s
Val Loss: 4.187（采用归一化后数据计算损失）
训练过程可视化
预测平面（与sk-learn对比）
其中蓝色为本算法的预测平面，灰色为sk-learn的预测平面
对于线性回归算法的实现打到了较好的性能，可以尝试调节学习率或者迭代次数来获得更好的性能。由于采用了矩阵运算代替了循环，所以训练时间大大缩短，但仍未到达sk-learn库函数的水平。