机器学习-->矢量化运算，矩阵运算，广播

在机器学习中，我们常常需要处理各式各样的数据，其中向量，矩阵形式的数据是经常遇见的。同时又经常需要对两个或者多个向量，矩阵形式的数据做元素级的处理。通常最简单无脑的办法就是循环的对数据一个个的进行处理。显然这种处理方式不仅费时费力而且写出来的代码给人感觉很繁琐。

既然数据呈向量或者矩阵形式，为什么我们不把这种形式加以利用呢？要知道矢量化的运算要比等价的纯python快上一两个数量级甚至更多。

首先来讲讲矢量运算。

矢量化运算函数

元素级数组函数

可以对ndarray（即向量形式）中的数据执行元素级运算，其结果是返回一个新的向量。

一元函数（参数只有一个向量）:

注：上面图片内函数是对向量中每一个元素进行操作，返回一个新的向量。

二元函数（参数为两个向量）：

注：上面图片内的函数可以对两个向量同时执行元素级操作，返回一个新的向量。

数学和统计方法

可以通过数组上的一组数学函数对整个数组（即矩阵形式）或者某个　轴向的数组（即向量形式）进行聚合计算。其结果返回一个标量或者是一个向量。

from numpy.random import randn
import numpy as np

arr=randn(5,4)
print "arr.mean:",arr.mean()##统计整个矩阵的平均值
##arr.mean: -0.432992161059 返回结果是一个标量

print "arr.mean(1):",arr.mean(1)##1表示行向，即统计每行的平均值
## arr.mean(1): [-1.02563423 -0.6250071  -0.40480331  0.28313662  0.4170792 ]返回结果是一个向量

print "arr.mean(0):",arr.mean(0)##1表示列向，即统计每列的平均值
## arr.mean(1): [-1.02563423 -0.6250071  -0.40480331  0.28313662  0.4170792 ]返回结果是一个向量

注：对矩阵而言返回标量或一个新的向量。

逻辑运算

注：返回一个新的向量。

矩阵运算

矩阵加法

这里我们假设有两个numpy.array（向量或是矩阵）数据

向量加法除了上面说到的np.add(A,B)，还可以直接，同样的表示两向量对应位相加得到一个新的矩阵。

import numpy as np

A=np.array([[1,2],[3,4]])
B=np.array([[2,3],[4,5]])

print A+B
'''
[[3 5]
 [7 9]]
'''

矩阵减法

与上面的矩阵加法同理，只不过是两矩阵对应元素作减法操作得到一个新的矩阵。

A=np.array([[1,2],[3,4]])
B=np.array([[2,3],[4,5]])

print A-B
'''
[[-1 -1]
 [-1 -1]]
'''

矩阵乘法

向量的乘法稍微麻烦一点。如果两矩阵是numpy.array类型。那么 表示两矩阵对应元素相乘得到一个新的矩阵。

而 表示线性代数上矩阵相乘，其结果返回一个矩阵（若是两向量相乘则返回标量）。注意这里是矩阵B的转置形式，这样矩阵相乘格式上才能对应上，否则会报错。

A=np.array([[1,2],[3,4]])
B=np.array([[2,3],[4,5]])

print A*B
'''
[[ 2  6]
 [12 20]]
'''

print A.dot(B)
'''
[[10 13]
 [22 29]]
'''

如果你觉得用dot表示矩阵乘法不太直观，那么你可以将矩阵由numpy.array类型转成numpy.matrix类型。

A=np.matrix(A)
B=np.matrix(B)
c=A*B

这个时候直接用点乘就可以了。

矩阵除法

跟上面矩阵点乘法类似，只不过改为对应元素作除法

import numpy as np

A=np.array([[1,2],[3,4]])
B=np.array([[2,3],[4,5]])

print (A/B.astype(float))##astype将矩阵内每个元素转换类型
'''
[[ 0.5         0.66666667]
 [ 0.75        0.8       ]]
'''


print A.dot(1/B.astype(float))
'''
[[ 1.          0.73333333]
 [ 2.5         1.8       ]]
'''

矩阵求导

向量偏导公式：

标量对向量的导数

A为n×n的矩阵， x为n×1的列向量，记

那么可得：

标量对方阵的导数

A为n×n的矩阵， |A|为A的行列式。

广播

广播是矢量化运算里面一个非常重要但却又非常难以理解的一种对数据的操作方式。

将标量与数组（向量）合并时发生最简单的广播：

import numpy as np

arr=np.arange(5)
print "arr:",arr ##arr: [0 1 2 3 4]
print "arr*4:",arr*4 ##arr*4: [ 0  4  8 12 16]

显然标量4被广播到向量中每一个元素上。

import numpy as np
from numpy.random import randn

np.random.seed(1)
arr=randn(4,3)
print arr
'''
[[ 1.62434536 -0.61175641 -0.52817175]
 [-1.07296862  0.86540763 -2.3015387 ]
 [ 1.74481176 -0.7612069   0.3190391 ]
 [-0.24937038  1.46210794 -2.06014071]]
'''

print "arr.mean(0):",arr.mean(0) ##沿着行方向计算每列的均值
print "arr.mean(0).reshape(1,3):",arr.mean(0).reshape(1,3) ##在列方向计算每列的均值
'''
arr.mean(0): [ 0.0923349   0.42184748  0.85644227]
arr.mean(0).reshape(1,3): [[ 0.51170453  0.23863806 -1.14270302]]
'''

demeaned=arr-arr.mean(0)##每个元素减去这个元素所在列的均值,去均值化

print demeaned
'''
[[ 1.11264083 -0.85039448  0.61453126]
 [-1.58467315  0.62676957 -1.15883568]
 [ 1.23310723 -0.99984496  1.46174211]
 [-0.76107491  1.22346987 -0.91743769]]
'''

print demeaned.mean(0)
'''
[0,0,0]
'''

下图形象的展示了上面广播过程。

一维数组在轴0上的广播。

我们以二维数组(4,3)为例来说明轴方向的意义：

• 所谓轴0方向就是沿着第一行到第二行再到第三行的方向，在二维数组中可记为行方向。在上图中就是竖直方向。
• 所谓轴1方向就是沿着第一列到第二列再到第三列的方向，在二维数组中可记为列方向。在上图中就是水平方向。

• 我们可以这样记忆，沿着行方向是对每列进行操作，沿着列方向是对每行进行操作。

  这里需要注意理解轴方向的概率，axis=0（轴０方向）即为行方向，那么何为行方向呢？例如上面４行３列的二维数组是由４个长度为３的行向量堆叠
  而成，那么所谓行方向就是沿着第一行到第二行再到第三行的方向，在上图中就是竖直方向。axis=1（轴１方向）即为列方向，也就是第一列到第二列
  再到第三列的方向，上图中的水平方向。轴方向是一个非常重要的概念，在许多对数据操作中都会遇到，务必要深刻理解。
  我们可以这样记忆，沿着行方向是对每列进行操作，沿着列方向是对每行进行操作。

广播的原则

如果两个数组的后缘维度（即从末尾开始算起的维度）的轴长度相符或者其中一方的长度为1（注：这里是指两个数组同时从低维开始一一的比较，不能把其中一个数组的末尾第一个维度和第二个数组的末尾第二维度比较，只能是同时比较各自的末尾第一个），则可认为他们是广播兼容的。广播会在缺失或长度为１的维度上进行。

例如：上面的 二维(4,3)矩阵，得到的广播矩阵 (3,)是一个 一维向量，长度为３，我们从后缘维度起进行对比，二维矩阵最低维为３，广播矩阵最低维度也为３，对应的上，故广播兼容，可进行广播操作。我们也可将广播矩阵从一个长度为３一维向量reshape成一个二维矩阵(1,3)进行广播，此时两矩阵从最末尾一一比较，３和３相符，４和１相符。故可兼容广播，并得到与(3,)同样的广播结果。那么可能会有人问能不能将shape为(3,)的向量reshape成(3,1)，我们按照广播原则从低维开始一一比较,两数组的最末尾3和1相符，末尾倒数第二位4和3不相符，故reshape成(3,1)不能广播兼容。

例如，我们要在一个二维数组(4,3)，对其每一个元素减去他所在行的平均值。

import numpy as np
from numpy.random import randn

np.random.seed(1)
arr=randn(4,3)
#print arr
'''
[[ 1.62434536 -0.61175641 -0.52817175]
 [-1.07296862  0.86540763 -2.3015387 ]
 [ 1.74481176 -0.7612069   0.3190391 ]
 [-0.24937038  1.46210794 -2.06014071]]
'''

print "arr.mean(1):",arr.mean(1) ##沿着列方向计算每行的均值

'''
arr.mean(1): [ 0.1614724  -0.83636656  0.43421465 -0.28246772]
'''

print arr.mean(1).shape##(4,)一维向量，长度为４

demeaned=arr-arr.mean(1)
print demeaned

报错：ValueError: operands could not be broadcast together with shapes (4,3) (4,) 

显然根据上面的广播原则，两个矩阵的shape分别为(4,3)和(4,)，我们从低维开始一一比较，３和４明显不一致故不能兼容广播。广播原则中说到从低维比较，相符或者其中一个长度为１，故我们可以把shape为(4,)的数组reshape成(4,1)，那么两数组从末尾比较3和１，４和４均相符，故兼容广播。

demeaned=arr-arr.mean(1).reshape(4,1)##每个元素减去这个元素所在列的均值,取均值化

print demeaned
'''
[[ 1.46287296 -0.77322881 -0.68964415]
 [-0.23660206  1.70177419 -1.46517213]
 [ 1.31059711 -1.19542155 -0.11517556]
 [ 0.03309734  1.74457565 -1.77767299]]
[ 0.64249134  0.36942487 -1.01191621]
'''

print demeaned.mean(1)
'''
[0,0,0,0]
'''

下图形象的展示沿着列方向（轴１方向）的广播过程

同理三维数组上在轴０方向上的广播

显然上述各种数组的广播形式在shape上都符合广播原则。

沿其他轴向广播

对于三维情况，在三维中任何一维上的广播其实就是将数据重塑为兼容的的形状而已。

于是就有了一个非常普遍的问题，即专门为了广播而添加一个长度为１的新轴。虽然reshape是一个办法，但插入轴需要构造一个表示新形状的元组。

numpy数组提供了一种通过索引机制插入轴的特殊语法。下面这段代码通过特殊的np.newaxis属性以及“全”切片来插入新轴：

import numpy as np
from numpy.random import randn


arr=np.zeros((4,4))
'''
[[ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]
 [ 0.  0.  0.  0.]]
'''

arr_3d=arr[:,np.newaxis,:]##在1位置上插入一个长度为1的新轴
#arr_3d.shape:(4, 1, 4)

arr_1d=np.random.normal(size=3)
'''
[ 1.08181214  1.31445112  0.602738  ]  arr_1d.shape:(3,)即长度为３的一维向量
'''

arr_1d[:,np.newaxis] ##在1位置处添加一个长度为1的新轴，使其shape变为(3,1)
'''
[[ 1.77574557]
 [ 0.16187092]
 [-0.0083977 ]]
'''

再尝试沿着其他轴方向进行广播操作：

import numpy as np
from numpy.random import randn


arr=randn(3,4,5)
'''
[[[ 0.11184253 -0.15279816  0.23680397 -1.22944714 -0.13064376]
  [ 0.8787086  -0.62844149 -0.60911339  0.03564705  2.96975448]
  [ 1.15338316  0.38345286  2.27896049 -0.36889857 -0.18188945]
  [ 0.67825372 -0.84459366 -0.33390995  0.75827916 -0.91231555]]

 [[ 0.18424474 -0.34040063 -0.25759413 -1.23163447  0.21883981]
  [ 0.23571121 -1.58454497 -0.40426192  0.38485706 -1.66662353]
  [ 0.0591205   1.69397703 -0.62524794 -0.65350386  1.51347344]
  [-0.04532048 -0.37944063 -1.08981511 -0.61746888  1.03732066]]

 [[ 0.23612507 -0.84245519 -0.86550754  1.69379564 -0.93494228]
  [ 0.61076779 -1.24393264  2.3993014   1.47599124 -0.9150117 ]
  [-1.62919265 -1.53302405  0.217103    0.81751423 -0.01095928]
  [ 0.16081681  0.29020805  0.3932391  -0.10560003 -0.31416152]]]

  显然是由３个二维(4,5)矩阵堆叠起来的
'''

depth_mean=arr.mean(2)##在shape为(3,4,5)的三维矩阵中，沿着轴２方向计算均值，会得出shape为(3,4)的二维矩阵。
'''
[[ 0.20483366 -0.80578108 -0.51208484 -0.49015926]
 [-0.39936945  0.08758814 -0.31015946 -0.08575578]
 [-0.76100035 -0.25742263  0.2735308  -0.0159981 ]]
'''

demaned=arr-depth_mean[:,:,np.newaxis]##以depth_mean为广播矩阵进行广播，为了符合广播原则，需要插入新轴
'''
[[[-0.01317216  0.14490089 -0.66325081 -0.16924289  0.70076497]
  [-1.14228656 -0.14249658  0.24861719  0.00235415  1.0338118 ]
  [ 0.70445284  0.337685   -0.62309796  0.67583763 -1.0948775 ]
  [ 1.76712889 -0.93199695 -0.18065161 -1.29275054  0.6382702 ]]

 [[ 0.07375443  0.82294686 -0.08986837 -0.48565195 -0.32118097]
  [ 0.21026802 -1.17308564 -0.41494633 -0.76211049  2.13987444]
  [ 0.63042099 -0.55465847  0.23111804 -0.58707723  0.28019666]
  [ 0.33952837 -0.54703888 -0.82572159  0.4372118   0.5960203 ]]

 [[-0.36926057  1.79907294 -1.93419447  1.48362004 -0.97923794]
  [ 0.46237558 -0.24116939 -0.08864153  0.7037832  -0.83634786]
  [ 1.52145771 -1.65248743  0.08076382 -0.62356982  0.67383571]
  [ 0.65008267 -0.09869905  0.26594781 -1.04780342  0.23047198]]]
'''

np.newaxis：其实就是在指定位置插入长度为１的新轴。

在三维或三维以上进行广播时，为了广播兼容，会根据广播原则，对广播矩阵进reshape和np.newaxis（插入新轴）操作，其实都是大同小异。只要在广播时，严格根据广播原则即可。

利用矩阵运算实现逻辑回归

逻辑回归的原理可以参考我的这篇博文： 机器学习->监督学习->logistic回归,softMax回归

下面讲怎么利用矩阵运算快速实现逻辑回归。

数据归一化（去均值）

def normalization(X):
    '''
    X标准化：以一列作为一个单位，标准化公式(x-该列均值）/（该列标准差）
    因为给的X有两列，Y一列，则对应的theta有三列，其中theta0对应的X默认是1，故从X的第二列开始标准化。
    '''
    for i in range(1,int(X.shape[1])):
        std,mean=np.std(X[:,i],ddof=1),np.mean(X[:i])##标准差公式里面默认(ddof=0)除以N，就是数组长度,ddof=1就是除以(N-1)
        X[:,i]=(X[:,i]-mean)/(std+1e-10)
    return X

上述代码利用了向量运算，每一次循环都是以向量为单位做整体的去均值化。因为除去第一列不需要作归一化，故不能以矩阵作为整体来归一化。

计算模型预测值

def sigmoid(X,T):
    '''
    :param X: 样本的特征矩阵
    :param T: 参数向量
    :return: 返回预测概率
    '''
    return 1.0/(1+np.exp(-(X.dot(T))))

上述代码中Ｘ矩阵dot乘Ｔ向量得到一个新的向量，np.exp对新的向量作元素级　操作，其结果返回一个新的向量。返回向量里的每个元素对应每个样本的预测结果。

计算预测值与真实值的均方误差

def MSError(y,y_hat):
    '''
    :param y:真实y向量
    :param y_hat: 预测向量
    :return: 返回均方误差
    '''
    #return ((y-y_hat).transpose().dot((y-y_hat)))/(2*y.shape[0])或
    return sum(np.power(y-y_hat,2))/(2*y.shape[0])

y是一个向量，y_hat也是一个向量，均方误差：

显然我们用y向量减去y_hat向量，在这两个向量之间做的是元素级操作，得到一个(y-y_hat)向量，np.power对 (y-y_hat) 向量作元素级操作，也就是对各元素作平方运算，得到一个新的向量。然后通过sum对向量作聚合运算，就是计算这个向量的总和。（这一步也可以用y-y_hat矩阵转置乘以矩阵本身）

梯度下降

def GD(X,Y):
    epsilon=1e-8
    alph=0.001
    error1=0
    cnt=0
    m=np.shape(Y)[0]
    T=np.ones((X.shape[1],1))
    while True:
        cnt=cnt+1
        y_hat=sigmoid(X,T)
        error=y_hat-Y
        T=T-alph*X.transpose().dot(error)## 根据error×feature形式更新参数矩阵
        error0=error1
        error1=MSError(Y,y_hat)##计算均方误差
        if abs(error1-error0)<epsilon:##如果上一次的均方误差和当前次的均方误差很小（可视为误差收敛），则停止
            break
        plot2D(X,Y,T,cnt,error1)
    print "iterations",cnt
    return T