Python：用Numpy解决线性代数的问题

# coding: utf-8# # Numpy and Linear Algebra# qcy# 2017年8月12日16:11:05import mathimport shutilimport pandas as pdimport numpy as npimport osfrom urllib.request import urlopen, quoteimport ffnimport matplotlib as mplimport matplotlib.pyplot as pltmpl.rcParams['axes.unicode_minus']=False #用来正常显示负号from matplotlib.ticker import FuncFormatterfrom matplotlib import rcParamsfrom matplotlib.font_manager import FontPropertiesimport scipy.optimize as scorcParams.update({   # 'font.family': 'sans-serif',    'font.sans-serif': ['SimHei'],    'font.size': 14})# # 矩阵A = np.mat([[1,2,3],[4,5,6]])print(A, type(A), A.shape)At = A.Tprint(At)print(At.shape)# # 二维数组# 其实和矩阵真的挺像的，我暂时还不知道有什么区别B = np.array([[1,2,3,4],[5,6,7,8]])print(type(B), B.shape)print(B)print(B+2)print(B*2)print(B/2)print(B-2)print(B.T)# # 一维数组# 就是向量？x = np.array([1,0,3])y = np.array([1,2,3])print(x,y)# 点对点地乘z1 = x*yprint(z1)# 内积z2 = x.dot(y)print(z2)z3 = y.dot(x)print(z3)# 常见特殊矩阵diagY = np.diag(y)print(diagY)I3 = np.eye(3)print(I3)# # 矩阵运算A = np.array([[1,2,3,4],[5,6,7,8]])print(A)B = Aprint(id(A),id(B))B = A.copy() # deep copy, not just ref.print(id(A),id(B))C1 = A*B # 点乘print(C1)C2 = A.dot(B.T) # 矩阵乘法print(C2)# Trace and Rankprint('Trace of C2 is', np.trace(C2))print('Rank of C2 is', C2.ndim) # np.rank() 已过时……# ?????#  这里会涉及到rank的概念，在线性代数(math)rank表示秩，但是必须明确的是在numpy里#  rank不是表示秩的概念，是表示维数的概念，这个理解的话需要看此文章：对于多维arrays的数据结构解释：# [多维arrays数据结构理解][1]print('np.linalg.matrix_rank(C2):',np.linalg.matrix_rank(C2))# # 逆矩阵、伴随矩阵np.random.seed(0)A = np.random.rand(3,3)print(A)# 行列式 Detprint('Det of C2 is', np.linalg.det(A))# 逆矩阵A_1 = np.linalg.inv(A)print(A_1)print(A.dot(A_1))# 伴随矩阵A_star = A_1 * np.linalg.det(A)print(A_star)# # 关于范数A = np.array([[1,2,3,4],[5,6,7,8]])print(A)print(np.linalg.norm(A))  # Default: Frobenius normprint(np.sqrt((A**2).sum())) # calc Fro. norm of a matrix by defnp.random.seed(87)x = np.random.rand(3)print(x)print(np.sqrt(x.dot(x)))# x = np.random.rand(10,1)# 向量的范数# 自己 dot 自己，然后开平方# # 线性方程组A1 = [[1,2,1],[2,-1,3],[3,1,2]]A = np.array(A1)print(A)# b是一行的话，也能解？？b = np.array([7,8,18])b = b.reshape(-1,1) # reshpae可以把np.array([7,8,18])转成一列，但transpose没有用print(type(b),b.shape)b = np.array([[7],[8],[18]]) # 列也可以print(type(b),b.shape)print(b)x = np.linalg.inv(A).dot(b)print(x)# 另一种方法，直接调用solvex2 = np.linalg.solve(A,b)print(x2)# Checkprint(A.dot(x))print(A.dot(x2))# # 矩阵的特征值A = np.diag([1,2,3])print(A)[eig_vals,eig_vecs] = np.linalg.eig(A)print('特征值：',eig_vals)print('特征向量：',eig_vecs[0],eig_vecs[1],eig_vecs[2])a0 = A.dot(eig_vecs[0]) # 矩阵乘以向量print(a0 - eig_vals[0]*eig_vecs[0])a0 = np.dot(A,eig_vecs[0]) # 矩阵乘以向量 也可以print(a0)# (numpy.dot(x,b[:][i])==a[i]*b[:][i]).all():# # 矩阵的奇异值分解# 注意，这里的svd出来的V不是V，而本来就是V.hermitian了吧？np.random.seed(10)A = np.random.rand(100,100)U,S,VH = np.linalg.svd(A) # print(U.dot(np.diag(S)).dot(VH))# print(A)plt.plot(S)plt.title('SVD of Matrix A')plt.grid()plt.show()# # 判定矩阵正定性np.random.seed(10)A = np.random.rand(3,3)print(A)# 构造一个对称矩阵 A+A.TA2 = A + A.Tprint(A2)# 若特征值全部大于0，则正定res = np.all(np.linalg.eigvals(A2)>0) # 全部大于0print(res)# Test 'any' & 'or'print(np.all(np.array([1,2,3,4,5])>0))print(np.any(np.array([-1,2,3,4,5])<0))# # 协方差矩阵、相关系数x1 = np.random.rand(100,2)# print(x1)C = np.cov(x1,rowvar=False)print(C)Rho = np.corrcoef(x1, rowvar=False)print(Rho)# # 随机变量线性组合的均值和方差np.random.seed(99)x = np.random.rand(100,2)w = np.array([1,1]).reshape(2,1)u1 = np.mean(x,axis=0)# x2 是 w[0] * x的第一列 + w[1] * x的第2列x2 = x.dot(w) # 随机变量的线性组合u2 = np.mean(x2)print('==均值==')print(sum(u1),u2) # 均值的线性性# 为什么会不一样，这是样本方差还是方差??print('==方差==')var2 = np.var(x2) # 变量线性组合的方差print(var2,type(var2))print('我算的方差',np.mean([(x-u2)**2 for x in x2]))C = np.cov(x,rowvar=False)var1 = (w.T).dot(C.dot(w)) # 这种算法居然算的是样本方差!!!二者差距很大!!!print(var1,type(var1))print('我算的样本方差',np.sum([(x-u2)**2 for x in x2])/(len(x2)-1))# 验证var是方差还是协方差x0 = [1,2,3,4,5]# np.mean(x0)var = np.mean([(x-3)**2 for x in x0]) # 这是方差S2 = np.sum([(x-3)**2 for x in x0])/(len(x0)-1) # 这是样本方差print(var, S2)print(np.var(x0))