내적 (2)

하우스홀더 방법을 이용한 QR 분해

A = [[1, -1, 4], [1, 4, -2], [1, 4, 2], [1, -1, 0]]

QR 분해할 행렬 A 선언

A1 = deepcopy(A)

a1 = transpose(A)[0]
print(a1)

>>> [1, 1, 1, 1]

deepcopy 함수를 이용해 행렬 A의 내용을 복사한 행렬 A1을 구함

첫번째 열 추출하여 a1이라 함

nm1 = norm(a1)
print(nm1)

>>> 2.0

벡터 a1의 norm 구하기

e1 = [1]
for i in range (0, len(a1)-1):
	e1.append(0)
    
tmp_e1 = []
for i in range (0, len(e1)):
	val = sign(a1[0])*nm1*e1[i]
    tmp_e1.append(val)
    
v1 = v_add(a1, tmp_e1)
print(v1)

>>> [3.0, 1.0, 1.0, 1.0]

e1 벡터 리스트 생성하고 벡터 v1 생성

🔽 v1 생성 공식

v1 = a1 + sign(a1)||a1||e1

# H1 생성
H1 = householder(v1)
print(H1)

행렬 H1 

H1 = I -2(v1v1^T)/(v1^Tv1)

tmp_res1 = matmul(H1, A1)

H1과 A1의 행렬 곱을 구하고 tmp_res1에 저장

matmul 함수 이용

A2 = []
for i in range(1, len(A1)):
	tmp_row = []
    for j in range(1, len(A1[0])):
    	tmp_row.append(tmp_res1[i][j]))
    A2.append(tmp_row)
print(A2)

tmp_res1 행렬의 1행과 1열을 제거한 행렬 A2를 구하는 과정

A2의 원소를 담을 비어있는 리스트를 생성하고 A2에 추가할 행에 대한 반복문 수행.

인덱스는 0이 아닌 1부터 시작해야함. 행렬 tmp_res1의 1행을 제외해야하기 때문

행렬 A2에 추가할 행이 될 비어 있는 리스트 tmp_row 정의하고

A2의 열에 대한 반복문 수행

a2 = transpose(A2)[0]
print(a2)

행렬 A2의 1열을 추출해 벡터 a2 생성

.

.

.

같은 과정 거쳐서 v2랑 v3 생성하고 H2와 H3 생성

tmp_H2 = identity(len(A))
for i in range(1, len(A)):
	for j in range(1, len(A)):
    	tmp_H2[i][j] = H2[i-1][j-1]
H2 = tmp_H2
print(H2)

3행 3열 크기의 행렬 H2를 4행 4열 크기의 행렬로 바꾸는 과정

기본 뼈대가 될 단위 행렬 tmp_H2를 identity 함수를 이용해 만듦

비슷한 과정 거쳐서 2행 2열이었던 H3 행렬을 4행 4열 형태로 변환하기 

 

여기까지 하면 H1, H2, H3를 모두 구할 수 있음

Q = H1*H2*H3
Q = matmul(matmul(H1, H2), H3)
print(Q)

Q = H1H2H3 

행렬 곱 이용해 행렬 Q 구할 수 있음

R = H3 * H2 * H1 * A
R = matmul(matmul(matmul(H3, H2), H1), A)
print(R)

R = H3H2H1A

행렬 R 구하는 식

 

> 함수 생성

여태까지의 과정 일반화하는 함수 생성하여 하우스홀더 방법 이용한 QR 분해 수행

def qr_householder(A):
'''
입력값: 행렬 A
출력값: 행렬 Q, R
'''

n = len(A)
p = len(A[0])

H_list = []

for i in range(0, p):

	if i==0:
    	A1 = deepcopy(A)
        exA = A1
    elif i < p-1:
    	Ai = []
        for j in range(1, len(exA):
        	row = []
            for k in range(1, len(exA[0])):
            Ai.append(row)
        exA = Ai
    elif i==p-1:
    	Ap = []
        for j in range(1, len(HA)):
        	Ap.append(HA[j][i])
        exA = Ap
   
    # 열 추출
    if i<p-1:
    	a = transpose(exA)[0]
    else:
    	a = exA
    nm = norm(a)
    
    # e 생성
    e = [1]
    for j in range(0, len(a)-1):
    	e.append(0)
    
    # v 생성
    tmp_e = []
    for j in range(0, len(e)):
    	val = sign(a[0])*nm*e[j]
        tmp_e.append(val)
    v = v_add(a, tmp_e)
    
    # H 생성
    H = householder(v)
    
    # H*A
    if i==p-1:
    	HA = []
        for j in range(0, len(H[0]):
        	val += H[j][k]*exA[k]
        HA.append(val)
    else: 
    	HA = matmul(H, exA)
    
    H_list.append(H)
    
    if i > 0:
    	tmp_H = identity(len(A))
        for j in range (i, len(A)):
        	for k in range(i, len(A)):
            	tmp_H[j][k] = H_list[-1][j-i][k-i]
        H_list[01] = tmp_H
        
    Q = deepcopy(H_list[0])
    for j in range(0, len(H_list)-1):
    	Q = matmul(Q, H_list[j+1])
    
    R = deepcopy(H_list[-1])
    for j in rnage(1, len(H_list)):
    	R = matmul(R, H_list[-(j+1)])
    R = matmul(R, A)
    
    return Q, R

qr_householder 함수

qr_householder 함수의 입력값은 행렬 A