最小二乘法的梯度梯度下降解法

考虑一个最简单的问题，假如构建了一个如图的散点，需要线性拟合。

首先，数据我们自己产生吧。

import numpy as np
import pylab as plt
x = np.linspace(0,20,40)          #0到20,40个点
y = 2*x+[np.random.randint(0,10) for i in range(40)]  #2×x加一个扰动

我们先用最小二乘法的标准公式，计算出w，公式见 百度百科。主要要计算矩阵的逆，代码如下：

A = np.mat(np.zeros((20,1)))        
B = np.mat(np.zeros((20,1)))
for i in range(20):
    A[i]=np.mat([x[i]])       #矩阵.这里考虑y=kx，因此没有考虑b。如果考虑则np.mat([x[i],1]
    B[i]=np.mat([y[i]])    #y矩阵
AT = np.transpose(A)       #转置
c = np.linalg.inv(AT*A)   #逆
w = c*AT*B          #套用百度百科公式
print(w)

得到结果：

[[ 2.62526316]]

接下来 我们用梯度下降法来计算w的值。相关理论见 维基百科-Solution of a linear system。
主要就是迭代，设置一个初始数，然后减去梯度导数得到新值，再进去迭代。。。

i = 0
start = -2
xishu = 0.001
while i<100:
    new = start - xishu*2*AT*(A*start-B)
    start = new
    i+=1
print(start)

得到的结果是：

[[ 2.62526316]]

这里需要注意的一点，迭代的系数要足够小，才能收敛，如果系数比较大，很容易发散，感兴趣可以自己试试。

缺点：

（1）系数的影响很大，需要很小的系数

（2）收敛速度慢，这种原始的代码，有些情况迭代1万多次情况才稍微好些。相较于直接计算得出来，差距很大。

y=kx+bz这种形式的，稍微改一下。

import numpy as np
import pylab as plt

x = np.linspace(0,20,100)
y = 2*x+[np.random.randint(0,1) for i in range(100)]

A = np.mat(np.zeros((100,2)))
B = np.mat(np.zeros((100,1)))
for i in range(100):
    A[i]=np.mat([x[i],1])
    B[i]=np.mat([y[i]])
AT = np.transpose(A)


i = 0
start = np.mat([[1],[1]])
xishu = 0.00001
while i<5000:
    new = start - xishu*2*AT*(A*start-B)

    start = new
    i+=1
print(start)

c = np.linalg.inv(AT*A)
w = c*AT*B
print(w)


plt.plot(x,y,'ro')
plt.show()