Power method

powermethod.py

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))