Try demmel's implementation.

src/afgl/util/lanczos.py

@@ -3,6 +3,7 @@ import numpy.linalg as LA
 
 """
 Classic Lanczos method (without re-orthogonalization)
+Using Demmel's book version.
 
 Arguments
 L : Real valued NxN symmetric matrix
@@ -23,22 +24,20 @@ beta : ndarray
 def lanczos(L, s, M):
     N = len(s)
     alp = np.zeros(M)
-    beta = np.zeros(M - 1)
+    beta = np.zeros(M)
-    V = np.zeros((N, M))
+    V = np.zeros((N, M + 1))
-    V[:, 0] = s / LA.norm(s)
+    V[:, 1] = s / LA.norm(s)
 
-    for j in range(M):
+    for j in range(1, M):
         w = L @ V[:, j]
         alp[j] = np.dot(V[:, j], w)
 
-        v_tilde = w - V[:, j] * alp[j]
+        w = w - V[:, j] * alp[j] - V[:, j - 1] * beta[j - 1]
-        if j > 0:
-            v_tilde = v_tilde - V[:, j - 1] * beta[j - 1]
 
-        if j < M - 1:
+        beta[j] = LA.norm(w)
-            beta[j] = LA.norm(v_tilde)
         if beta[j] == 0:
+            print("Breakdown")
             break
-            V[:, j + 1] = v_tilde / beta[j]
+        V[:, j + 1] = w / beta[j]
 
-    return [V, alp, beta]
+    return [V[:, 1:], alp, beta[1:]]