mirror of
https://github.com/Visualize-ML/Book4_Power-of-Matrix.git
synced 2026-05-05 12:00:17 +08:00
Add files via upload
This commit is contained in:
Visualize-ML
108
Book4_Ch15_Python_Codes/Bk4_Ch15_01.py
Normal file
108
Book4_Ch15_Python_Codes/Bk4_Ch15_01.py
Normal file
@@ -0,0 +1,108 @@
|
||||
|
||||
###############
|
||||
# Authored by Weisheng Jiang
|
||||
# Book 4 | From Basic Arithmetic to Machine Learning
|
||||
# Published and copyrighted by Tsinghua University Press
|
||||
# Beijing, China, 2022
|
||||
###############
|
||||
|
||||
# Bk4_Ch15_01.py
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
def visualize(X_circle,X_vec,title_txt):
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
|
||||
plt.plot(X_circle[:,0], X_circle[:,1],'k',
|
||||
linestyle = '--',
|
||||
linewidth = 0.5)
|
||||
|
||||
plt.quiver(0,0,X_vec[0,0],X_vec[0,1],
|
||||
angles='xy', scale_units='xy',scale=1,
|
||||
color = [0, 0.4392, 0.7529])
|
||||
|
||||
plt.quiver(0,0,X_vec[1,0],X_vec[1,1],
|
||||
angles='xy', scale_units='xy',scale=1,
|
||||
color = [1,0,0])
|
||||
|
||||
plt.axvline(x=0, color= 'k', zorder=0)
|
||||
plt.axhline(y=0, color= 'k', zorder=0)
|
||||
|
||||
plt.ylabel('$x_2$')
|
||||
plt.xlabel('$x_1$')
|
||||
|
||||
ax.set_aspect(1)
|
||||
ax.set_xlim([-2.5, 2.5])
|
||||
ax.set_ylim([-2.5, 2.5])
|
||||
ax.grid(linestyle='--', linewidth=0.25, color=[0.5,0.5,0.5])
|
||||
ax.set_xticks(np.linspace(-2,2,5));
|
||||
ax.set_yticks(np.linspace(-2,2,5));
|
||||
plt.title(title_txt)
|
||||
plt.show()
|
||||
|
||||
|
||||
theta = np.linspace(0, 2*np.pi, 100)
|
||||
|
||||
circle_x1 = np.cos(theta)
|
||||
circle_x2 = np.sin(theta)
|
||||
|
||||
X_vec = np.array([[1,0],
|
||||
[0,1]])
|
||||
|
||||
X_circle = np.array([circle_x1, circle_x2]).T
|
||||
|
||||
# plot original circle and two vectors
|
||||
visualize(X_circle,X_vec,'Original')
|
||||
|
||||
A = np.array([[1.6250, 0.6495],
|
||||
[0.6495, 0.8750]])
|
||||
|
||||
# plot the transformation of A
|
||||
|
||||
visualize(X_circle@A.T, X_vec@A.T,'$A$')
|
||||
|
||||
|
||||
#%% SVD
|
||||
# A = U @ S @ V.T
|
||||
|
||||
U, S, V = np.linalg.svd(A)
|
||||
S = np.diag(S)
|
||||
V[:,0] = -V[:,0] # reverse sign of first vector of V
|
||||
U[:,0] = -U[:,0] # reverse sign of first vector of U
|
||||
|
||||
print('=== U ===')
|
||||
print(U)
|
||||
print('=== S ===')
|
||||
print(S)
|
||||
print('=== V ===')
|
||||
print(V)
|
||||
|
||||
# plot the transformation of V
|
||||
|
||||
visualize(X_circle@V, X_vec@V,'$V^T$')
|
||||
|
||||
# plot the transformation of V @ S
|
||||
|
||||
visualize(X_circle@V@S, X_vec@V@S,'$SV^T$')
|
||||
|
||||
# plot the transformation of V @ S @ U.T
|
||||
|
||||
visualize(X_circle@V@S@U.T, X_vec@V@S@U.T,'$USV^T$')
|
||||
|
||||
e1 = np.array([[1],
|
||||
[0]])
|
||||
|
||||
e2 = np.array([[0],
|
||||
[1]])
|
||||
|
||||
# Calculate step by step from e1 and e2
|
||||
VT_e1 = V.T@e1
|
||||
VT_e2 = V.T@e2
|
||||
|
||||
S_VT_e1 = S@VT_e1
|
||||
S_VT_e2 = S@VT_e2
|
||||
|
||||
U_S_VT_e1 = U@S_VT_e1
|
||||
U_S_VT_e2 = U@S_VT_e2
|
||||
225
Book4_Ch15_Python_Codes/Bk4_Ch15_02.py
Normal file
225
Book4_Ch15_Python_Codes/Bk4_Ch15_02.py
Normal file
@@ -0,0 +1,225 @@
|
||||
|
||||
###############
|
||||
# Authored by Weisheng Jiang
|
||||
# Book 4 | From Basic Arithmetic to Machine Learning
|
||||
# Published and copyrighted by Tsinghua University Press
|
||||
# Beijing, China, 2022
|
||||
###############
|
||||
|
||||
# Bk4_Ch15_02_A
|
||||
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
plt.rcParams['image.cmap'] = 'RdBu_r'
|
||||
import seaborn as sns
|
||||
PRECISION = 3
|
||||
|
||||
def svd(X):
|
||||
|
||||
full_matrices = True
|
||||
|
||||
U, s, Vt = np.linalg.svd(X,full_matrices = full_matrices)
|
||||
|
||||
# Put the vector singular values into a padded matrix
|
||||
|
||||
if full_matrices:
|
||||
S = np.zeros(X.shape)
|
||||
np.fill_diagonal(S, s)
|
||||
else:
|
||||
S = np.diag(s)
|
||||
|
||||
# Rounding for display
|
||||
return np.round(U, PRECISION), np.round(S, PRECISION), np.round(Vt.T, PRECISION)
|
||||
|
||||
|
||||
def visualize_svd(X,title_X,title_U,title_S,title_V, fig_height=5):
|
||||
|
||||
# Run SVD, as defined above
|
||||
U, S, V = svd(X)
|
||||
|
||||
all_ = np.r_[X.flatten(order='C'),U.flatten(order='C'),
|
||||
S.flatten(order='C'),V.flatten(order='C')]
|
||||
|
||||
# all_max = max(all_.max(),all_.min())
|
||||
# all_min = -max(all_.max(),all_.min())
|
||||
all_max = 6
|
||||
all_min = -6
|
||||
# Visualization
|
||||
fig, axs = plt.subplots(1, 7, figsize=(12, fig_height))
|
||||
|
||||
plt.sca(axs[0])
|
||||
ax = sns.heatmap(X,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title(title_X)
|
||||
|
||||
plt.sca(axs[1])
|
||||
plt.title('@')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[2])
|
||||
ax = sns.heatmap(U,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title(title_U)
|
||||
|
||||
plt.sca(axs[3])
|
||||
plt.title('@')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[4])
|
||||
ax = sns.heatmap(S,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title(title_S)
|
||||
|
||||
plt.sca(axs[5])
|
||||
plt.title('@')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[6])
|
||||
ax = sns.heatmap(V.T,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title(title_V)
|
||||
return X, U, S, V
|
||||
|
||||
# Repeatability
|
||||
np.random.seed(1)
|
||||
|
||||
# Generate random matrix
|
||||
X = np.random.randn(6, 4)
|
||||
|
||||
# manipulate X and reduce rank to 3
|
||||
# X[:,3] = X[:,0] + X[:,1]
|
||||
|
||||
X, U, S, V = visualize_svd(X,'$X$','$U$','$S$','$V^T$', fig_height=3)
|
||||
|
||||
X_2, U_2, S_2, V_2 = visualize_svd(X.T@X,'$X^TX$','$V$','$S^TS$','$V^T$', fig_height=3)
|
||||
|
||||
X_3, U_3, S_3, V_3 = visualize_svd(X@X.T,'$XX^T$','$U$','$SS^T$','$U^T$', fig_height=3)
|
||||
|
||||
#%%
|
||||
|
||||
# Bk4_Ch15_02_B
|
||||
|
||||
#%% U*U.T = I
|
||||
|
||||
all_max = 6
|
||||
all_min = -6
|
||||
|
||||
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
|
||||
|
||||
plt.sca(axs[0])
|
||||
ax = sns.heatmap(U,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$U$')
|
||||
|
||||
|
||||
plt.sca(axs[1])
|
||||
ax = sns.heatmap(U.T,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$U^T$')
|
||||
|
||||
plt.sca(axs[2])
|
||||
ax = sns.heatmap(U@U.T,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$I$')
|
||||
|
||||
#%%
|
||||
|
||||
# Bk4_Ch15_02_C
|
||||
|
||||
#%% V*V.T = I
|
||||
|
||||
all_max = 6
|
||||
all_min = -6
|
||||
|
||||
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
|
||||
|
||||
plt.sca(axs[0])
|
||||
ax = sns.heatmap(V,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$V$')
|
||||
|
||||
plt.sca(axs[1])
|
||||
ax = sns.heatmap(V.T,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$V^T$')
|
||||
|
||||
plt.sca(axs[2])
|
||||
ax = sns.heatmap(V@V.T,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('$I$')
|
||||
|
||||
#%%
|
||||
|
||||
# Bk4_Ch15_02_D
|
||||
|
||||
#%% analysis of singular value matrix
|
||||
|
||||
fig, axs = plt.subplots(1, 4, figsize=(12, 3))
|
||||
|
||||
for j in [0, 1, 2, 3]:
|
||||
X_j = S[j,j]*U[:,j][:, None]@V[:,j][None, :];
|
||||
plt.sca(axs[j])
|
||||
ax = sns.heatmap(X_j,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
title_txt = '$s_'+ str(j+1) + 'u_'+ str(j+1) + 'v_'+ str(j+1) + '^T$'
|
||||
plt.title(title_txt)
|
||||
|
||||
#%% projection
|
||||
|
||||
for j in [0, 1, 2, 3]:
|
||||
|
||||
fig, axs = plt.subplots(1, 7, figsize=(12, 3))
|
||||
|
||||
v_j = V[:,j]
|
||||
v_j = np.matrix(v_j).T
|
||||
s_j = S[j,j]
|
||||
s_j = np.matrix(s_j)
|
||||
u_j = U[:,j]
|
||||
u_j = np.matrix(u_j).T
|
||||
|
||||
plt.sca(axs[0])
|
||||
ax = sns.heatmap(X,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('X')
|
||||
|
||||
plt.sca(axs[1])
|
||||
plt.title('@')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[2])
|
||||
ax = sns.heatmap(v_j,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('v_'+ str(j+1))
|
||||
|
||||
plt.sca(axs[3])
|
||||
plt.title('=')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[4])
|
||||
ax = sns.heatmap(s_j,cmap='RdBu_r',vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('s_'+ str(j+1))
|
||||
|
||||
plt.sca(axs[5])
|
||||
plt.title('@')
|
||||
plt.axis('off')
|
||||
|
||||
plt.sca(axs[6])
|
||||
ax = sns.heatmap(u_j,cmap='RdBu_r', vmax = all_max,vmin = all_min,
|
||||
cbar_kws={"orientation": "horizontal"})
|
||||
ax.set_aspect("equal")
|
||||
plt.title('u_'+ str(j+1))
|
||||
Reference in New Issue
Block a user