Untitled
unknown
python
a year ago
3.3 kB
1
Indexable
Never
import numpy as np from numpy import array , eye , outer import scipy from scipy.linalg import norm import math from matplotlib.image import imread from matplotlib import pyplot as plt def proj(v,U): vp = np.zeros(np.size(v)) for u in U: vp += (np.dot(u,v)/np.dot(u,u))*u return vp def gramschmidt(A): n,m = np.shape(A) u1 = A[:,0] e1 = u1/math.sqrt(np.dot(u1,u1)) U = [e1] # A[:,1:] har alla kolumner förutom första for v in A[:,1:].T: u = v - proj(v,U) e = u/math.sqrt(np.dot(u,u)) U.append(e) Q = np.zeros(np.shape(A)) for i,e in enumerate(U): Q[:,i] = e return Q def tests(Q): m,n = np.shape(Q) print(f'2Norm of (Q^TQ) = {norm(Q,ord=2)}') print(f'2Norm of (Q^TQ) - I = {norm(Q.T@Q-np.identity(n),ord=2)}') # eigenvärden också av Q^TQ print(f'Eigenvalues of (Q^TQ) = {np.linalg.eigvals(Q.T@Q)}') def Householder_step(A): m,n = np.shape(A) a = A[:,0] ahat=array([norm(a)]+(m-1)*[0.]) v = a - ahat v=v/norm(v) H = eye(m)-2*outer(v,v) return H def Givens_step(A): m,n = np.shape(A) a = A[:,0] # kanske behöver börja nerifrån om fel a0 = a[0] G = np.eye(m) for i in range(m-1,0,-1): ahat = a[i] theta = math.atan(ahat/a0) Gi = np.block([[math.cos(theta), np.zeros( (1,i-1) ), -math.sin(theta), np.zeros( ( 1, m-(i+1) ) )], [np.zeros((i-1,1)), np.eye(i-1), np.zeros((i-1,m-i)) ], [math.sin(theta), np.zeros((1,i-1)), math.cos(theta), np.zeros( (1,m-(i+1)) )], [np.zeros((m-(i+1),i+1)), np.eye(m-(i+1))] ]) G = Gi@G return G def QRH(A): m,n = np.shape(A) tempA = A H = np.eye(m) for i in range(n): # byt ut denna för reflection vs rotation Hi = Householder_step(tempA) tempA = Hi@tempA # ta bort första rad första kolumn tempA = np.delete(tempA,0,0) tempA = np.delete(tempA,0,1) H = H@np.block([[np.eye(i), np.zeros((i,m-i))],[np.zeros((m-i,i)), Hi]]) R = H.T@A Q = H return Q,R def QRG(A): m,n = np.shape(A) tempA = A H = np.eye(m) for i in range(n): # byt ut denna för reflection vs rotation Hi = Givens_step(tempA) tempA = Hi@tempA # ta bort första rad första kolumn tempA = np.delete(tempA,0,0) tempA = np.delete(tempA,0,1) H = H@np.block([[np.eye(i), np.zeros((i,m-i))],[np.zeros((m-i,i)), Hi]]) R = H.T@A Q = H return Q,R n = 4 m = n + 2 A = np.random.rand(m,n) Q1,_ = scipy.linalg.qr(A) Q2,R2 = QRG(A) #print(f'Q1Q1^T = {Q1@Q1.T}') print(f'Det (Q1) = {np.linalg.det(Q1)}') #print(f'Q2Q2^T = {Q2@Q2.T}') print(f'Det (Q2) = {np.linalg.det(Q2)}') print(f'A = {A}') print(f'Q2R2 = {Q2@R2}') R2[R2<10**(-10)] = 0 print(f'R2 = {R2}') # B = imread('kvinnagr.jpeg') # print(B) # U,s,VT = scipy.linalg.svd(B) # m,n = np.shape(B) # S = scipy.linalg.diagsvd(s,m,n) #S = np.zeros((m,n)) #for i in range(0,min(m,n)): # S[i,i] = s[i] # # b = np.round(U@S@VT) # b = b.astype(int) # # print(f'Difference of B and USVT = {norm(B-b,ord=2)}') # plt.imshow(b,cmap='gray', vmin=0, vmax=255.) #Q2 = gramschmidt(A) #tests(Q1) #tests(Q2)