Editing Arnoldi iteration (section)

==The Arnoldi iteration==

The Arnoldi iteration uses the [[Gram-Schmidt#Numerical stability|modified Gram–Schmidt process]] to produce a sequence of orthonormal vectors, ''q''<sub>1</sub>, ''q''<sub>2</sub>, ''q''<sub>3</sub>, ..., called the ''Arnoldi vectors'', such that for every ''n'', the vectors ''q''<sub>1</sub>, ..., ''q''<sub>''n''</sub> span the Krylov subspace <math>\mathcal{K}_n</math>. Explicitly, the algorithm is as follows:

 '''Start''' with an arbitrary vector ''q''<sub>1</sub> with norm 1.
 '''Repeat for''' ''k'' = 2, 3, ...
   ''q''<sub>''k''</sub> := ''A'' ''q''<sub>''k''−1</sub>
   '''for''' ''j'' '''from''' 1 '''to''' ''k'' − 1
     ''h''<sub>''j'',''k''−1</sub> :=  ''q''<sub>''j''</sub><sup>*</sup> ''q''<sub>''k''</sub>
     ''q''<sub>''k''</sub> := ''q''<sub>''k''</sub> − h<sub>''j'',''k''−1</sub> ''q''<sub>''j''</sub>
   ''h''<sub>''k'',''k''−1</sub> := {{norm|''q''<sub>''k''</sub>}}
   ''q''<sub>''k''</sub> := ''q''<sub>''k''</sub> / ''h''<sub>''k'',''k''−1</sub>

The ''j''-loop projects out the component of <math>q_k</math> in the directions of <math>q_1,\dots,q_{k-1}</math>.  This ensures the orthogonality of all the generated vectors.

The algorithm breaks down when ''q''<sub>''k''</sub> is the zero vector. This happens when the [[Minimal polynomial (linear algebra)|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 4''mk'' floating point operations.

In the programming language [[Python (programming language)|Python]] with support of the [[NumPy]] library:
<syntaxhighlight lang="numpy">
import numpy as np

def arnoldi_iteration(A, b, n: int):
    """Compute a basis of the (n + 1)-Krylov subspace of the matrix A.

    This is the space spanned by the vectors {b, Ab, ..., A^n b}.

    Parameters
    ----------
    A : array_like
        An m × m array.
    b : array_like
        Initial vector (length m).
    n : int
        One less than the dimension of the Krylov subspace, or equivalently the *degree* of the Krylov space. Must be >= 1.
    
    Returns
    -------
    Q : numpy.array
        An m x (n + 1) array, where the columns are an orthonormal basis of the Krylov subspace.
    h : numpy.array
        An (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].conj(), 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
</syntaxhighlight>