inverse of a triangular matrix (2 heuristics)

 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
import numpy as np
from numpy import identity as I
from numpy import tril, triu, dot
from collections import Counter
from operator import itemgetter
def triangular_inversion(triang_arg):
"""Counts inversion of a triangular matrix (lower or upper).
NOT TESTED PROPERLY
function trying to predict form of inversion and count it without heavy computations
SEE gauss atomic triangular matrix and it's block analog
NOT TESTED PROPERLY
function also trying to predict power of nilpotence (atomic_number)
of nilpotent matrix inside of an algorithm
Args:
triang_arg (np.array): triangular matrix for inversion
Returns:
np.array: inverse of triangular matrix
Raises:
Exception: An error occurred while passing non-square matrix
Exception: An error occurred while passing non-triangular matrix
Exception: An error occurred while passing singular matrix
"""
if len(triang_arg.shape) != 2 or triang_arg.shape[0] != triang_arg.shape[1]:
raise Exception('Matrix is non-square')
if not np.array_equal(triang_arg, tril(triang_arg)) and not np.array_equal(triang_arg, triu(triang_arg)):
raise Exception('Matrix is not triangular')
if not len(triang_arg.diagonal().nonzero()[0]):
raise Exception('Matrix is singular')
triang = triang_arg.copy()
n = triang.shape[0]
unitriang_maker = I(n) / triang.diagonal()
unitriang = dot(unitriang_maker, triang)
nilpotent = unitriang - I(n)
# check for possibility of simple inversion
z = zip(nilpotent.nonzero()[0], nilpotent.nonzero()[1])
if z[0][0] > z[0][1]:
# lower triangular case
i_min = min(z, key=itemgetter(0))[0]
j_max = max(z, key=itemgetter(1))[1]
if i_min > j_max:
return dot(I(n) - nilpotent, unitriang_maker)
else:
# upper triangular case
i_max = max(z, key=itemgetter(0))[0]
j_min = min(z, key=itemgetter(1))[1]
if j_min > i_max:
return dot(I(n) - nilpotent, unitriang_maker)
# nilpotence power prediction
atomic_number = len(Counter([column for (row, column) in z]))
unitriang_inverse = I(n)
for i in xrange(atomic_number):
unitriang_inverse = I(n) - dot(nilpotent, unitriang_inverse)
return dot(unitriang_inverse, unitriang_maker)
# test
# make random matrix sample, take it low triangle and up triangle
n = 7
sample = mtx(np.random.random((n, n)))
print sample
L = np.tril(sample)
U = np.triu(sample)
# find inversions
L_I = triangular_inversion(L)
U_I = triangular_inversion(U)
# check result
print 'Check: L * L_I'
print L * L_I
print 'Check: L_I * L'
print L_I * L
print 'Check: U * U_I'
print U * U_I
print 'Check: U_I * U'
print U_I * U