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)
```