Exercise9-Numpy

最新推荐文章于 2024-05-29 18:45:06 发布

ponyooo3

最新推荐文章于 2024-05-29 18:45:06 发布

阅读量193

点赞数

分类专栏： python作业文章标签： python numpy exercise

本文链接：https://blog.csdn.net/ponyooo3/article/details/80394777

版权

python作业专栏收录该内容

21 篇文章 0 订阅

订阅专栏

本文通过一系列Exercise介绍了如何使用numpy库进行矩阵操作，包括矩阵的生成、线性方程组求解、范数计算、幂法求特征值、奇异值分解以及最近邻计算。通过实例代码详细解析了numpy的使用方法和内置函数。

摘要由CSDN通过智能技术生成

Premise:

Premise

分析：

要生成两个矩阵，A为高斯（正态分布）矩阵，可以numpy.random.normal获得；B为Toeplitz矩阵，可用scipy.linalg 中的toeplitz获得。

内置函数等介绍：

numpy.random.normal(loc=0.0, scale=1.0, size=None)
对应于高斯分布的概率密度函数，第一，二个参数分别为均值和标准差，最后为输出的形状。例如numpy.random.normal(loc=0, scale=1, size=(2，3))正态分布的矩阵。
scipy.linalg.toeplitz(c, r=None)
用于生成toeplitz矩阵，第一个参数为矩阵的第一列，第二个为矩阵的第一行（默认为c的转置）。

代码实现：

import numpy as np
import time
from scipy.linalg import toeplitz
#N(0,1)gaussian
mu,sigma = 0,1 #mean and standard deviation
n,m = 200,500 #the dimension of A
A = np.random.normal(loc=mu, scale=sigma, size=(n,m))
c = list(range(1,m+1))
B = toeplitz(c,c)

Exercise 9.1: Matrix operations:

9-1

分析：

矩阵之间的简单运算，使用简单的运算符或numpy的相关方法即可实现，I为单位矩阵，可用方法numpy.eye获得。

内置函数等介绍：

基本运算
获得矩阵A的转置：A.T
共轭矩阵：A.I
矩阵相加减：A+A A-A
逆矩阵：numpy.linalg.inv(a)
矩阵相乘：a.dot(b.T) numpy.dot(a,b) a*b
求迹： np.trace(a)
生成对角矩阵
numpy.eye(N,M=None, k=0, dtype=(type ‘float’))
第一个参数为输出方阵（行数=列数）的规模，即行数或列数，第二个参数k默认情况下输出的是对角线全“1”，其余全“0”的方阵，如果k为正整数，则在右上方第k条对角线全“1”其余全“0”，k为负整数则在左下方第k条对角线全“1”其余全“0”。例如：

>>> np.eye(2, dtype=int)
array([[1, 0],
       [0, 1]])
>>> np.eye(3, k=1)
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])

代码实现：

其中函数以lambda为2调用

#9.1
print("---9.1:Matrix operations")
print("calculate A+A:")
dA = A+A
print(dA)
print("calculate AA(T):")
mA = np.dot(A,A.T)
print(mA)
print("calculate A(T)A:")
mA2 = np.dot(A.T,A)
print(mA2)
print("calculate AB:")
mAB = np.dot(A,B)
print(mAB)

def func(A, B, lamb):
    C = B-lamb*(np.eye(m))
    print(np.dot(A,C))
print("calculate A(B-I) for  = 2:")
func(A,B,2)

Exercise 9.2: Solving a linear system

这里写图片描述

分析：

numpy.linalg中的函数solve可以求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量。

内置函数等介绍：

np.linalg.solve(a,b)
solve函数有两个参数a和b。a是一个N*N的二维数组，而b是一个长度为N的一维数组，solve函数找到一个长度为N的一维数组x，使得a和x的矩阵乘积正好等于b，数组x就是多元一次方程组的解。
Python列表和Numpy数组的区别：
Numpy使用ndarray对象来处理多维数组，该对象是一个快速而灵活的大数据容器。使用Python列表可以存储一维数组，通过列表的嵌套可以实现多维数组，那么为什么还需要使用Numpy呢？Numpy是专门针对数组的操作和运算进行了设计，所以数组的存储效率和输入输出性能远优于Python中的嵌套列表，数组越大，Numpy的优势就越明显。通常Numpy数组中的所有元素的类型都是相同的，而Python列表中的元素类型是任意的，所以在通用性能方面Numpy数组不及Python列表，但在科学计算中，可以省掉很多循环语句，代码使用方面比Python列表简单的多。

代码实现：

用numpy.ones生成全1的向量b

#9.2
print("\n---9.2:Solving a linear system ")
b = np.ones(m)
x = np.linalg.solve(B,b)
print("The solution to Bx=b, b is a unit vector is:")
print(x)

Exercise 9.3: Norms

9-3

分析：

使用numpy中计算范数的方法，并根据范数与特征值的对应关系求特征值。（谱范数对应最大）

内置函数等介绍：

求范数
列范数numpy.linalg.norm(b,1)
谱范数,为x里最大值开平方numpy.linalg.norm(b,2)
无穷范数，行范数numpy.linalg.norm(b,np.inf)
Frobenius范数numpy.linalg.norm(b,”fro”)
2-范数：
║A║2 = A的最大奇异值 = (max{ λi(AH*A) }) 1/2 （谱范数，即A^H*A特征值λi中最大者λ1的平方根，其中AH为A的转置共轭矩阵）

代码实现：

#9.3
print("\n---9.3:Norms ")
print("The frobenius norm of A is: ")
norm_A = np.linalg.norm(A,'fro')
print(norm_A)
print("The infinity norm of B is: ")
norm_B = np.linalg.norm(B,np.inf)
print(norm_B)
print("The biggest singular value of B is: ")
lar_B = np.linalg.norm(B,2)
print(lar_B)
print("The smallest singular value of B is: ")
sma_B = np.linalg.norm(B,-2)
print(sma_B)

Exercise 9.4: Power iteration

这里写图片描述

分析：

参考数值计算中幂法求主特征值和主特征向量的算法，判断误差小于0.0005时停止迭代。Z的生成同上面的A矩阵。使用time.clock()记录开始的时间和结束的时间，相减。

内置函数等介绍：

代码实现：

#9.4
print("\n---9.4:Power iteration ")
def power_iteration(n):
    Z = np.random.normal(loc=mu, scale=sigma, size=(n,n))
    u = np.ones(n)
    lamb1 = 0
    ite_num = 0
    begin = time.clock()
    while True:
        v = Z*u
        temp = lamb1
        lamb1 = np.linalg.norm(b,1)
        u = v/lamb1
        ite_num += 1
        if abs(lamb1-temp)<0.0005:
            break;
    end = time.clock()

    print("the largest eigenvalue: ", end="")
    print(lamb1)
    print("corresponding eigenvector: ", end="")
    print(lamb1)
    print("n: ", end="")
    print(n)
    print("time taken: ", end="")
    print(end-begin)
    print("iteration num: ", end="")
    print(ite_num)
power_iteration(5)
power_iteration(50)
power_iteration(500)

Exercise 9.5: Singular values

这里写图片描述

分析：

已知numpy中实现二项分布的方法，而这个相当于对于每一个C的entry，有p的可能性是1，也就是0-1分布，适当地赋值该二项分布的方法参数即可，通过SVD分解的方法算出特征值，并求最大数。

内置函数等介绍：

np.random.binomial(n,p,size=N)
函数的返回值表示n中成功的次数，且以Cn^x*p^x*(1-p)^(n-x)的概率选择成功x次。这里令n=1，即为0-1分布，size为矩阵的（n,n）
numpy.linalg.svd(a, full_matrices=1, compute_uv=1)
奇异值分解，将矩阵A表示为u * np.diag(s) * v，s为一维向量（特征值）

代码实现：

#9.5
print("\n---9.5:Singular values ")
def cal_singu(n,p):
    C = np.random.binomial(1,p,(n,n))
    U, s, V = np.linalg.svd(C, full_matrices=True) #svd
    singu_value = max(s)
    print(singu_value)

Exercise 9.6: Nearest neighbor

这里写图片描述

分析：

根据提示分析，可以先将A中的每个值都减z，获得一个新矩阵，再使用numpy.argmin求其中最小值（绝对值最小），再加上z即为A中原值。

内置函数等介绍：

numpy.argmin(a, axis=None, out=None)
返回沿轴axis最小值的索引

代码实现：

#9.6
print("\n---9.6:Nearest neighbor ")
def neighbor(z, A):
    sub = A-z
    result = abs(sub)
    print("nearest neighbor:")
    print(np.argmin(result))
A = np.arange(6).reshape(2,3)

完整代码：

import numpy as np
import time
from scipy.linalg import toeplitz
#N(0,1)gaussian
mu,sigma = 0,1 #mean and standard deviation
n,m = 200,500 #the dimension of A
A = np.random.normal(loc=mu, scale=sigma, size=(n,m))
c = list(range(1,m+1))
B = toeplitz(c,c)
#9.1
print("---9.1:Matrix operations")
print("calculate A+A:")
dA = A+A
print(dA)
print("calculate AA(T):")
mA = np.dot(A,A.T)
print(mA)
print("calculate A(T)A:")
mA2 = np.dot(A.T,A)
print(mA2)
print("calculate AB:")
mAB = np.dot(A,B)
print(mAB)

def func(A, B, lamb):
    C = B-lamb*(np.eye(m))
    print(np.dot(A,C))
print("calculate A(B-I) for  = 2:")
func(A,B,2)
#9.2
print("\n---9.2:Solving a linear system ")
b = np.ones(m)
x = np.linalg.solve(B,b)
print("The solution to Bx=b, b is a unit vector is:")
print(x)
#9.3
print("\n---9.3:Norms ")
print("The frobenius norm of A is: ")
norm_A = np.linalg.norm(A,'fro')
print(norm_A)
print("The infinity norm of B is: ")
norm_B = np.linalg.norm(B,np.inf)
print(norm_B)
print("The biggest singular value of B is: ")
lar_B = np.linalg.norm(B,2)
print(lar_B)
print("The smallest singular value of B is: ")
sma_B = np.linalg.norm(B,-2)
print(sma_B)
#9.4
print("\n---9.4:Power iteration ")
def power_iteration(n):
    Z = np.random.normal(loc=mu, scale=sigma, size=(n,n))
    u = np.ones(n)
    lamb1 = 0
    ite_num = 0
    begin = time.clock()
    while True:
        v = Z*u
        temp = lamb1
        lamb1 = np.linalg.norm(b,1)
        u = v/lamb1
        ite_num += 1
        if abs(lamb1-temp)<0.0005:
            break;
    end = time.clock()

    print("the largest eigenvalue: ", end="")
    print(lamb1)
    print("corresponding eigenvector: ", end="")
    print(lamb1)
    print("n: ", end="")
    print(n)
    print("time taken: ", end="")
    print(end-begin)
    print("iteration num: ", end="")
    print(ite_num)
power_iteration(5)
power_iteration(50)
power_iteration(500)
#9.5
print("\n---9.5:Singular values ")
def cal_singu(n,p):
    C = np.random.binomial(1,p,(n,n))
    U, s, V = np.linalg.svd(C, full_matrices=True) #svd
    singu_value = max(s)
    print("largest singular value of C:")
    print(singu_value)
print("n=5, p=0.5:")
cal_singu(5,0.5)
print("n=5, p=0.8:")
cal_singu(5,0.8)
print("n=50, p=0.5:")
cal_singu(50,0.5)
print("n=50, p=0.8:")
cal_singu(50,0.8)
print("n=500, p=0.5:")
cal_singu(500,0.5)
print("n=500, p=0.8:")
cal_singu(500,0.8)
#9.6
print("\n---9.6:Nearest neighbor ")
def neighbor(z, A):
    sub = A-z
    result = abs(sub)
    print("nearest neighbor:")
    print(np.argmin(result))
A = np.arange(6).reshape(2,3)
print("z=3.3, a=[[0 1 2][3 4 5]]:")
neighbor(3.3, A)