forked from tomer279/NumericalAnalysisMethods
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Householder method.py
66 lines (57 loc) · 2.21 KB
/
Householder method.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
import numpy as np
np.set_printoptions(suppress=True)
# Householder method for matrices to obtain a matrix that is similar to that matrix and "simple".
# if the matrix is symmetric, than the method returns a tridiagonal symmetric matrix that is similar to it.
# if not, the method returns a hessenberg matrix which is similar to the matrix as well.
# the code is based on algorithm 9.5 from the book "Numerical Analysis" by Burden, 10th Ed.
# PARAMETERS:
# A -> a matrix.
def householder(A: np.matrix):
n = A.shape[0]
for k in range(0, n - 2):
q = 0
for j in range(k + 1, n):
q += A[j, k] ** 2
if A[k + 1, k] == 0:
alph = -(q ** 0.5)
else:
alph = - (q ** 0.5) * A[k + 1, k] / np.abs(A[k + 1, k])
rsq = alph ** 2 - alph * A[k + 1, k]
v = np.zeros(n)
v[k + 1] = A[k + 1, k] - alph
for j in range(k + 2, n):
v[j] = A[j, k]
if np.array_equal(A, A.T):
u = (1 / rsq) * np.ravel(A.dot(v))
z = u - (1/(2*rsq)) * np.ravel(v.dot(u)) * v
A = A - np.outer(v, z) - np.outer(z, v)
else:
u = (1 / rsq) * np.ravel(A.dot(v))
y = np.zeros(n)
for j in range(0, n):
sig = 0
for i in range(k + 1, n):
sig += A[i, j] * v[i]
y[j] = sig / rsq
prod = np.ravel(u.dot(v))
z = np.zeros(n)
for j in range(0, n):
z[j] = u[j] - (prod / rsq) * v[j]
B = A.copy()
for l in range(k + 1, n):
for j in range(0, k + 1):
A[j, l] = B[j, l] - z[j] * v[l]
A[l, j] = B[l, j] - y[j] * v[l]
for j in range(k + 1, n):
A[j, l] = B[j, l] - z[j] * v[l] - y[l] * v[j]
return A
# input parameters
A_mat = np.matrix([[2, -1, 3],
[2, 0, 1],
[-2, 1, 4]], dtype=float)
B_mat = np.matrix([[4,1,-2,2],
[1,2,0,1],
[-2,0,3,-2],
[2,1,-2,-1]], dtype=float)
print(householder(A_mat))
print(householder(B_mat))