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

# coding: utf-8


# # Numpy and Linear Algebra
# qcy
# 2017年8月12日16:11:05

import math
import shutil
import pandas as pd
import numpy as np
import os
from urllib.request import urlopen, quote
import ffn
import matplotlib as mpl
import matplotlib.pyplot as plt


mpl.rcParams['axes.unicode_minus']=False #用来正常显示负号
from matplotlib.ticker import FuncFormatter
from matplotlib import rcParams
from matplotlib.font_manager import FontProperties


import scipy.optimize as sco

rcParams.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.T
print(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*y
print(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 = A
print(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 Rank
print('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)


# 行列式 Det
print('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 norm
print(np.sqrt((A**2).sum())) # calc Fro. norm of a matrix by def


np.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)


# 另一种方法，直接调用solve
x2 = np.linalg.solve(A,b)
print(x2)



# Check
print(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.T
A2 = A + A.T
print(A2)


# 若特征值全部大于0，则正定
res = np.all(np.linalg.eigvals(A2)>0) # 全部大于0
print(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))