pydata

Keep Looking, Don't Settle

Tensorflow简介--05: Multivariate Regression with Stochastic Gradient Descent

尽管TenforFlow(tf)最主要的功能是进行深度学习,包括卷积神经网络(CNN),循环神经网络(RNN)深度神经网络(DNN)等等。但是作为一个深度学习框架,其也提供了进行一般统计学习和优化的工具。前面(1, 2, 3, 4)我们已经介绍了tf的基本设置和各个命令的含义。下面我们使用一个例子来展示怎么使用tf来进行多元回归拟合数据。具体的说,用一个例子来说明怎么借助tfStochastic Gradient Descent做多项式回归。下面的动图展示了拟合函数的参数(回归系数)在tf经过多次迭代以后,拟合曲线渐渐逼近真实值的过程。

alt text

在使用tf的时候,我们需要定义下面几个东西:

1. 输入和输出数据
2. 定义参数
3. 定义损失函数,以及使用什么办法来优化(比如下面使用随机梯度下降的办法)

首先我们需要调用将要使用的一些packages:

import math,sys,os,numpy as np, pandas as pd, tensorflow as tf
from matplotlib import pyplot as plt, rcParams, animation, rc, style
from __future__ import print_function, division
from ipywidgets import interact, interactive, fixed
from ipywidgets.widgets import *

%precision 4
np.set_printoptions(precision=4, linewidth=100)

style.use('ggplot')
%matplotlib inline
plt.rcParams['figure.figsize'] = (12, 9)
rc('animation', html='html5')

我们的目的是使用tf来进行多元回归(这个例子是多项式)。首先我们需要模拟产生一个样本数据:\(X\)是0到9上的均匀分布的点,\(Y\)\(X\)的多项式函数, \(y = -50 + 30 \cdot x - 12 \cdot x^2 + x^3 + \epsilon\). \(X\)\(Y\)的关系见下图。

np.random.seed(99)
x = np.linspace(0, 9, 100)[:, np.newaxis]
y =  np.polyval([1, -12, 30, -50], x) + np.random.randn(100, 1)

plt.plot(x, y)
plt.xlabel('X')
h = plt.ylabel('Y')
h.set_rotation(0)
plt.title("$y=x^3 - 12 x^2 + 30 x - 50 + \epsilon$")

png

做多项式回归的时候选择基函数为\([1, \space x, \space x^2, \space x^3 ]\),python中可以使用np.power(x, range(4))得到。因为随机梯度下降的时候梯度跟\(X\)的值有关系,所以我们需要对\(X\)的每一列进行标准化让所有的\(X\)的值处在同样的尺度下。

数据标准化的一个办法是所有的值(我们的例子中是每一列数据)减去均值,除以方差

$$ X^{normalized} = \frac{X - \mu(X)}{\sigma(X)}$$

numpy对应的code是(x - x.mean(axis = 0)) / x.std(axis = 0).

另一个办法是每一列数据除以最大值。因为我们的\(X\)都大于0,所以这样做的结果就是把所有的\(X\)变换到\([0, 1]\)区间上。numpy对应的code是x / np.max(x, axis = 0).

pol_order = 4
x = np.power(x, range(pol_order))

def normalized_(x):
    return (x - x.mean(axis = 0)) / x.std(axis = 0)

def normalized(x):
    l = len(x)
    return np.concatenate((np.ones(l).reshape(l, 1), normalized_(x[:, 1:])), axis = 1)

x = normalized(x)

normalized过的\(X\)\(Y\)的关系为

png

把上面定义好的data分为training和test两部分:

np.random.seed(99)
order = np.random.permutation(len(x))
portion = 20
test_x = x[order[:portion]]
test_y = y[order[:portion]]
train_x = x[order[portion:]]
train_y = y[order[portion:]]

下面开始定义模型。tf通过数据流程图来定义模型,然后编译执行。

我们从定义输入(实际值)输出开始。如前所述,这是通过占位符(placehold)来定义.placehold有两个必须的参数,第一个是数据类型,通常为浮点类型(float).第二个是输入输出的数据的维度。假如我们的输入X是一个(m, n)维的矩阵,我们可以把第二个参数定义为[m, n]。我们用[None n]来表示待定的维度,这样我们可以使用batch来部分的输入数据,而不是全部m个数据。这样做的原因是如果m很大,机器的内存可能处理不了。这时候把数据分为很多batch,每次输入部分的数据(一个batch)来拟合,保证内存不会溢出。我们还可以添加第三个参数name,用来在TensorBoard中展示名字。

with tf.name_scope("Input/Output"):
    inputs = tf.placeholder(tf.float32, [None, pol_order], name="X")
    outputs = tf.placeholder(tf.float32, [None, 1], name="Yhat")

接下来定义模型表达式。模型的参数变量通过tf.Variable来定义,跟tf.placehold的定义很相似,不同的地方在于tf.Variable的值会随着拟合而发生变化。线性模型的表达式为\(Y=X \times W\)。所以参数W的维度与X的列数相等。矩阵X和W相乘由tf.matmul定义。

with tf.name_scope("Model"):
    W = tf.Variable(tf.zeros([pol_order, 1], dtype=tf.float32), name="W")
    y = tf.matmul(inputs, W)

最后定义损失函数和它的最优化函数。通常对线性回归模型我们都是要最小化均方误差。损失函数定义为平方损失\(L = \sum_{i=1}^{n}(y_i - f(x_{i}|W)) ^2\)。 把回归函数写成矩阵的形式为\(Y = X \cdot W\),那么损失函数也可以表达为矩阵相乘\(L = (Y - XW)'(Y-XW)\). 要最小化损失函数,我们首先求解出梯度为

$$ \begin{aligned} \frac{\partial L}{\partial W} & = \frac{\partial ( (Y - XW)'(Y-XW))}{\partial W} \\\\ & = \frac{\partial (Y'Y - 2W'X'Y + W'X'XW)}{\partial W} \\\\ &= -2X'(Y - XW) \end{aligned}$$

tf(包括sklearn)都是通过数值办法来最小化损失函数的,而不是通过解析解,尽管我们知道对线性模型回归参数\(W\)参在解析解\(\hat{W} = (X'X)^{-1}X'Y\). 具体就是让W每次沿着它的梯度方向变化一点点(learning rate),然后反复迭代,直到达到我们的收敛条件或者停止条件。在第一篇中我们介绍过梯度下降的原理:假如损失函数是一座山,我们现在站在山坡上朝四周看看,然后找一个方向能让我们走一点点路但是下降的最快,那个方向就是梯度的方向。

tf提供了很多最优化(不同的梯度下降)的办法,参见Optimizers.对某些数据使用梯度下降的时候会出现鞍点,导致GradientDescentOptimizer可能会找不到最优点或者只是local的最优点。这时应该尝试别的最优化的办法。在我们的例子中,GradientDescentOptimizer找不到最优解, 但是AdamOptimizer可以找到最优解。

learning rate的参数设置为trainable=False表示每次迭代的时候learning rate本身不变化。learning rate太小会导致迭代步骤太多,计算所需要的时间就会很长。learning rate太大也会有问题:在你将要靠近最优点的时候,因为下一步跨的太大,导致跨过了最优点(overshoot)。比较好的办法是:刚开始迭代的时候learning rate可以大一些,然后随着渐渐靠近最优点,learning rate应该渐渐减小。

with tf.name_scope("SGD"):
    learning_rate = tf.Variable(0.5, trainable=False)
    cost_op = tf.reduce_mean(tf.pow(y-outputs, 2))
    train_op = tf.train.AdamOptimizer(learning_rate).minimize(cost_op)

到此为止我们已经准备好了计算流程图,下面我们开始定义对话(session)来执行数据流程图。数据流程图中的每个节点(operation or op)在计算过程当中会获得tensor,每个tensor是一个多为数组。对话将图的op分到cpu或者gpu进行计算。计算完后以tensor返回,在python中其类型为np.ndarray

1. `sess.as_default()`定义了默认对话, `sess.run()` 和 `sess.eval()` 只能在这个对话下面执行
2. `sess.run()`如果执行`operation`的话(比如`cost_op`),需要提供输入值`feed_dict={inputs: train_x, outputs: train_y}`
3. 对tensor来说,`tensor.eval()` == `sess.run(tensor)`
tolerance = 1e-5

epochs = 1
last_cost = 0
alpha = 0.3
max_epochs = 50000

sess = tf.Session() # default session

# make default session, run() and eval() should be executed in this session
with sess.as_default(): 
    init = tf.initialize_all_variables()
    sess.run(init)
    sess.run(tf.assign(learning_rate, alpha))
    while True:
        sess.run(train_op, feed_dict={inputs: train_x, outputs: train_y})

        if epochs%100==0:
            cost = sess.run(cost_op, feed_dict={inputs: train_x, outputs: train_y})
            print ("Epoch: %d - Error: %.4f" %(epochs, cost))

            if abs(last_cost - cost) < tolerance or epochs > max_epochs:
                print ("Stop Criteria Reached. Break!")
                break
            last_cost = cost

        epochs += 1

    w = W.eval()  # fetch the value of tensor
    print ("w =", w)
    print ("Test Cost =", sess.run(cost_op, feed_dict={inputs: test_x, outputs: test_y}))
Epoch: 100 - Error: 1110.4116
Epoch: 200 - Error: 321.6223
Epoch: 300 - Error: 167.4508
Epoch: 400 - Error: 144.2539
Epoch: 500 - Error: 135.3392
Epoch: 600 - Error: 126.8727
Epoch: 700 - Error: 118.1265
Epoch: 800 - Error: 109.2254
Epoch: 900 - Error: 100.3067
Epoch: 1000 - Error: 91.4902
Epoch: 1100 - Error: 82.8799
Epoch: 1200 - Error: 74.5638
Epoch: 1300 - Error: 66.6154
Epoch: 1400 - Error: 59.0937
Epoch: 1500 - Error: 52.0441
Epoch: 1600 - Error: 45.4990
Epoch: 1700 - Error: 39.4793
Epoch: 1800 - Error: 33.9944
Epoch: 1900 - Error: 29.0440
Epoch: 2000 - Error: 24.6189
Epoch: 2100 - Error: 20.7023
Epoch: 2200 - Error: 17.2708
Epoch: 2300 - Error: 14.2958
Epoch: 2400 - Error: 11.7446
Epoch: 2500 - Error: 9.5815
Epoch: 2600 - Error: 7.7693
Epoch: 2700 - Error: 6.2698
Epoch: 2800 - Error: 5.0450
Epoch: 2900 - Error: 4.0584
Epoch: 3000 - Error: 3.2750
Epoch: 3100 - Error: 2.6623
Epoch: 3200 - Error: 2.1907
Epoch: 3300 - Error: 1.8338
Epoch: 3400 - Error: 1.5684
Epoch: 3500 - Error: 1.3748
Epoch: 3600 - Error: 1.2362
Epoch: 3700 - Error: 1.1391
Epoch: 3800 - Error: 1.0725
Epoch: 3900 - Error: 1.0279
Epoch: 4000 - Error: 0.9986
Epoch: 4100 - Error: 0.9800
Epoch: 4200 - Error: 0.9685
Epoch: 4300 - Error: 0.9615
Epoch: 4400 - Error: 0.9574
Epoch: 4500 - Error: 0.9551
Epoch: 4600 - Error: 0.9538
Epoch: 4700 - Error: 0.9531
Epoch: 4800 - Error: 0.9528
Epoch: 4900 - Error: 0.9526
Epoch: 5000 - Error: 0.9525
Epoch: 5100 - Error: 0.9525
Epoch: 5200 - Error: 0.9525
Epoch: 5300 - Error: 0.9525
Stop Criteria Reached. Break!
w = [[ -56.3677]
 [  78.7995]
 [-293.1844]
 [ 209.9559]]
Test Cost = 1.12171

最终,损失函数在training上面收敛到0.9525,循环达到了停止的条件。这时候损失函数在test上的值为1.12177

最后的拟合函数为

$$ y = -56.3677 + 78.8013 \cdot x - 293.1888 \cdot x^2 + 209.9586 \cdot x^3$$

用于拟合的数据的前5行为

y x_0 x_1 x_2 x_3
0 -32.780156 1.0 1.645531 2.074123 2.392560
1 -40.784224 1.0 -1.576245 -1.106391 -0.877834
2 -59.082880 1.0 -0.155892 -0.426134 -0.551511
3 -73.715893 1.0 0.155892 -0.124437 -0.313780
4 -32.524330 1.0 -0.848747 -0.900180 -0.822070

拟合函数的图像为

png

alt text

最后说明的是,我们用多元线性回归的例子来演示怎么使用tf. 对于这个问题,python的机器学习package sklearn 可以提供更快捷的结果。

from sklearn import linear_model
from sklearn.metrics import mean_squared_error

reg = linear_model.LinearRegression()
reg.fit(train_x, train_y)
print ("The regression coefficient is %s" %(reg.coef_))

print ("The MSE from sklearn is: %s" %(mean_squared_error(test_y, reg.predict(test_x))))
The regression coefficient is [[   0.       78.8177 -293.2297  209.9838]]
The MSE from sklearn is: 1.12222708205

Reference

  1. TensorFlow Mechanics 101
  2. Session management
  3. what is the difference between Session.run() and Tensor.eval()