generated vectors. The algorithm breaks down when qk is the zero vector. This happens when the minimal polynomial of A is of degree k. In most applications of the Arnoldi iteration, including the eigenvalue algorithm below and GMRES, the algorithm has converged at this point. Every step of the k-loop takes one matrix-vector product and approximately 4mk floating point operations. In the programming language Python with support of the NumPy library: import numpy as np def arnoldi_iteration(A, b, n: int): """Computes a basis of the (n + 1)-Krylov subspace of A: the space spanned by {b, Ab, ..., A^n b}. Arguments A: m × m array b: initial vector (length m) n: dimension of Krylov subspace, must be >= 1 Returns Q: m x (n + 1) array, the columns are an orthonormal basis of the Krylov subspace. h: (n + 1) x n array, A on basis Q. It is upper Hessenberg. """ eps = 1e-12 h = np.zeros((n+1,n)) Q = np.zeros((A.shape[0],n+1)) # Normalize the input vector Q[:,0] = b / np.linalg.norm(b,2) # Use it as the first Krylov vector for k in range(1,n+1): v = np.dot(A, Q[:,k-1]) # Generate a new candidate vector for j in range(k): # Subtract the projections on previous vectors h[j,k-1] = np.dot(Q[:,j].T, v) v = v - h[j,k-1] * Q[:,j] h[k,k-1] = np.linalg.norm(v,2) if h[k,k-1] > eps: # Add the produced vector to the list, unless Q[:,k] = v/h[k,k-1] else: # If that happens, stop iterating. return Q, h return Q, h Properties of the Arnoldi iteration Let Qn denote the m-by-n matrix formed by the first n Arnoldi vectors q1, q2, ..., qn, and let Hn be the (upper Hessenberg) matrix formed by the numbers hj,k computed by the algorithm: The orthogonalization method has to be specifically chosen such that the lower Arnoldi/Krylov components are removed from higher Krylov vectors. As can be expressed in terms of q1, ..., qi+1 by construction, they are orthogonal to qi+2, ..., qn, We then have The matrix Hn can be viewed as A in the subspace with the Arnoldi vectors as an orthogonal basis; A is orthogonally projected onto . The matrix Hn can be characterized by the following optimality condition. The characteristic polynomial of Hn minimizes ||p(A)q1||2 among all monic polynomials of degree n. This optimality problem has a unique solution if and only if the Arnoldi iteration does not break down. The relation