Bunch-kaufmann

 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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
import numpy as np
from operator import itemgetter
#--------------------------------------------------------
from math import factorial
def binomial(n, k):
if k < 0 or k > n:
return 0
if k == 0 or k == n:
return 1
return factorial(n) // (factorial(k) * factorial(n-k))
def hilb(n, m=0):
if n < 1 or m < 0:
raise ValueError("Matrix size must be one or greater")
elif n == 1 and (m == 0 or m == 1):
return np.array([[1]])
elif m == 0:
m = n
return 1. / (np.arange(1, n + 1) + np.arange(0, m)[:, np.newaxis])
def invhilb(n):
H = np.empty((n, n))
for i in range(n):
for j in range(n):
H[i, j] = ((-1)**(i + j)) * (i + j + 1) * binomial(n + i, n - j - 1) * \
binomial(n + j, n - i - 1) * binomial(i + j, i) ** 2
return H
#--------------------------------------------------------
def trans_matrix(size, idx1, idx2):
permutation = np.array([[0 for i in xrange(size)] for j in xrange(size)])
for i in xrange(size):
permutation[i][i] = 1
if idx1 != idx2:
permutation[idx1][idx2], permutation[idx2][idx1] = 1
permutation[idx1][idx1], permutation[idx1][idx1] = 0
return np.matrix(permutation)
def bunch_kaufmann(mtx, alpha):
n_th = []
sum = 0
while sum < mtx.shape[0]:
# i-th step of algorithm: we work with n[i]-th mtx_shiffer
mtx_shiffer = mtx[sum: mtx.shape[0], sum: mtx.shape[1]]
# find index with largest diagonal abs
idx = max([(abs(mtx_shiffer[j][j]), j) for j in xrange(mtx_shiffer.shape[0])], key=itemgetter(0))[1]
# make permutation matrix for (0, idx) transposition
permutation = trans_matrix(mtx_shiffer.shape[0], 0, idx)
# conjugate M' with permutation matrix
mtx_shiffer = np.array(permutation*np.matrix(mtx_shiffer)*permutation)
# find index for larger column abs and this abs
[lambda_val, idx] = max([(abs(mtx_shiffer[j][0]), j) for j in xrange(mtx_shiffer.shape[0])], key=itemgetter(0))
n_k = 0
if abs(mtx_shiffer[0][0]) >= alpha*lambda_val:
n_k = 1
# make identity permutation matrix
permutation = trans_matrix(mtx_shiffer.shape[0], 0, 0)
else:
# same as above, for idx-th column of mtx_shiffer (without diagonal element)
[sigma_val, j_idx] = max([(abs(mtx_shiffer[j][idx]), j) for j in xrange(mtx_shiffer.shape[0]) if j != idx], key=itemgetter(0))
if sigma_val*abs(mtx_shiffer[0][0]) >= alpha*lambda_val**2:
n_k = 1
# make identity permutation matrix
permutation = trans_matrix(mtx_shiffer.shape[0], 0, 0)
else:
if abs(mtx_shiffer[idx][idx]) >= alpha*sigma_val:
n_k = 1
# make non-identity permutation matrix (0, idx) for mtx_shiffer
permutation = trans_matrix(mtx_shiffer.shape[0], 0, idx)
# conjugate M' with permutation matrix
mtx_shiffer = permutation*np.matrix(mtx_shiffer)*permutation
else:
n_k = 2
permutation = trans_matrix(mtx_shiffer.shape[0], 2, idx)*trans_matrix(1, j_idx)
# conjugate M' with permutation matrix
mtx_shiffer = permutation*np.matrix(mtx_shiffer)*permutation
sum += n_k
n_th.append(n_k)
return mtx
#mtx = np.loadtxt('matrix.txt')
mtx = hilb(10)
f = open("result.txt", "w")
f.write(str(mtx))
print bunch_kaufmann(mtx, 0.5)