首页 > 代码库 > 『科学计算』线性代数部分作业
『科学计算』线性代数部分作业
最小二乘法求解垂足
from matplotlib import pyplot as plt import numpy as np from mpl_toolkits.mplot3d import Axes3D A=np.array([[1],[2],[3]]) B=np.array([[1],[1],[1]]) x=np.linspace(-0.5,1,10) x.shape=(1,10) xx=A.dot(x) C=A.T.dot(B) AA=np.linalg.inv(A.T.dot(A)) P=A.dot(AA).dot(C) E=B-P fig = plt.figure() ax = fig.add_subplot(111,projection=‘3d‘) ax.plot(xx[0,:],xx[1,:],xx[2,:],label="lineA") ax.plot(A[0],A[1],A[2],‘ko‘,label="A") ax.plot([0,B[0]],[0,B[1]],[0,B[2]],‘m-o‘,label="0B") ax.plot([B[0][0],P[0][0]],[B[1][0],P[1][0]],[B[2][0],P[2][0]],‘r-o‘,label="BP") ax.plot([0,E[0]],[0,E[1]],[0,E[2]],‘y-o‘,label="0E") ax.legend() ax.axis(‘equal‘) plt.show()
因为是三维图么,所以扭了一下多截了一张(笑):
最小二乘法拟合函数
import numpy as np import matplotlib.pyplot as plt import copy # 产生一个方波(x,y) x = np.linspace(-10,10,300) y=[] for i in np.cos(x): if i>0: y.append(0) else: y.append(2) y=np.array(y) def projection(A,b): #### # return A*inv(AT*A)*AT*b #### AA = A.T.dot(A) w=np.linalg.inv(AA).dot(A.T).dot(b) print(w) return A.dot(w) def fourier(x,y,N): A = [] for i in x: A.append(copy.deepcopy([])) for j in range(N): A[-1].append(np.sin((j + 1) * i)) A[-1].append(np.cos((j + 1) * i)) A[-1].append(1) b = y.T yw = projection(np.array(A),b) return yw # write Your code, Fourier function plt.plot(x,y,color=‘g‘,label=‘origin‘) plt.plot(x,fourier(x,y,3),color=‘r‘,label=‘3‘) plt.plot(x,fourier(x,y,8),color=‘b‘,label=‘8‘) plt.plot(x,fourier(x,y,23),color=‘k‘,label=‘23‘) plt.legend() plt.axis(‘equal‘) plt.show()
从输出图可以直观看出来项数越多拟合效果越好:
使用动画模拟矩阵变换和特征值特征向量的关系
import numpy as np import copy import matplotlib.pyplot as plt from matplotlib import animation A=np.array([[3,1],[2,4]])/4.0 # write Your code fig = plt.figure() ax = fig.add_subplot(111) line1, = ax.plot([],[],c=‘r‘) line2, = ax.plot([],[],c=‘b‘) ax.axis(‘equal‘) ax.set_xlim(-2,2) ax.set_ylim(-2,2) ax.set_xticks(np.linspace(-2,2,5)) ax.set_yticks(np.linspace(-2,2,5)) ax.grid(True) def data(): point_num = 100 x = np.cos(np.linspace(0,2 * np.pi,point_num)) y = np.sin(np.linspace(0,2 * np.pi,point_num)) rot_x,rot_y = copy.deepcopy([]),copy.deepcopy([]) int_x,int_y = copy.deepcopy([]),copy.deepcopy([]) for i in range(point_num): rot_x.append(A.dot([x[i],y[i]])[0]) rot_y.append(A.dot([x[i],y[i]])[1]) int_x.append(x[i]) int_y.append(y[i]) yield (rot_x,rot_y,int_x,int_y) def update(data): line1.set_xdata(data[0]) line1.set_ydata(data[1]) line2.set_xdata(data[2]) line2.set_ydata(data[3]) return line1,line2 def init(): line1.set_data([],[]) line2.set_data([],[]) return line1,line2 ani = animation.FuncAnimation(fig,update,data,init_func=init,interval=100,repeat=False) plt.show()
这是个动画,所以我截了两张图来表示,可以看到同一变换矩阵对不同方向的单位向量放缩长度并不相同,是个椭圆形:
使用矩阵求解差分方程(发散情况)
import numpy as np import matplotlib.pyplot as plt import time def time_cost(f): def _f(*arg, **kwarg): start = time.clock() a=f(*arg,**kwarg) end = time.clock() print(f.__name__,"run cost time is ",end-start) return a return _f @time_cost def fib_opt_seq(seq): return [fib_opt(i) for i in seq] def fib_opt(n): a,b,i=0,1,0 while i<n: a,b=b,a+b i+=1 else: #print(b) return b import random #seq = [random.randint(800,1000) for i in xrange(1000)] seq = range(1000) a=fib_opt_seq(seq) # write Your code fib_eig_seq function A = np.array([[1,1],[1,0]]) eigval,eigvect = np.linalg.eig(A) S_inv = np.linalg.inv(eigvect) @time_cost def fib_eig_seq(seq): return [fib_eig(i) for i in seq] def fib_eig(n): af = (np.diag(eigval)**(n+1)).dot(S_inv) #print((eigvect.dot(af)).dot(np.array([[1],[0]]))[1]) return (eigvect.dot(af)).dot(np.array([[1],[0]]))[1] b=fib_eig_seq(seq)
Python 3.5.2 |Anaconda 4.2.0 (64-bit)| (default, Jul 5 2016, 11:41:13) [MSC v.1900 64 bit (AMD64)]
Type "copyright", "credits" or "license" for more information.
IPython 5.1.0 -- An enhanced Interactive Python.
? -> Introduction and overview of IPython‘s features.
%quickref -> Quick reference.
help -> Python‘s own help system.
object? -> Details about ‘object‘, use ‘object??‘ for extra details.
PyDev console: using IPython 5.1.0fib_opt_seq run cost time is 0.08637829184750435
fib_eig_seq run cost time is 0.01634134003747284
上面对比了传统的递归方式生成斐波那契数列方式,一般来说随着数量的上升,使用矩阵速度更快。
使用矩阵求解差分方程(收敛情况)
import numpy as np A = np.array([[ 0.8 , 0.1], [ 0.2 , 0.9]]) N_year = 1000 x=np.array([3200, 4000]) def hw_3_5(a,x,n): eigval, eigvact = np.linalg.eig(a) res = (eigvact.dot((np.diag(eigval)**n).dot(np.linalg.inv(eigvact)).dot(x.T))) print(res) return res hw_3_5(A,x,N_year)
IPython 5.1.0 -- An enhanced Interactive Python.
? -> Introduction and overview of IPython‘s features.
%quickref -> Quick reference.
help -> Python‘s own help system.
object? -> Details about ‘object‘, use ‘object??‘ for extra details.
PyDev console: using IPython 5.1.0
Running C:/Projects/python3_5/homework/hw_3-5_demo.py
[ 2400. 4800.]
虽然求解都类似,但是这是个收敛的差分方程,把100年换成1000年结果还是2400和4800。
实践:使用PCA给Mnist图片降维
再写程序的时候发现作业资料给的数据载入包并不能用(也许是python版本的问题),对debug不是很感兴趣,所以索性使用了tensorflow中提供的Mnist数据集调用方法了。
思路有一点偏差,我看了答案,其原意是把N幅数据组成N*(28*28)的二维矩阵,对这个矩阵进行降维,然后在绘图时还原整个矩阵,再在矩阵中进行子图分割;我理解成把每幅小图独自降维保存了,不过除了使我的程序麻烦了一点之外没什么其他影响:
import os import numpy as np import matplotlib.pyplot as plt from tensorflow.examples.tutorials.mnist import input_data print(os.getcwd()) os.environ[‘TF_CPP_MIN_LOG_LEVEL‘] = ‘3‘ batch_size = 10 mnist = input_data.read_data_sets(‘../../Mnist_data/‘, one_hot=True) batch_xs,batch_ys = mnist.train.next_batch(batch_size) print(batch_xs.shape) # write Your code def pca(data,k=0.5): # data 原始的图片 # k是保留特征值的百分比 # return 返回降维后重建的图片 res = np.empty_like(data).reshape(batch_size, 28, 28) rec = np.empty_like(res) data = http://www.mamicode.com/data.reshape(batch_size, 28, 28)>
保留占比50%的特征值后压缩(上行是原图,下行是压缩后的图片):
保留占比90%的特征值后压缩(上行是原图,下行是压缩后的图片):
直观来看90%对50%似乎提升不大,不过查看90%和50%保留的特征个数就卡以发现,其实两者的特征数目相差不大,基本上都在5个以下(总数应该是28个左右),也就是说PCA压缩是有道理的——实际上图片的大量信息被保存在极少的几个特征上了。
『科学计算』线性代数部分作业