Gram-Schmidt orthonormalization & QR decomposition

<Gram-Schmidt orthonormalization & QR decomposition>

    1. Install Packages

    2. Gram-Schmidt orthonormalization

    3. QR decomposition

 

---

Install Packages

# visualization을 위한 helper code입니다.
from urllib.request import urlretrieve
URL = 'https://go.gwu.edu/engcomp4plot'
urlretrieve(URL, 'plot_helper.py')

import sys
sys.path.append('../scripts/')

# 다음 세 custom function (1)plot_vector, (2)plot_linear_transformation, (3) plot_linear_transformations
# 을 사용할 것입니다.
from plot_helper import *

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import scipy as sp
import scipy.linalg
import sympy as sy

sy.init_printing() 
np.set_printoptions(precision=3)
np.set_printoptions(suppress=True)

---

Gram-Schmidt orthonormalization

 

Gram-Schmidt orthonormalization:

주어진 basis를 orthonormal basis로 만들기 위한 method

 

 

{x1, x2, x3}를 R3의 basis라고 하고, x1= (3, 6, 2), x2 = (1, 2, 4), x3 = (2, -2, 1)라고 할 시,

x1, x2, x3은 서로 orthogonal이 아니다. (x1·x2 != 0)

 

→ Gram-Schmidt orthonormalization을 이용해 이를 orthonormal basis u1, u2, u3로 변환

 

s = np.linspace(-1, 1, 10)
t = np.linspace(-1, 1, 10)
S, T = np.meshgrid(s, t)

vec = np.array([[[0,0,0,3, 6, 2]],
             [[0,0,0,1, 2, 4]],
             [[0,0,0,2, -2, 1]]])

# x1과 x_2의 linear combination, 즉 span
X = vec[0,:,3] * S + vec[1,:,3] * T
Y = vec[0,:,4] * S + vec[1,:,4] * T
Z = vec[0,:,5] * S + vec[1,:,5] * T

fig = plt.figure(figsize = (7, 7))
ax = fig.add_subplot(projection='3d')

# x1, x2가 span하는 공간을 격자로 나타냅니다. 
ax.plot_wireframe(X, Y, Z, linewidth = 1.5, alpha = .3)

############################# x1 and x2 ##############################
colors = ['r','b','g']
s = ['$x_1$', '$x_2$', '$x_3$']
for i in range(vec.shape[0]):
    X,Y,Z,U,V,W = zip(*vec[i,:,:])
    ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False,
              color = colors[i], alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
              linestyles = 'solid',linewidths = 3)
    ax.text(vec[i,:,3][0], vec[i,:,4][0], vec[i,:,5][0], s = s[i], size = 15)

ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')
plt.show()

파란색 격자 = x1, x2가 span하는 2-dimensional subspace

 

x1 vector를 고정 (v1 = x1)

 

span {v1}을 W1이라고 할 시, v2를 구하기 위해 x2를 W1에 orthogonal projection한 vector

x1 = np.array([3, 6, 2])
x2 = np.array([1, 2, 4])
x3 = np.array([2, -2, 1])

v1 = x1 
u1 = v1 / (np.sqrt(v1.dot(v1)))
x2_hat = ((x2.dot(u1)) / (u1.dot(u1))) * u1
v2 = x2 - x2_hat
u2 = v2 / (np.sqrt(v2.dot(v2)))

print(u1 @ u2)
print(x2_hat @ u2)

s = np.linspace(-1, 1, 10)
t = np.linspace(-1, 1, 10)
S, T = np.meshgrid(s, t)

# x1과 x_2의 linear combination, 즉 span
X = x1[0] * S + x2[0] * T
Y = x1[1] * S + x2[1] * T
Z = x1[2] * S + x2[2] * T

fig = plt.figure(figsize = (10, 10))
ax = fig.add_subplot(projection='3d')

# x1, x2가 span하는 공간을 격자로 나타냅니다. 
ax.plot_wireframe(X, Y, Z, linewidth = 1.5, alpha = .3)

############################# x1, x2, x3, x2_hat, v2 ##############################

vec = np.array([[0, 0, 0, x1[0], x1[1], x1[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'red', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x2[0], x2[1], x2[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'blue', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x3[0], x3[1], x3[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'green', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x2_hat[0],x2_hat[1], x2_hat[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'blue', alpha = .6,arrow_length_ratio = .12, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, v2[0], v2[1], v2[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'purple', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)



ax.text(x1[0], x1[1], x1[2], '$\mathbf{x}_1 = \mathbf{v}_1 $', size = 15)
ax.text(x2[0], x2[1], x2[2], '$\mathbf{x}_2$', size = 15)
ax.text(x3[0], x3[1], x3[2], '$\mathbf{x}_3$', size = 15)
ax.text(x2_hat[0], x2_hat[1], x2_hat[2], '$\hat{\mathbf{x_2}}$', size = 15)
ax.text(v2[0], v2[1], v2[2], '$\mathbf{v}_2$', size = 15)

ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')

################################# Dashed Line ##################################

point1 = [x2_hat[0], x2_hat[1], x2_hat[2]]
point2 = [x2[0], x2[1], x2[2]]
line1 = np.array([point1, point2])
ax.plot(line1[:,0], line1[:,1], line1[:, 2], c = 'b', lw = 3.5,alpha =0.5, ls = '--')

point1 = [v2[0], v2[1], v2[2]]
point2 = [x2[0], x2[1], x2[2]]
line1 = np.array([point1, point2])
ax.plot(line1[:,0], line1[:,1], line1[:, 2], c = 'b', lw = 3.5,alpha =0.5, ls = '--')
plt.show()

 

 

v1, v2를 구했으니, 세 번째 basis인 v3를 구하려면,

span {v1, v2}를 W2라고 하고, x3를 W2에 orthogonal projection 한 vector

 

x3_hat = (x3.dot(u1)) / (u1.dot(u1)) * u1 + (x3.dot(u2)) / (u2.dot(u2)) * u2
v3 = x3 - x3_hat
u3 = v3 / (np.sqrt(v3.dot(v3)))

print(u1 @ u3)
print(u2 @ u3)
print(x3_hat @ u3)

s = np.linspace(-1, 1, 10)
t = np.linspace(-1, 1, 10)
S, T = np.meshgrid(s, t)

# x1과 x_2의 linear combination, 즉 span
X = x1[0] * S + x2[0] * T
Y = x1[1] * S + x2[1] * T
Z = x1[2] * S + x2[2] * T

fig = plt.figure(figsize = (9, 9))
ax = fig.add_subplot(projection='3d')
# x1, x2가 span하는 공간을 격자로 나타냅니다. 
ax.plot_wireframe(X, Y, Z, linewidth = 1.5, alpha = .3)

############################# x1, x2, x2_hat, v2, x3, x3_hat, v3 ##############################

vec = np.array([[0, 0, 0, x1[0], x1[1], x1[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'red', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x2[0], x2[1], x2[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'red', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x2_hat[0],x2_hat[1], x2_hat[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'blue', alpha = .6,arrow_length_ratio = .12, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, v2[0], v2[1], v2[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'purple', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)

vec = np.array([[0, 0, 0, x3[0], x3[1], x3[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'red', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)


vec = np.array([[0, 0, 0, x3_hat[0], x3_hat[1], x3_hat[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'black', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)


vec = np.array([[0, 0, 0, v3[0], v3[1], v3[2]]])
X, Y, Z, U, V, W = zip(*vec)
ax.quiver(X, Y, Z, U, V, W, length=1, normalize=False, color = 'purple', alpha = .6,arrow_length_ratio = .08, pivot = 'tail',
          linestyles = 'solid',linewidths = 3)


ax.text(x1[0], x1[1], x1[2], '$\mathbf{x}_1 = \mathbf{v}_1 $', size = 15)
ax.text(x2[0], x2[1], x2[2], '$\mathbf{x}_2$', size = 15)
ax.text(x2_hat[0], x2_hat[1], x2_hat[2], '$\hat{\mathbf{x}}_2$', size = 15)
ax.text(v2[0], v2[1], v2[2], '$\mathbf{v}_2$', size = 15)
ax.text(x3[0], x3[1], x3[2], '$\mathbf{x}_3$', size = 15)
ax.text(x3_hat[0], x3_hat[1], x3_hat[2], '$\hat{\mathbf{x}}_3$', size = 15)
ax.text(v3[0], v3[1], v3[2], '$\mathbf{v}_3$', size = 15)

ax.set_xlabel('x-axis')
ax.set_ylabel('y-axis')
ax.set_zlabel('z-axis')

################################# Dashed Line ##################################


point1 = [x2_hat[0], x2_hat[1], x2_hat[2]]
point2 = [x2[0], x2[1], x2[2]]
line1 = np.array([point1, point2])
ax.plot(line1[:,0], line1[:,1], line1[:, 2], c = 'b', lw = 3.5,alpha =0.5, ls = '--')

point1 = [v2[0], v2[1], v2[2]]
point2 = [x2[0], x2[1], x2[2]]
line1 = np.array([point1, point2])
ax.plot(line1[:,0], line1[:,1], line1[:, 2], c = 'b', lw = 3.5,alpha =0.5, ls = '--')


point1 = [x3_hat[0], x3_hat[1], x3_hat[2]]
point2 = [x3[0], x3[1], x3[2]]
line1 = np.array([point1, point2])
ax.plot(line1[:,0], line1[:,1], line1[:, 2], c = 'b', lw = 3.5,alpha =0.5, ls = '--')


################################ Axes ######################################
ax.set_xlim3d(-5, 5)
ax.set_ylim3d(-5, 5)
ax.set_zlim3d(-5, 5)
plt.show()

orthonormal basis {v1, v2, v3}

 

• u1은 x1의 Linear Combination (u1 = x1)
• u2는 x1, x2의 Linear Combination
• u3는 x1, x2, x3의 Linear Combination

→ Span의 subspace가 같다.

→ 세 vector의 Linear Combination으로 만들 수 있는 vector가 같다.

→ span {x1, x2, x3} = span {u1, u2, u3}

 

 

 

따라서, 세 vector u1, u2, u3를 각 column으로 하는 Matrix U는 orthogonal Matrix가 된다.

U = np.vstack((u1,u2,u3)).T
U.T@U

U@U.T  # UU^T는 identity matrix가 됩니다.

---

QR decomposition

 

Gram-Schmidt를 하면 자연스럽게 Matrix A의 QR decomposition을 수행할 수 있다.

A = np.array([[3,1,2],[6,2,-2],[2,4,1]])

Q = np.vstack((u1,u2,u3)).T

r1 = np.array([np.sqrt(v1.dot(v1)), 0, 0])
r2 = np.array([x2.dot(u1) / u1.dot(u1), np.sqrt(v2.dot(v2)), 0])
r3 = np.array([x3.dot(u1) / u1.dot(u1), x3.dot(u2) / u2.dot(u2), np.sqrt(v3.dot(v3))])
R = np.vstack((r1, r2, r3)).T

print(Q)
print(R)
print('')
print(A)
print(Q@R)

 

 

Numpy의 linalg.qr 함수 이용

Q,R = np.linalg.qr(A)

print(Q)
print(R)