Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Power method
class LinearMap:
def __init__(self, matrix):
self.domain_dimension = len(matrix[0])
self.range_dimension = len(matrix)
self.matrix = matrix
# def shifted(self, c):
# return LinearMap(tuple(tuple(r - c if i == j else r
# for j, r in enumerate(row))
# for i, row in enumerate(self.matrix)))
def __call__(self, v):
return tuple(sum(a * b for a, b in zip(row, v)) for row in self.matrix)
def __sub__(self, other):
return LinearMap(tuple(tuple(x - y for x, y in zip(x_row, y_row))
for x_row, y_row in zip(self, other)))
def __iter__(self):
yield from self.matrix
def __len__(self):
return self.range_dimension
def Diag(*xs):
return LinearMap(tuple(tuple(xs[i] if i == j else 0 for i in range(len(xs)))
for j in range(len(xs))))
def powermethod(f, x=None, precision=1e-3):
x = x if x is not None else tuple(1 for _ in range(f.domain_dimension))
m0, m1 = float('-inf'), float('inf')
while abs(m1 - m0) > precision:
m0, m1 = m1, max(x, key=abs)
x = f(tuple(a / m1 for a in x))
return m1
def shifted_powermethod(f, shift, *args, **kwargs):
shifted = f - Diag(*(shift for _ in range(len(f))))
return shift + powermethod(shifted, *args, **kwargs)
#L = LinearMap(((91/10, -2/10), (-5/10, 0)))
#print(u := powermethod(L))
#print(101 + 1/u)
#print(l1 := powermethod(L))
#print(shifted_powermethod(L, l1))
#L = LinearMap(((367/16, 1, -1/16), (1, 0, 0), (-1/16, 0, -1/16)))
#print(20+1/powermethod(L))
L = LinearMap(((1, 1), (1, 0)))
print(u := powermethod(L))
print(shifted_powermethod(L, u))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment