(A, maxiter, symmetric=None, initial_guess=None)
return (Vects, Eigs, H, V, breakdown_flag)
Apprixmate eigenvalues.

    Used by approximate_spectral_radius and condest.

    Returns [W, E, H, V, breakdown_flag], where W and E are the eigenvectors
    and eigenvalues of the Hessenberg matrix H, respectively, and V is the
    Krylov space.  breakdown_flag denotes whether Lanczos/Arnoldi suffered
    breakdown.  E is therefore the approximate eigenvalues of A.

    To obtain approximate eigenvectors of A, compute V*W.
Apprixmate eigenvalues.
Apprixmate eigenvalues .
def _approximate_eigenvalues(A, maxiter, symmetric=None, initial_guess=None):
    """Apprixmate eigenvalues.

    Used by approximate_spectral_radius and condest.

    Returns [W, E, H, V, breakdown_flag], where W and E are the eigenvectors
    and eigenvalues of the Hessenberg matrix H, respectively, and V is the
    Krylov space.  breakdown_flag denotes whether Lanczos/Arnoldi suffered
    breakdown.  E is therefore the approximate eigenvalues of A.

    To obtain approximate eigenvectors of A, compute V*W.

    """
    A = aslinearoperator(A)  # A could be dense or sparse, or something weird

    # Choose tolerance for deciding if break-down has occurred
    breakdown = set_tol(A.dtype)
    breakdown_flag = False

    if A.shape[0] != A.shape[1]:
        raise ValueError('expected square matrix')

    maxiter = min(A.shape[0], maxiter)

    if initial_guess is None:
        v0 = np.random.rand(A.shape[1], 1)
        if A.dtype == complex:
            v0 = v0 + 1.0j * np.random.rand(A.shape[1], 1)
    else:
        v0 = initial_guess

    v0 /= norm(v0)

    # Important to type H based on v0, so that a real nonsymmetric matrix, can
    # have an imaginary initial guess for its Arnoldi Krylov space
    H = np.zeros((maxiter+1, maxiter),
                 dtype=np.find_common_type([v0.dtype, A.dtype], []))

    V = [v0]

    beta = 0.0
    for j in range(maxiter):
        w = A * V[-1]

        if symmetric:
            if j >= 1:
                H[j-1, j] = beta
                w -= beta * V[-2]

            alpha = np.dot(np.conjugate(w.ravel()), V[-1].ravel())
            H[j, j] = alpha
            w -= alpha * V[-1]  # axpy(V[-1],w,-alpha)

            beta = norm(w)
            H[j+1, j] = beta

            if (H[j+1, j] < breakdown):
                breakdown_flag = True
                break

            w /= beta

            V.append(w)
            V = V[-2:]  # retain only last two vectors

        else:
            # orthogonalize against Vs
            for i, v in enumerate(V):
                H[i, j] = np.dot(np.conjugate(v.ravel()), w.ravel())
                w = w - H[i, j]*v

            H[j+1, j] = norm(w)

            if (H[j+1, j] < breakdown):
                breakdown_flag = True
                if H[j+1, j] != 0.0:
                    w = w/H[j+1, j]
                V.append(w)
                break

            w = w/H[j+1, j]
            V.append(w)

            # if upper 2x2 block of Hessenberg matrix H is almost symmetric,
            # and the user has not explicitly specified symmetric=False,
            # then switch to symmetric Lanczos algorithm
            # if symmetric is not False and j == 1:
            #    if abs(H[1,0] - H[0,1]) < 1e-12:
            #        #print("using symmetric mode")
            #        symmetric = True
            #        V = V[1:]
            #        H[1,0] = H[0,1]
            #        beta = H[2,1]

    # print("Approximated spectral radius in %d iterations" % (j + 1))

    Eigs, Vects = eig(H[:j+1, :j+1], left=False, right=True)

    return (Vects, Eigs, H, V, breakdown_flag)
def _approximate_eigenvalues ( A , maxiter , symmetric = None , initial_guess = None ) : A = aslinearoperator ( A ) # A could be dense or sparse, or something weird # Choose tolerance for deciding if break-down has occurred breakdown = set_tol ( A . dtype ) breakdown_flag = False if A . shape [ 0 ] != A . shape [ 1 ] : raise ValueError ( 'expected square matrix' ) maxiter = min ( A . shape [ 0 ] , maxiter ) if initial_guess is None : v0 = np . random . rand ( A . shape [ 1 ] , 1 ) if A . dtype == complex : v0 = v0 + 1.0j * np . random . rand ( A . shape [ 1 ] , 1 ) else : v0 = initial_guess v0 /= norm ( v0 ) # Important to type H based on v0, so that a real nonsymmetric matrix, can # have an imaginary initial guess for its Arnoldi Krylov space H = np . zeros ( ( maxiter + 1 , maxiter ) , dtype = np . find_common_type ( [ v0 . dtype , A . dtype ] , [ ] ) ) V = [ v0 ] beta = 0.0 for j in range ( maxiter ) : w = A * V [ - 1 ] if symmetric : if j >= 1 : H [ j - 1 , j ] = beta w -= beta * V [ - 2 ] alpha = np . dot ( np . conjugate ( w . ravel ( ) ) , V [ - 1 ] . ravel ( ) ) H [ j , j ] = alpha w -= alpha * V [ - 1 ] # axpy(V[-1],w,-alpha) beta = norm ( w ) H [ j + 1 , j ] = beta if ( H [ j + 1 , j ] < breakdown ) : breakdown_flag = True break w /= beta V . append ( w ) V = V [ - 2 : ] # retain only last two vectors else : # orthogonalize against Vs for i , v in enumerate ( V ) : H [ i , j ] = np . dot ( np . conjugate ( v . ravel ( ) ) , w . ravel ( ) ) w = w - H [ i , j ] * v H [ j + 1 , j ] = norm ( w ) if ( H [ j + 1 , j ] < breakdown ) : breakdown_flag = True if H [ j + 1 , j ] != 0.0 : w = w / H [ j + 1 , j ] V . append ( w ) break w = w / H [ j + 1 , j ] V . append ( w ) # if upper 2x2 block of Hessenberg matrix H is almost symmetric, # and the user has not explicitly specified symmetric=False, # then switch to symmetric Lanczos algorithm # if symmetric is not False and j == 1: #    if abs(H[1,0] - H[0,1]) < 1e-12: #        #print("using symmetric mode") #        symmetric = True #        V = V[1:] #        H[1,0] = H[0,1] #        beta = H[2,1] # print("Approximated spectral radius in %d iterations" % (j + 1)) Eigs , Vects = eig ( H [ : j + 1 , : j + 1 ] , left = False , right = True ) return ( Vects , Eigs , H , V , breakdown_flag )