approximation.py 文件源码

def power_inverse(A, x, tol = 0.0001, N =5000):
    '''
    inverse power method
    faster convergence rate

        tol: tolerance
        N: maximum number of iterations

    '''

    q = np.dot(x.T, np.dot(A, x))/np.dot(x.T, x)

    p = -1
    u0, u1 = 0, 0
    for i in range(A.shape[1]):
        p += 1
        if la.norm(x, np.inf) == x[i]:
            break

    x = x / x[p]

    for k in range(N):
        try:
            y = la.solve(A - q * np.eye(A.shape[0]), x)
        except LinAlgError:
            return q, x
        else:
            u = y[p]

            p = -1

            for i in range(len(y)):
                p += 1
                if la.norm(y, np.inf) == y[i]:
                    break

            err = la.norm(x - (y/y[p]), np.inf)
            x = y/y[p]
            if err < tol:
                u = 1/u + q
                return u, x