其实所谓的多变量的线性回归（Linear Regression with multiple variables ）本质上将与单变量的线性回归没啥差别。因此我们完全可以用上一节中的梯度下降算法来解决，只需要在每一次迭代的时候多考虑几个变量而已。所以这一节就稍微介绍一下了，不再用例子分析。不过毕竟多了一些变量，在对多变量跑梯度下降算法时，显然对参数的调节就更加重要了，因此我们首先得学会一些参数调节的技巧。这些技巧在实际的操作过程中尤为重要。
调节参数
特征值缩放（Feature Scaling）
我们知道在进行梯度下降算法时，经常会有这样一个问题，就是如果两个特征变量x_i和x_j的范围差别很大，那么每一次迭代时，他们下降的相对速率就不一样，这就直接导致了范围大的那个变量下降的过慢。为了解决这个问题，我们引入了特征值缩放这个技巧。其实也很简单，就是把每一个参数的变化范围大概控制在一起，这样他们的下降速率就不会受到相互的牵连。具体的操作就是将他所有参数的值，都缩放在大小差不多的一个区域内（通常是(-1,1)）。为了进行这样的缩放，我们就需要接下来的方法。
标准化（Mean Normalization）
为了缩放到一个区域（(-1,1)）内，我们要做两件事：第一，把平均值对应到0点；第二，把最大最小值对应到-1和1这两点。想想也很好理解，用公式表示就是：
这里的x_i表示需要标准化的数据，u表示在数据集中这种参数的平均值，sigma表示数据集中这种参数的标准差。看着很眼瘦？没错，这个就是我们求正态分布的那个标准化转换函数。经过这样的转换，每一个变量都会将他的值缩放到（-1,1）中了，进而方便我们进行梯度下降。
学习因子$alpha$的选取
我们知道学习因子既不能过大也不能过小。实际上，对于一个学习因子是好是坏我们可以这样判断。对于一个特定的学习因子alpha，我们可以画出J函数与迭代次数的关系图。对于一个好的alpha，他的函数图像应该是一个类似y=frac1x的图像。如果J函数下降的速度过慢，显然我们需要增加学习因子使他下降的更快；如果J函数是增函数，或者在某些区间内有上升的趋势，那么显然我们需要减小学习因子防止学习过度。
多项式回归（Polynomial Regression ）
对于某些不能用线性回归的问题，我们有时候可以试着用多项式来进行回归拟合。其实多项式回归完全可以看成是多变量的线性回归问题，因为我们完全可以把其中的x^i看成是第i个独立的变量，只不过他的值是由x推出来的而已。原理很简单，但是如果想不到那就头大了0.0。
公式法（Normal equation）
介绍
对于多变量的线性回归，除了用我们之前学的GD算法，我们其实还有另外一个直接套公式的算法（卧槽早说）。这牵涉到线性代数的知识，我们需要做的只是将数据集组合成几个矩阵，然后运算一个公式即可，这个公式就叫 Normal equation （觉得叫成“正规方程”好难听）：
其实这个方法就是矩阵最小二乘法了。
其中X^T表示X的转置，X^{-1}表示X的逆。
这个公式，得到的结果就是对整个数据集最好的线性逼近。下面介绍下这个公式怎么用。
用法
假设我们的数据集有m个数据，我们的变量有n个（即每一个数据有n个特征值），那么我们可以把所有的输入数据写成一个m*（n 1）的矩阵。为什么是n 1列，因为我们知道，对n个参数的回归需要用到n 1个theta。只是这当中的theta_0一直取1罢了。因此这个矩阵的第一列全部都是1。
上面的这个矩阵，就是我们公式中的X了。
对于公式中的y，其实是由所有的输出数据组成的列向量，与输入一一对应。
最后，我们只需要计算上面这个公式的结果就能计算出theta向量的值了，（很明显这个向量应该是一个n 1维的列向量）
至于这个公式的证明。。。貌似很麻烦，百度normal equation 可以找到需要的证明过程。
说明
这里有一个问题，就是如果X^TX没有逆元怎么办?
事实上正常情况下是不会没有逆元的，没有逆元的情况一般有两种：
1、参数中有线性相关的两组量。比如一个参数是用厘米表示的高度，另一个参数是用米表示的同样的那个高度。显然，这样选择参数是没有意义的，对于这种情况我们只要删掉其中一个参数就好了。
2、数据集的数量m小于参数的个数n。通常情况下当然是不会发生这种事的，毕竟都没几个数据还怎么好意思来拟合？不过万一出现这种状况，也是有对应的办法的，这种情况以后再说。
性能
我们注意到这里需要求矩阵的逆，而X^TX是一个n*n的矩阵，因此他的复杂度通常至少是O(n^3)（虽然可以优化，但是也差不多这么大）。这在n比较小的时候效果还挺不错的，而且比GD算法更加方便准确。但是当参数的个数变多了的时候，效率就会特别慢，相比之下GD算法的优势就显现出来了。
实现
接下来我就结合上一节的例子试着实现一下，虽然是单变量，但是道理是一样的。
代码语言：javascript复制import numpy as np
import theano.tensor as T
from pylab import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#read file
file=open("data","r")
brain=[]
body=[]
for line in file:
num,_brain,_body=line.split()
brain.append(_brain)
body.append(_body)
#calculate
size=len(brain)
X=np.matrix(np.vstack((np.ones(size),brain))).reshape(2,size).T.astype(np.float)
Y=np.matrix(body).T.astype(np.float)
out=(X.T*X).I*X.T*Y
print out
#draw
out=np.array(out)
px=np.linspace(-100,600,100)
py=px*out[1] out[0]
plt.plot(px,py)
plt.scatter(brain,body,s=100,alpha=.5)
plt.xlim(-100,600)
plt.ylim(-200,1400)
plt.show()（这里写代码的时候要非常注意np.matrix的用法以及强制指定数据类型）
最后算得：theta_0=55.95185719,theta_1=1.24271276，从图上看还是拟合的挺不错的。
