Bidiagonalizace¶
příprava
In [1]:
import numpy as np;
import math;
matice pro příklad:
In [4]:
A = np.array([[3, -5, 15, -4],
[7, 7, 9, -1],
[0, 8, -1, 6],
[-10, 7, 8, -10]])
funkce na vytvoření Householderova reflektoru
In [10]:
def make_householder(x,y):
w = (x - y) / np.linalg.norm(x - y);
return (np.identity(len(x)) - 2 * np.dot(w[:,None],w[None,:]))
ukázka použití Householderova reflektoru.
In [11]:
x = A[:,0]
y = np.array([np.linalg.norm(x),0,0,0])
H = make_householder(x,y)
np.dot(H,A)
Out[11]:
pro zpřehlednění výsledku, zaokrouhlíme
In [13]:
np.around(np.dot(H,A),2)
Out[13]:
pokračujeme pro submatici
In [15]:
A1 = np.transpose(np.dot(H,A))[1:,:]
A1
Out[15]:
In [20]:
x = A1[:,0]
y = np.array([np.linalg.norm(x),0,0]);
W1 = make_householder(x,y)
WT = np.pad(W1,((1,0),(1,0)),'constant')
WT[0,0] = 1
WT
Out[20]:
In [21]:
np.around(np.dot(np.dot(H,A),WT),2)
Out[21]:
uděláme bidiagonalizační proceduru
In [25]:
def bidiagonalize_step(A):
x = A[:,0]
y = np.zeros(np.size(x))
y[0] = np.linalg.norm(x);
H = make_householder(x,y)
return(H)
def bidiagonalize(A):
colnum = np.size(A,1)
H = np.identity(colnum)
W = np.identity(colnum)
for i in range(colnum-1):
xlist = [0]
if (i!=colnum-2): xlist = [0,1]
for col in xlist:
A1 = A
if (col): A1 = np.transpose(A);
A1 = A1[i+col:,i:]
HH = bidiagonalize_step(A1)
HH = np.pad(HH,((i+col,0),(i+col,0)),'constant')
for k in range(i+col):
HH[k,k]=1;
if (col):
W = np.dot(W,HH)
A = np.dot(A,HH)
else:
H = np.dot(HH,H)
A = np.dot(HH,A)
return [np.transpose(H),A,np.transpose(W)]
testneme
In [30]:
[H,B,W] = bidiagonalize(A)
In [31]:
np.around(H,2)
Out[31]:
In [32]:
np.around(B,2)
Out[32]:
In [33]:
np.around(W,2)
Out[33]:
In [34]:
np.around(np.dot(np.dot(H,B),W),2)
Out[34]:
In [35]:
A
Out[35]:
In [41]:
T = np.dot(B,np.transpose(B))
np.around(T)
Out[41]:
QR algorithmus pro tridiagonální matici¶
Budeme potřebovat QR rozklad tridiagonalní matice
In [37]:
def QRdecomp(T):
siz = np.size(T,0);
Q = np.identity(siz)
for i in range(siz-1):
alpha = math.atan(T[i+1,i]/T[i,i])
rot = np.identity(siz)
s = math.sin(alpha)
c = math.cos(alpha)
rot[i,i] = c
rot[i+1,i+1] = c
rot[i,i+1] = s
rot[i+1,i] = -s
T = np.dot(rot,T)
Q = np.dot(rot,Q)
return([np.transpose(Q),T])
test
In [42]:
[Q,R]=QRdecomp(T)
In [43]:
np.around(Q)
Out[43]:
In [44]:
np.around(R)
Out[44]:
In [45]:
np.around(np.dot(Q,R))
Out[45]:
In [46]:
np.around(np.dot(R,Q))
Out[46]:
nyní můžeme vytvořit QR algoritmus
In [49]:
def QR(T):
siz = np.size(T,0);
QQ = np.identity(siz)
for i in range(100):
[Q,R] = QRdecomp(T)
T = np.dot(R,Q)
QQ = np.dot(QQ,Q)
return(QQ,T)
In [50]:
[Q,TT] = QR(T)
In [51]:
np.around(Q)
Out[51]:
In [52]:
np.around(TT)
Out[52]:
In [54]:
BB = np.dot(np.dot(Q,TT),np.transpose(Q));
np.around(BB)
Out[54]:
In [55]:
AA = np.dot(np.dot(H,BB),np.transpose(H));
np.around(AA,4)
Out[55]:
In [56]:
np.dot(A,np.transpose(A))
Out[56]:
Wilkinsonovův posun
In [69]:
def QR_wilkshift(T):
I = np.identity(np.size(T,0))
n = np.size(T,0);
for i in range(5):
mu=wilkshift(T[n-2:n,n-2:n]);
[Q,R]=QRdecomp(T-mu*I)
T= np.dot(R,Q) +mu*I
print(np.around(T,2))
return(T)
def wilkshift(T):
l,ll = np.linalg.eig(T);
if (abs(l[0]-T[1,1]) < abs(l[1]-T[1,1])):
return(l[0]);
else:
return(l[1]);
In [71]:
QR_wilkshift(T)
Out[71]:
Toto není celé, .. Dodělejte za domácí úkol.
In [72]:
U = np.dot(H,Q)
Sigma = np.diag(np.sqrt(np.diag(TT)))
In [73]:
AA
Out[73]:
In [74]:
U_Sigma = np.dot(U,Sigma)
In [78]:
np.around(np.dot(U_Sigma,np.transpose(U_Sigma)))
Out[78]:
In [79]:
[Q,TT2] = QR(np.dot(np.transpose(B),B))
In [80]:
Q
Out[80]:
In [82]:
np.diag(TT2)
Out[82]:
In [86]:
V = np.dot(np.transpose(Q),W)
In [87]:
V
Out[87]:
In [91]:
np.around(np.dot(U,np.dot(Sigma,V)))
Out[91]:
In [95]:
[_U,_Sigma,_V] = np.linalg.svd(A, full_matrices=True)
print(_U)
print(_Sigma)
print(_V)
In [96]:
U
Out[96]:
In [97]:
V
Out[97]:
In [102]:
my_data = np.genfromtxt('iris.csv', delimiter=',')[:,:-1]
In [103]:
my_data
Out[103]:
In [146]:
[_U,_Sigma,_V] = np.linalg.svd(my_data,full_matrices=False)
In [147]:
_Sigma
Out[147]:
In [148]:
np.shape(_U)
Out[148]:
In [149]:
_Sigma
Out[149]:
In [150]:
_Sigma = np.diag(_Sigma[:2])
_Sigma
Out[150]:
In [151]:
transformed_data = np.dot(_U[:,0:2],_Sigma)
In [130]:
from matplotlib import pyplot as plt
In [153]:
plt.figure(figsize=(5, 4))
plt.scatter(transformed_data[:,0],transformed_data[:,1])
plt.show()
