为了显示方便,n=4,m=10
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
n=4
m=10
A = np.random.normal(0,1,n*m)
A = A.reshape(n,m)
B = toeplitz(list(range(1,m+1)))
print("A=",A)
print("B=",B)
print("A+A=",A+A) #A+A
print("A*A.T=",A@A.T) #A*A.T
print("A.T*A=",A.T@A) #A.T*A
print("AB=",A@B) #AB
def cal(A,B,lamb):
return A@(B-lamb*np.eye(m)) #A(B−λI)
print("A(B−λI)=",cal(A,B,random.random()))
输出:
A= [[-0.14127543 -0.18536895 -0.3919833 0.3917631 0.27335622 -1.95794558
1.1441752 -0.63235459 0.1720679 -0.97279904]
[-1.38390472 -0.16494005 -0.85358231 -0.31362912 -0.44616255 0.26444861
0.08132372 1.16820247 -1.91757158 1.07926125]
[-0.2395814 -0.35344047 0.27283204 1.05136257 -1.504381 -0.19087067
-1.41278342 0.17866861 -0.49925697 0.81361346]
[ 0.38664661 -0.01007491 1.69917109 1.72546369 -0.28064128 0.15060185
0.84807271 1.42510812 -1.42309945 -0.50420029]]
B= [[ 1 2 3 4 5 6 7 8 9 10]
[ 2 1 2 3 4 5 6 7 8 9]
[ 3 2 1 2 3 4 5 6 7 8]
[ 4 3 2 1 2 3 4 5 6 7]
[ 5 4 3 2 1 2 3 4 5 6]
[ 6 5 4 3 2 1 2 3 4 5]
[ 7 6 5 4 3 2 1 2 3 4]
[ 8 7 6 5 4 3 2 1 2 3]
[ 9 8 7 6 5 4 3 2 1 2]
[10 9 8 7 6 5 4 3 2 1]]
A+A= [[-0.28255086 -0.37073789 -0.78396659 0.78352621 0.54671244 -3.91589116
2.28835041 -1.26470919 0.3441358 -1.94559809]
[-2.76780944 -0.32988011 -1.70716461 -0.62725824 -0.89232511 0.52889721
0.16264745 2.33640494 -3.83514316 2.15852249]
[-0.4791628 -0.70688094 0.54566408 2.10272513 -3.00876199 -0.38174135
-2.82556684 0.35733721 -0.99851394 1.62722692]
[ 0.77329322 -0.02014982 3.39834218 3.45092739 -0.56128256 0.3012037
1.69614542 2.85021623 -2.84619889 -1.00840059]]
A*A.T= [[ 6.95467871 -2.22745554 -2.24005637 -0.09962909]
[-2.22745554 9.25155372 2.37724554 1.55859419]
[-2.24005637 2.37724554 6.60081763 1.93879692]
[-0.09962909 1.55859419 1.93879692 11.14503533]]
A.T*A= [[ 2.14204587 0.33523173 1.82826745 0.79394409 0.83074066 0.01459668
0.39219294 -1.0191372 2.19880351 -1.74603612]
[ 0.33523173 0.18658854 0.09990256 -0.40986868 0.55745488 0.38526826
0.26528249 -0.15297102 0.47518353 -0.28017082]
[ 1.82826745 0.09990256 3.84387337 3.33284709 -0.61361549 0.74557598
0.53765399 1.72096451 -0.9849453 -1.17466007]
[ 0.79394409 -0.40986868 3.33284709 4.33442976 -1.81886575 -0.79080588
0.40071121 2.03270229 -2.31159037 -0.73417109]
[ 0.83074066 0.55745488 -0.61361549 -1.81886575 2.61570636 -0.40832655
2.16384413 -1.36279608 2.05303761 -1.82993184]
[ 0.01459668 0.38526826 0.74557598 -0.79080588 -0.40832655 3.96259649
-1.82134659 1.72756672 -0.96302661 1.95886828]
[ 0.39219294 0.26528249 0.53765399 0.40071121 2.16384413 -1.82134659
4.03093476 0.32765339 -0.46061807 -2.60234112]
[-1.0191372 -0.15297102 1.72096451 2.03270229 -1.36279608 1.72756672
0.32765339 3.82742496 -4.4661919 1.30277685]
[ 2.19880351 0.47518353 -0.9849453 -2.31159037 2.05303761 -0.96302661
-0.46061807 -4.4661919 5.98115767 -1.92562321]
[-1.74603612 -0.28017082 -1.17466007 -0.73417109 -1.82993184 1.95886828
-2.60234112 1.30277685 -1.92562321 3.02732762]]
AB= [[-15.73079285 -13.71297924 -12.06590353 -11.20279441 -9.55615909
-7.36281132 -9.08535471 -8.51954769 -9.21844986 -9.57321624]
[ -2.72381528 -3.00507044 -3.6162057 -5.93450558 -8.88006369
-12.71794692 -16.02693294 -19.1732715 -19.98320512 -24.62828191]
[ -9.4069582 -8.00228375 -7.30449023 -6.06103263 -2.7148499
-2.37742916 -2.42174976 -5.29163721 -7.80418744 -11.31525162]
[ 11.35374551 8.10999059 4.84608584 4.98052327 8.56588809
11.58997036 14.91525632 19.93668771 27.80833533 32.83378406]]
A(B−λI)= [[-15.62625848 -13.57581863 -11.77586213 -11.49267289 -9.7584244
-5.91406263 -9.93196778 -8.05164762 -9.34576859 -8.85341005]
[ -1.69981841 -2.88302584 -2.98461192 -5.70244102 -8.54993327
-12.91362119 -16.08710705 -20.03766314 -18.56433049 -25.426863 ]
[ -9.229684 -7.74076146 -7.50636768 -6.83897056 -1.60170864
-2.23619763 -1.37638458 -5.42384002 -7.4347707 -11.91727113]
[ 11.0676529 8.11744534 3.58881296 3.70379562 8.77354386
11.47853507 14.28773927 18.88220308 28.86133367 33.20685853]]
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
n=4
m=10
A = np.random.normal(0,1,n*m)
A = A.reshape(n,m)
B = toeplitz(list(range(1,m+1)))
b = np.array(range(0,m))
print("A=",A)
print("B=",B)
print("b=",b)
print("B\\b=",np.linalg.solve(B,b))
print("Bx=",B@np.linalg.solve(B,b))
输出:
A= [[ 1.10072277 -1.74155593 -0.87292138 2.36295768 1.06313563 -1.09282835
-1.65235618 -1.29662255 1.02735804 0.92605323]
[-0.20442559 -1.13772172 -0.09905539 -2.16642607 0.05472793 1.26894958
0.02942699 0.06011498 1.61455627 -0.53479065]
[ 0.09988037 0.78840479 -0.22336113 1.17640885 0.65951768 -0.35056981
0.27483197 -1.74387682 -0.21052275 -0.51640518]
[-0.82507688 -0.2172729 -1.12598322 1.04701727 0.46405799 -1.51862459
1.44115554 -1.92401456 0.08674978 -0.15373378]]
B= [[ 1 2 3 4 5 6 7 8 9 10]
[ 2 1 2 3 4 5 6 7 8 9]
[ 3 2 1 2 3 4 5 6 7 8]
[ 4 3 2 1 2 3 4 5 6 7]
[ 5 4 3 2 1 2 3 4 5 6]
[ 6 5 4 3 2 1 2 3 4 5]
[ 7 6 5 4 3 2 1 2 3 4]
[ 8 7 6 5 4 3 2 1 2 3]
[ 9 8 7 6 5 4 3 2 1 2]
[10 9 8 7 6 5 4 3 2 1]]
b= [0 1 2 3 4 5 6 7 8 9]
B\b= [ 9.09090909e-01 5.04646829e-17 -1.11022302e-16 1.11022302e-16
-2.77555756e-16 1.11022302e-16 3.60822483e-16 5.55111512e-17
-4.71844785e-16 -9.09090909e-02]
Bx= [4.4408921e-16 1.0000000e+00 2.0000000e+00 3.0000000e+00 4.0000000e+00
5.0000000e+00 6.0000000e+00 7.0000000e+00 8.0000000e+00 9.0000000e+00]
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
n=4
m=10
A = np.random.normal(0,1,n*m)
A = A.reshape(n,m)
B = toeplitz(list(range(1,m+1)))
print("A=",A)
print("B=",B)
print("||A||F=",np.linalg.norm(A,'fro')) #2-范数
print("||B||∞=",np.linalg.norm(B,np.inf)) #无穷范数
u,sigma,v = np.linalg.svd(B) #奇异值分解
print("Max singular value = ",max(sigma))
print("Min singular value = ",min(sigma))
输出:
A= [[ 1.76739889 -0.61062597 0.02890403 -0.73851563 -0.67586851 0.63724716
0.08316141 -1.5604136 1.55553185 0.22876358]
[-0.35941143 -0.80682195 0.09212975 -1.36546868 1.06930886 1.71811841
-0.82373042 -0.83853979 -1.7913538 -0.26047562]
[ 0.90643659 -0.74680685 -1.00215034 0.45642739 1.76741504 1.23191851
2.03571596 0.2671933 1.85416188 0.49919575]
[-1.46066207 0.28280865 -0.41617868 1.73892589 -0.66836749 0.06630803
-0.03678402 -0.0268019 1.28013338 0.27449861]]
B= [[ 1 2 3 4 5 6 7 8 9 10]
[ 2 1 2 3 4 5 6 7 8 9]
[ 3 2 1 2 3 4 5 6 7 8]
[ 4 3 2 1 2 3 4 5 6 7]
[ 5 4 3 2 1 2 3 4 5 6]
[ 6 5 4 3 2 1 2 3 4 5]
[ 7 6 5 4 3 2 1 2 3 4]
[ 8 7 6 5 4 3 2 1 2 3]
[ 9 8 7 6 5 4 3 2 1 2]
[10 9 8 7 6 5 4 3 2 1]]
||A||F= 6.628749769033752
||B||∞= 55.0
Max singular value = 44.076292527938
Min singular value = 0.5125428154684579
题目是想算出主特征值和主特征向量。算法是用幂法求解。
算法:
输入:A;输出:λ和X
u:=随机向量
While 不满足判停准则 do
v:=Au
λ:=max(v)
u:=v/λ
End
x:=u
max(v)为向量v的绝对值最大的分量,即max(v)=vj,其中j满足|vj|=max(1≤k≤n)|vk|,若j的值不唯一,则取最小的那个,并且称u=v/max(v)为向量v的规格化向量。
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
import time
n=100
Z = np.random.normal(0,1,n*n)
Z = Z.reshape(n,n)
def powerIteration(A):
begintime = time.clock()
u = np.random.normal(0,1,n)
lamb=0
lastlamb=0
count=0
while abs(lastlamb-lamb)>=1e-6 or count==0:
v = A@u
lastlamb=lamb
j = np.argmax(abs(v))
lamb = v[j]
u = v/lamb
count += 1
x = u
endtime=time.clock()
return x,lamb,count,endtime-begintime
x,lamb,count,time=powerIteration(Z)
print("主特征值=",lamb)
#print("主特征向量=",x)
print("迭代次数=",count)
print("花费时间=",time)
可以通过修改n来测试算法
输出:
精度要求6位小数
由于数据是随机的,因此时间上没有可比性
矩阵大小= 10
主特征值= -3.6392824996637034
主特征向量= [ 1. 0.14361591 -0.19157489 0.03119712 0.52552813 -0.36625306
0.94472046 -0.72621135 -0.14788755 -0.44033934]
迭代次数= 69
花费时间= 0.00036187682913724775
矩阵大小= 50
主特征值= 7.766544319679567
迭代次数= 238
花费时间= 0.0013120010366427944
矩阵大小= 100
主特征值= 10.500415332054505
迭代次数= 320
花费时间= 0.0225031288913611
矩阵大小= 1000
主特征值= 32.54980745398664
迭代次数= 447
花费时间= 0.14958747621726515
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
import time
for n in (10,50,100):
for p in (0.1,0.5,0.9):
print("n=%d,p=%.1f"%(n,p))
C = np.random.binomial(1,p,n*n)
C = C.reshape(n,n)
u,sigma,v = np.linalg.svd(C)
print("最大特异值=",max(sigma))
输出:
n=10,p=0.1
最大特异值= 1.9021130325903073
n=10,p=0.5
最大特异值= 4.859185150253102
n=10,p=0.9
最大特异值= 9.101742974321443
n=50,p=0.1
最大特异值= 6.189969662651178
n=50,p=0.5
最大特异值= 25.767586355786325
n=50,p=0.9
最大特异值= 44.688202529805025
n=100,p=0.1
最大特异值= 10.733560377216394
n=100,p=0.5
最大特异值= 49.55564032221308
n=100,p=0.9
最大特异值= 90.75654035220786
猜测关系是最大特异值=n*p
argmin返回最小值的下标,因此求出两个数的距离的绝对值再调用argmin就知道是哪个数了。
代码:
import numpy as np
from scipy.linalg import toeplitz
import random
import time
z = 5
A = np.random.normal(5,5,10)
print("A=",A)
print("z=",z)
print("|A-z|=",abs(A-z))
print(np.argmin(abs(A-z)))
输出:
A= [ 3.40970621 3.3574148 7.01400765 -1.56122245 7.15815137 7.22773429
11.59661275 3.49448966 -2.19594911 1.23482191]
z= 5
|A-z|= [1.59029379 1.6425852 2.01400765 6.56122245 2.15815137 2.22773429
6.59661275 1.50551034 7.19594911 3.76517809]
7