In [1]:
import itertools as itt
import numpy as np
import numpy.random as nprnd
import scipy.stats as spstat
import scipy.special as ssp
import pandas as pd
import bokeh
import bokeh.plotting as bplt
import bokeh.palettes as bpal
In [2]:
# 工具配置
TOOLS = "pan,wheel_zoom,box_zoom,box_select"
# 作图配置
plot_config = dict(plot_width=800, plot_height=500, tools=TOOLS)
bplt.output_notebook()
我们先了解一下矩阵最重要也是最基本的工程用途：通过矩阵求解多元一次方程组，这已经成功应用于农业（？）和金融业。
本节对于矩阵的相关概念、计算细节都不做任何解释或者证明。
鸡兔同笼问题是经典的小学竞赛题：
问鸡和兔子分别有多少？
把$y$写成关于$x$的表达式并代入的方法的本意，是用来得到一个消去$y$的式子从而求得$x$，再根据$x,y$的关系求得$y$
第一步： $$ x+y = 20\\ 4x+2y = 62 $$
第二步，用第二个式子减去2倍的第一个式子，然后扔掉第一个式子： $$ 2x = 22\\ 4x+2y = 62 $$
后面的故事还是很熟悉，但是我们为什么可以扔掉第一个式子呢？因为 新来的第一个式子和旧有的第二个式子可以恢复出旧有的第一个式子，因此我们能够利用的信息并没有变多，也没有变少。
第三步，第一个式子缩水成为原来的1/2： $$ x = 11\\ 4x+2y = 62 $$
最后一步，第二个式子先减去第一个式子，再缩水到1/2：
根据上式推出
$$\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 11\\ 9 \end{bmatrix}$$左边方框包住的方方正正的东西就是一个2*2的系数矩阵了。中间未知数代表的矩阵只有1列所以也可以叫列向量，右边的1列也是列向量。简化一下：
$$Ax = b$$为了求得未知数列向量$x$，我们可以直接使用：
$$x = A^{-1}b$$来得到，这就简单多了！
In [3]:
A = np.mat([[1,1],[4,2]])
A_Inv = A.I
b = np.array([20,62])
A_Inv.dot(b)
Out[3]:
In [4]:
A = np.mat([[1,3],[4,12]])
A_Inv = A.I
这种情况下，矩阵的某些行通过加减乘除的组合能够完全得到另一些行，称这个矩阵是线性相关的，也可以称作奇异（singular）
最终的结果是：去掉那些多余的行之后，
假设我们要建立一个关于某只股票X的组合，要求这个组合对于X股价的一阶以及二阶变化都不敏感。组合包含：
a份股票 $\Delta = 1.0, \Gamma = 0.0$
b分某种看涨期权 $\Delta = 0.1, \Gamma = 0.05$
c份某种看空期权 $\Delta = -0.2, \Gamma = 0.15$
$$\begin{bmatrix} 1.0 & 0.1 & -0.2 \\ 0.0 & 0.05 & 0.15 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0 \end{bmatrix}$$
In [5]:
# 先补齐最后一行
A = np.array([[1,0.1,-0.2],[0,0.05,0.15],[0,0,0]])
d,U = np.linalg.eig(A)
D = np.diag(d)
print D
print U
这时我们有：$AU = UD$
其中，$U$的每一个列向量$v$在左边乘以了$A$之后都变成了$v$自己的倍数。
$D$是$A$这个矩阵特征值组成的对焦矩阵，$D$里面的每一项，都顺序地对应着$U$的每一个列向量。我们称$U$的列向量是
我们要把$Av = 0*v = 0$的向量找到，也就是：
In [6]:
U[:,2]
Out[6]:
至此我们知道，15.6份股票、做空93.7份看涨期权与买入31.2份看跌期权可以实现我们构建市场中性组合的目的。
简称 $Ax = b$
可以看到，模型很粗糙，甚至还有矛盾，比如(3),(4)式的平均无法得出(1)式。
不过不要紧，我们给出一个$x = (A^TA)^{-1}A^T b$的解
In [7]:
A = np.array([[1.0,1.0],[1.0,0.7],[1.0,0.9],[1.0,1.1]])
b = np.array([3500,2700,3600,3350])
x = (np.mat(A.T.dot(A)).I).dot(A.T).dot(b.reshape(4,1))
b_ = A.dot(x)
print x
print b_
这使得：
其建模方法论可以简单的利用情况3的代码亲自实践。
In [8]:
A = np.array([[5,0,-2],[3,1,1]])
print A[0,1],A[1,0]
我们平日里习惯称矩阵$A$的第1行第2列元素为0，第2行第1列元素为3，就这样被取出来了。上图我们还可以看出使用了双重方括号来表示矩阵。
分别看一下矩阵A按照行列范围切出来的部分与直接取出一行/一列有什么区别：
In [9]:
print A[0:1,:]
print A[0,:]
In [10]:
print A[:,0:1]
print A[:,0]
In [11]:
B = np.array([[0,3,3],[-2,7,-1]])
print A
print B
print A+B
print A+A
print 2*A
print 1.5*A
如上，读者可以自己验证矩阵加法满足交换律与结合律：
$A+B=B+A,(A+B)+C=A+(B+C)$
矩阵的乘法虽然看起来不那么易懂，但就是从方程组抽象而来：
中俄边境有一些农贸市场，也有很多风格的餐馆。
农贸市场中商品的标价是以两种货币的计价（汇率有微妙的不同，因为两边商品的保存与生产成本不一样）
两家中俄连锁饭店菜品不一样，因此每个月进货也是不一样的：
假设餐馆A要以人民币计价进行采购，我们分别列出价格矩阵和数量矩阵：
$$P = \begin{bmatrix} 10 & 3 & 8 \end{bmatrix} , Q = \begin{bmatrix} 500 \\ 200 \\ 200 \end{bmatrix}$$
In [12]:
P = np.array([[10,3,8]])
Q = np.array([[500],[200],[200]])
print P.dot(Q)
print 10*500+3*200+8*200
总结一下矩阵乘法的规律：当$R = PQ$时，有
$$R_{ij} = \sum_{k} P_{ik}Q_{kj}$$矩阵之所以很方便，就是因为，我们可以一次性的计算出2家餐馆分别以2种货币进货时所需要付出的代价：
我们刚才知道以人民币计价餐馆A每个月进货需要7200元。
In [13]:
P = np.array([[60,30,120],[10,3,8]])
Q = np.array([[500,100],[200,100],[200,700]])
print P.dot(Q)
如果我们想知道以卢布计价餐馆B每个月的进货成本，直接可以从矩阵里看出是93000卢布。
In [58]:
A = np.array([[1,2],[6,-2]])
B = np.array([[2,3],[-5,1]])
In [59]:
print A.dot(B)
In [60]:
print B.dot(A)
In [61]:
print A
print A.T
同时根据矩阵加法和乘法的定义，能够知道：$$(A+B)^T = A^T+B^T\\ (AB)^T = B^TA^T$$
In [62]:
print (A+B).T
print A.T+B.T
In [63]:
print (A.dot(B)).T
print B.T.dot(A.T)
还是从鸡兔同笼的问题开始看：假设有x只兔子，y只鸡：
$$\begin{bmatrix} 1 & 1 \\ 4 & 2 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 22\\ 62 \end{bmatrix}$$如果我们把如上的矩阵写成$Ax = b$，可以看出，矩阵$A$中，第一列代表兔子的头/腿数，第二列代表鸡的头/腿数。改变写法
$$A[:,0]x+A[:,1]y = b$$可以看到，$b$是$A$中列向量乘以标量之和。此时，
In [17]:
A = np.array([[1,1,22],[4,2,62]])
np.linalg.matrix_rank(A)
Out[17]:
In [18]:
A = np.array([[4,1],[1,2]])
e1 = np.array([[1],[0]])
e2 = np.array([[0],[1]])
In [19]:
c1 = A.dot(e1)
print c1
In [20]:
c2 = A.dot(e2)
print c2
In [21]:
p1 = bplt.figure(title='向量与面积的映射',background_fill='#FFFFFF',**plot_config)
p1.line(*zip(np.array([0,0]),e1.flatten()),line_width=3)
p1.line(*zip(np.array([0,0]),e2.flatten()),line_width=3)
p1.patch(*zip(np.array([0,0]),e1.flatten(),(e1+e2).flatten(),e2.flatten()),alpha=0.3)
p1.line(*zip(np.array([0,0]),c1.flatten()),line_width=3,line_color='darkred')
p1.line(*zip(np.array([0,0]),c2.flatten()),line_width=3,line_color='darkred')
p1.patch(*zip(np.array([0,0]),c1.flatten(),(c1+c2).flatten(),c2.flatten()),alpha=0.3,color='red')
bplt.show(p1)
经过简单的切割法，我们知道向量$e1,e2$经过矩阵$A$的映射之后，形成的红色平行四边形的面积是蓝色正方形的7倍。
这就涉及到一个概念：行列式。
In [22]:
A1 = np.array([[4,1],[1,2]])
print np.linalg.det(A1)
In [23]:
A2 = np.array([[1,4],[2,1]])
print np.linalg.det(A2)
行列式的意思就是：一个矩阵把一个平行四边形映射之后放大的面积。
我们注意到如果交换了$A1$的两列变成了$A2$，则行列式变成了负数。
原本$e1,e2$是逆时针排列的，映射后$c1,c2$变为顺时针排列，因此行列式要带负号。
而且从这种映射关系我们了解到，尽管矩阵乘法无法对于任意$A,B$满足$AB=BA$，但是可以满足
$$det(AB) = det(A)det(B)=det(B)det(A)=det(BA)$$尽管刚才我们都在计算2阶矩阵(2*2矩阵)的行列式，这个规则同样也适用于3阶甚至高阶的矩阵
In [49]:
D = np.array([[5,0],[0,3]])
e1 = np.array([[1],[0]])
e2 = np.array([[0],[1]])
c1 = D.dot(e1)
c2 = D.dot(e2)
p1 = bplt.figure(title='对角矩阵对向量与面积的映射',background_fill='#FFFFFF',**plot_config)
p1.line(*zip(np.array([0,0]),e1.flatten()),line_width=3)
p1.line(*zip(np.array([0,0]),e2.flatten()),line_width=3)
p1.patch(*zip(np.array([0,0]),e1.flatten(),(e1+e2).flatten(),e2.flatten()),alpha=0.3)
p1.line(*zip(np.array([0,0]),c1.flatten()),line_width=3,line_color='darkred')
p1.line(*zip(np.array([0,0]),c2.flatten()),line_width=3,line_color='darkred')
p1.patch(*zip(np.array([0,0]),c1.flatten(),(c1+c2).flatten(),c2.flatten()),alpha=0.3,color='red')
bplt.show(p1)
In [50]:
np.linalg.det(D)
Out[50]:
In [24]:
T1 = np.mat(A1)
S1 = T1.I
print T1.dot(S1)
print S1.dot(T1)
print np.eye(2)
不难知道，单位矩阵$I$能把空间中任何的平行多面体映射到自身，因此放大倍数永远是1，$det(I)=1$永远成立。
In [25]:
A = np.array([[4,3],[1,2]])
v1 = np.array([[3],[1]])
v2 = np.array([[1],[-1]])
In [26]:
print A.dot(v1)
print A.dot(v2)
我们把刚才的关系理顺一下，用更简单的方式表达：
如果我们把特征值写到一个矩阵的对角线里，再把特征值排成一排写到另一个矩阵里：
$$D = \begin{bmatrix} 5 & 0 \\ 0 & 1 \end{bmatrix},V = \begin{bmatrix} 3 & 1 \\ 1 & -1 \end{bmatrix}$$则有关系：
$$AV = VD$$如果能找到一个可逆的$V$，则
$$V^{-1}AV = D$$
In [27]:
A = np.array([[4,3],[1,2]])
V = np.array([[3,1],[1,-1]])
D = np.array([[5,0],[0,1]])
print A.dot(V)
print V.dot(D)
这样说来矩阵$V$,$D$的性质相当不错。如何找到$V$和$D$呢？
可以看到，特征分解函数eig可以轻松帮我们找到：
In [77]:
d,V = np.linalg.eig(A)
print d
print np.diag(d)
print V
这时我们发现，算法给出的$V$矩阵每一列都是我们自己找到的$V$矩阵中对应列的倍数，所以这个新来的$V$的每一列仍然是特征向量：
$$Av = dv \rightarrow A(kv) = k(Av) = k(dv) = d(kv)$$问题的提出：某个游戏推出一个升级系统
问，当整个系统玩家上升到一定数量的时候，各个等级之间的比例将是怎样的？
我们要用更好的方式来描述问题：
而所谓根据一定的几率进行跳转，这个几率我们仍然可以用一个状态转移概率矩阵来描述：$P_{ij}$表示从状态$i$跳转到状态$j$的概率：
$$P = \begin{bmatrix} 0.85 & 0.15 & 0 & 0 & 0\\ 0.40 & 0.45 & 0.15 & 0 & 0 \\ 0 & 0.40 & 0.45 & 0.15 & 0 \\ 0 & 0 & 0.40 & 0.45 & 0.15 \\ 0 & 0 & 0 & 0.40 & 0.6 \end{bmatrix}$$容易发现这个5*5的概率矩阵每一行之和为1。不加证明的给出：
对于一个用行向量表示的概率状态$v^T$, $v^TP$将表示它在下一步跳转发生后的概率状态。
In [53]:
P = np.array([[0.85,0.15,0,0,0],[0.4,0.45,0.15,0,0],[0,0.4,0.45,0.15,0],[0,0,0.4,0.45,0.15],[0,0,0,0.4,0.6]])
v = np.array([0,0,1,0,0])
print v.dot(P)
由上式我们可以看出，当$v^T$表示一个确定的2级状态的时候，$v^TP$能够表示它在下一步的状态概率分布：
我们最初的问题变成了$v^TP = v^T$，求$v^T$，因为对于上述的式子转置后得到
$$P^T v = v$$因此$v$就是$P^T$这个矩阵对应特征值为1时候的特征向量
In [55]:
d,V = np.linalg.eig(P.T)
print np.diag(d)
print V
我们取出$V[:,2]$，并且对它进行归一化：
(归一化的意思是对它进行缩放，使得这个向量中的元素之和为1，以满足概率分布的要求)
In [57]:
v = V[:,2]/np.sum(V[:,2])
print v
print v.dot(P)
由此我们知道，如果一个系统中武器等级的概率分布为$v$，状态概率转移矩阵为$P$，那么下一步跳转后，状态分布仍然还是$v$。
这个$v$被称作平稳分布，我们也可以很轻松的知道：
在用户足够多的时候，武器满级的玩家大概占所有人的1.245%。
这样，游戏的数值策划们就可以通过微调矩阵$P$来让拥有各个武器等级比例更加合理了。
在我们了解了对角化之后，矩阵的乘方会变的特别容易，因为
$$A = VDV^{-1} \\ A^{N} = AAA...A = VDV^{-1}VDV^{-1}...VDV{-1} = VD^{N}V^{-1}$$我们把曾经介绍过的Fibonacci数列问题再一次提出来，并寄希望于矩阵法能够更快的求解出结果：
$$\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\begin{bmatrix} F_{n+1}\\ F_{n} \end{bmatrix} = \begin{bmatrix} F_{n+2}\\ F_{n+1} \end{bmatrix}$$反复应用上式，容易使用一个简洁的公式表达递推关系： $$\begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^N \begin{bmatrix} F_{1}\\ F_{0} \end{bmatrix} = \begin{bmatrix} F_{N+1}\\ F_{N} \end{bmatrix}$$
In [2]:
FibMat = np.array([[1,1],[1,0]])
d,V = np.linalg.eig(FibMat)
print np.diag(d)
print V
又看到了我们熟悉的1.618和0.618（黄金分割点）。
如果矩阵能够对角化的话，计算起来几乎仅仅耗费对角化的时间，比如我们想求$F_{17}$：
In [3]:
# N=16，先求FibMat的16次方
FibMat16 = V.dot(np.diag(d**16)).dot(np.mat(V).I)
print FibMat16
In [4]:
# 再乘以向量(1,1)
print FibMat16.dot(np.array([[1],[1]]))
上面的数字2584就是$F_{17}$了，可以看到求值的速度那是相当之快！
退一万步讲，就算一个矩阵不能成功的对角化，通过矩阵反复的乘方也比Fibonacci数列强行循环来的快。工具好就是不一样。
In [44]:
# 循环跟自己乘方四次：
Mat = FibMat
for i in xrange(4):
Mat = Mat.dot(Mat)
print Mat
有一类特殊的矩阵，满足条件：$T^{-1} = T^T$，被称为正交矩阵。也就是这个矩阵的逆就是它自己的转置。
例如，以下形式的二阶矩阵就都满足正交矩阵的性质：
$$T = \begin{bmatrix} \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{bmatrix}$$根据矩阵乘法的定义，容易知道这类矩阵有几个很好的性质：
In [70]:
theta = np.pi/6
T = np.array([[np.cos(theta),np.sin(theta)],[-np.sin(theta),np.cos(theta)]])
print T
In [71]:
e1 = np.array([[1],[0]])
e2 = np.array([[0],[1]])
c1 = T.dot(e1)
c2 = T.dot(e2)
p1 = bplt.figure(title='正交矩阵对向量与面积的映射',background_fill='#FFFFFF',**plot_config)
p1.line(*zip(np.array([0,0]),e1.flatten()),line_width=3)
p1.line(*zip(np.array([0,0]),e2.flatten()),line_width=3)
p1.patch(*zip(np.array([0,0]),e1.flatten(),(e1+e2).flatten(),e2.flatten()),alpha=0.3)
p1.line(*zip(np.array([0,0]),c1.flatten()),line_width=3,line_color='darkred')
p1.line(*zip(np.array([0,0]),c2.flatten()),line_width=3,line_color='darkred')
p1.patch(*zip(np.array([0,0]),c1.flatten(),(c1+c2).flatten(),c2.flatten()),alpha=0.3,color='red')
bplt.show(p1)
In [72]:
np.linalg.det(T)
Out[72]:
In [64]:
S = np.array([[3,1],[1,3]])
print S
print S.T
要说明对阵矩阵有什么好处，就要先说明其他矩阵有什么不好：
In [66]:
M1 = np.array([[1,1],[0,1]])
d,V = np.linalg.eig(M1)
D = np.diag(d)
print D
print V
In [67]:
print M1.dot(V)
print V.dot(D)
尽管$M_1$是一个维数为2的矩阵，而且$M_1V=VD$也成立，可是$V$并不是可逆矩阵，不能满足$V^{-1}M_1V = D$，这种情况我们称$M_1$无法对角化。
In [69]:
M2 = np.array([[4,8],[-5,2]])
d,V = np.linalg.eig(M2)
D = np.diag(d)
print D
print V
In [75]:
print np.float64(S)
d,V = np.linalg.eig(S)
D = np.diag(d)
print D
print V
在电影评论平台上，后台记录了形形色色的用户对电影评价的记录。
可以用有一个大矩阵X来粗略的保存"用户-电影"之间的关系：
如果我们希望发现这个大矩阵中存在某种模式，使得一个低维数的矩阵能够在数值上尽量接近原来的矩阵，那么我们需要SVD算法。
In [84]:
X = np.random.randint(-1,2,(30,10))
print X
我们随机生成了30位用户对10部电影的偏好。SVD算法的用处在于，能够找到两个正交矩阵$U,V$，使得： $$X = UDV\\ U^TXV^T = D $$
这里很明显，$X$并不是方阵，因此$U$是一个30阶方阵，$V$是一个10阶方阵，而$D$是一个非方形的对角矩阵
In [97]:
U,d,V = np.linalg.svd(X)
D = U.T.dot(X).dot(V.T)
In [100]:
print np.round(D,2)
当我们抛弃掉$D$矩阵对角元素末位的一些"信息"之后，很可能发现原来的矩阵并没有太大的损失.
而且我们可以假定，越小的对角元素所对应的信息越可能是"噪声"，也就是原来矩阵$X$中体现用户对电影随机偏好的那部分。
In [101]:
np.sum((d**2)[:6])/np.sum(d**2)
Out[101]:
我们可以获得一个修正后的modD矩阵，这个矩阵仅有6维：
In [103]:
modD = D
for i in xrange(6,10):
modD[i,i] = 0.0
print np.round(modD,2)
使用modD矩阵重构出来的modX矩阵可以被视作去掉了噪声的X矩阵：
In [108]:
modX = U.dot(modD).dot(V)
print np.round(modX,2)
这时我们分别从X矩阵和modX矩阵中拿出第0位用户的电影评价数据来做比较：
In [111]:
print np.round(np.float64(X[0,:]),2)
print np.round(modX[0,:],2)
用户0对于0号、2号、4号电影的评价都是0，假设这3部电影他都没看过。
而重构出来的矩阵中，他对这三部电影的评价的评分则是0.4,-0.42,0.8。这时你就可以试探性的优先推荐4号电影，而不要急着推荐2号电影了。
In [ ]: