Search code examples
pythonnumpylinear-algebraeigenvector

Eigen vectors in python giving seemingly random element-wise signs


I'm running the following code:

import numpy as np
import matplotlib
matplotlib.use("TkAgg")
import matplotlib.pyplot as plt

N = 100
t = 1

a1 = np.full((N-1,), -t)
a2 = np.full((N,), 2*t)
Hamiltonian = np.diag(a1, -1) +  np.diag(a2) + np.diag(a1, 1)

eval, evec = np.linalg.eig(Hamiltonian)
idx = eval.argsort()[::-1]
eval, evec = eval[idx], evec[:,idx]

wave2 = evec[2] / np.sum(abs(evec[2]))
prob2 = evec[2]**2 / np.sum(evec[2]**2)

_ = plt.plot(wave2)
_ = plt.plot(prob2)
plt.show()

And the plot that comes out is this: enter image description here

But I'd expect the blue line to be a sinoid as well. This has got me confused and I can't find what's causing the sudden sign changes. Plotting the function absolutely shows that the values associated with each x are fine, but the signs are screwed up.

Any ideas on what might cause this or how to solve it?


Solution

  • Here's a modified version of your script that does what you expected. The changes are:

    • Corrected the indexing for the eigenvectors; they are the columns of evec.
    • Use np.linalg.eigh instead of np.linalg.eig. This isn't strictly necessary, but you might as well use the more efficient code.
    • Don't reverse the order of the sorted eigenvalues. I keep the eigenvalues sorted from lowest to highest. Because eigh returns the eigenvalues in ascending order, I just commented out the code that sorts the eigenvalues.

    (Only the first change is a required correction.)


    import numpy as np
    import matplotlib
    import matplotlib.pyplot as plt
    
    N = 100
    t = 1
    
    a1 = np.full((N-1,), -t)
    a2 = np.full((N,), 2*t)
    Hamiltonian = np.diag(a1, -1) +  np.diag(a2) + np.diag(a1, 1)
    
    eval, evec = np.linalg.eigh(Hamiltonian)
    #idx = eval.argsort()[::-1]
    #eval, evec = eval[idx], evec[:,idx]
    
    k = 2
    wave2 = evec[:, k] / np.sum(abs(evec[:, k]))
    prob2 = evec[:, k]**2 / np.sum(evec[:, k]**2)
    
    _ = plt.plot(wave2)
    _ = plt.plot(prob2)
    plt.show()
    

    The plot:

    plot