Python Numpy matrix计算 初试

准备

• numpy.random.normal(local,scale,size) size是二元组，代表生成的随机高斯矩阵大小，local为均值，scale为方差
• numpy.dot(A,B) 矩阵乘法 A * B
• numpy.identity(size) size为单位矩阵大小
• from scipy.linalg import toeplitz。 Toeplitz 矩阵

scipy.linalg.toeplitz(c, r=None)
 c : array_like
First column of the matrix. Whatever the actual shape of c, it will be converted to a 1-D array.

r : array_like, optional
First row of the matrix. If None, r = conjugate(c) is assumed; in this case, if c[0] is real, the result is a Hermitian matrix. r[0] is ignored; the first row of the returned matrix is [c[0], r[1:]]. Whatever the actual shape of r, it will be converted to a 1-D array.  

• http://ergodic.ugr.es/cphys/LECCIONES/FORTRAN/power_method.pdf 题四的计算方法

  题目

  “`
  Exercise 9.1: Matrix operations
  Calculate A + A, AA>;A>A and AB. Write a function that computes A(B �� I) for any .

Exercise 9.2: Solving a linear system
Generate a vector b with m entries and solve Bx = b.

Exercise 9.3: Norms
Compute the Frobenius norm of A: kAkF and the in nity norm of B: kBk1. Also nd the largest and
smallest singular values of B.

Exercise 9.4: Power iteration
Generate a matrix Z, n n, with Gaussian entries, and use the power iteration to nd the largest
eigenvalue and corresponding eigenvector of Z. How many iterations are needed till convergence?
Optional: use the time.clock() method to compare computation time when varying n.

Exercise 9.5: Singular values
Generate an nn matrix, denoted by C, where each entry is 1 with probability p and 0 otherwise. Use
the linear algebra library of Scipy to compute the singular values of C. What can you say about the
relationship between n, p and the largest singular value?

Exercise 9.6: Nearest neighbor
Write a function that takes a value z and an array A and nds the element in A that is closest to z. The
function should return the closest value, not index.
Hint: Use the built-in functionality of Numpy rather than writing code to nd this value manually. In
particular, use brackets and argmin.
“`

解答

import numpy as np
from scipy.linalg import toeplitz


""" #1 """ 
#a = np.random.normal(0,1,(2,2))
#print(a)

def calcu_lam(A,B,lam):
    return np.dot(A,B- lam*np.identity(B.shape[0]));

n = 200
m = 500
np.random.normal()
A = np.random.normal(0,1,(n,m))
B = toeplitz(A[0,0:])

print('##1\n',A)
print('##1\n',B)
print('##1\n',A+A)
print('##2\n',np.dot(A,A.transpose()))
print('##3\n',np.dot(A.transpose(),A))
print('##4\n',np.dot(A,B))

print('##5\n',calcu_lam(A,B,10))

"""" #2 """
b = np.random.normal(0,1,m)
x = np.linalg.solve(B,b)
print("####2\n",x)

""" #3 """
print(np.linalg.norm(A,'fro'))
print(np.linalg.norm(B,np.inf))
print(svdvals(B).max())
print(svdvals(B).min())

# -*- coding: utf-8 -*- 
import numpy as np   
import time
n = 0   
""" #4 """
for n in [5, 20, 100, 200]:  
    s_time = time.clock()  

    A = np.random.normal(0,1,(n,n))

    v = u = np.ones([n, 1])  
    epson = 0  
    it = 0

    while (True):
        v = np.dot(A,u)
        pre_e = epson
        epson = v.max()
        u = v / epson

        if (abs(pre_e - epson) < 1e-10):
            break;
        it+=1

    print('epson:',str(epson))
    print('iter :',str(it))

    e_time = time.clock()

    print("n: ",n,"time:",e_time - s_time)

# -*- coding: utf-8 -*-

import numpy as np   
import scipy

""" #5 """
for n in [10, 50, 100, 200]:  
    for p in [0, 0.2, 0.4, 0.6, 0.8, 1]:  
        #generate metrix C   
        C = np.zeros([n*n, 1])  
        for i in range(len(C)):
            C[i] = 1 if np.random.sample(1) < p else 0  

        sv = scipy.linalg.svdvals(C).max()  
        print('n:' ,n, " p:" ,p, '\n\t', 'largest value:',sv) 


""" #6 """ 
def nearest(A,z):
    return A[np.argmin(np.abs(A-z))]

问题

第四题按照pdf上做的一直没收敛。。。