Unter (numerischen) Verfahren zur Lösung des Eigenwertproblems \begin{eqnarray} Ax=\lambda x\quad \text{order}\quad \left(A-\lambda I \right)x=0 \end{eqnarray} mit A ∈ ℝ^n×n versteht man Verfahren, welche die Lösung (λ, x), λ ∈ ℂ, x ∈ ℂⁿ als Grenzwert einer unendlichen Folge von Näherungslösungen berechnen.

Zur numerischen Lösung des Eigenwertproblems wurden zahlreiche Verfahren entwickelt, welche sich in zwei Verfahrensklassen unterteilen lassen. Welches Verfahren das geeignete zur Lösung des gegebenen Problems ist, hängt stark von der Problemstellung und der Größe der Matrix A ab.

Transformationsmethoden sind Methoden, bei denen die Matrix A selbst in jedem Iterationsschritt verändert wird. Ziel der Iteration ist es, A in eine Matrix zu transformieren, von der die Eigenwerte abgelesen werden können. Sie berechnen so alle Eigenwerte einer Matrix und sind hauptsächlich für Matrizen kleiner bis mittlerer Größe geeignet. Die zweite Verfahrensklasse sind Iterationsverfahren, bei denen A nicht verändert wird und man lediglich Matrix-Vektor-Produkte Ax benötigt. Hier wird häufig versucht, einzelne Eigenwert-Eigenvektor-Paare zu approximieren. Diese Verfahren sind besonders gut für große sparse Matrizen geeignet.

Da ein Eigenvektor x ungleich Null ist, muß, damit das homogene Gleichungssystem (A − λI)x = 0 eine nichttriviale Lösung x hat, die Matrix A − λI singular sein. Mit anderen Worten es gilt \begin{eqnarray} \det \left(A-\lambda I \right)=0. \end{eqnarray} Folglich muß der Eigenwert λ eine Nullstelle des charakteristischen Polynoms \begin{eqnarray}\begin{array}{lll}{p}_{n}(\lambda) & = & \det (A-\lambda I)\\ & = & {(-\lambda)}^{n}+{a}_{n-1}{\lambda}^{n-1}+{a}_{n-2}{\lambda}^{n-2}\\ & & +\cdots +{a}_{1}\lambda +{a}_{0}\end{array}\end{eqnarray} sein. Das charakteristische Polynom einer (n × n)- Matrix A hat stets den Grad n und damit genau n Nullstellen, wenn man diese entsprechend ihrer Vielfachheit zählt. Jede Matrix A hat also genau n Eigenwerte λ₁,…,λ_n; die Vielfachheit eines Eigenwertes als Nullstelle des charakteristischen Polynoms wird algebraische Vielfachheit des Eigenwerts genannt. Eigenvektoren zu einem Eigenwert sind nicht eindeutig bestimmt, sie bilden einen invarianten Unterraum von ℂⁿ, dessen Dimension die geometrische Vielfachheit des Eigenwerts ist. Es gilt: Die geometrische Vielfachheit eines Eigenwerts ist kleiner oder gleich der algebraische Vielfachheit. Zu einem ℓ-fachen Eigenwert β (d. h., einer ℓ-fachen Wurzel p_n(λ)) brauchen also nicht notwendigerweise ℓ linear unabhängige Eigenvektoren x gehören. In einem solchen Falle ist man typischerweise nicht an den einzelnen Eigenvektoren interessiert, sondern daran, einen invarianten Unterraum S ∈ ℂ^n×ℓ so zu bestimmen, daß \begin{eqnarray} AS=S\Delta \end{eqnarray} gilt, wobei Λ eine (ℓ × ℓ)-Matrix ist, welche den ℓ-fachen Eigenwert β besitzt.

Die Charakterisierung der Eigenwerte von A als Wurzeln des charakteristischen Polynoms legt es nahe, die Koeffizienten a_j von p_n durch Berechnen der Determinante det(A − λI) zu bestimmen und danach die Eigenwerte λ₁,…,λ_n als Nullstellen von p_n zu berechnen. Die Eigenvektoren ergeben sich dann als nichttriviale Lösungen der Gleichungen \begin{eqnarray} \left(A-{{\lambda}_{k}}I \right){{x}_{k}}=0. \end{eqnarray} Dieser Weg führt jedoch zu einem extrem instabilen Algorithmus, da jede numerische Methode zur Berechnung der Koeffizienten a_j des charakteristischen Polynoms Rundungsfehlern unterworfen ist. Anstelle des exakten charakteristischen Polynoms p_n erhält man dabei ein Polynom r_n(λ), welches sich häufig nur sehr wenig von p_n unterscheidet.

Die Nullstellen von r_n seien μ₁,…,μ_η. Leider gilt nicht immer μ_j ≈ λ_j (wobei angenommen sei, daß die Nullstelle von r_n und p_n der Größe nach geordnet sind: \begin{eqnarray} \left. \left| {{\lambda}_{1}} \right|\ge \left| {{\lambda}_{2}} \right|\ge \cdots \ge \left[ {{\lambda}_{n}} \right),\left| {{\mu}_{2}} \right|\ge \cdots \ge \left| {{\mu}_{n}} \right| \right). \end{eqnarray} Schon kleine Änderungen der Koeffizienten a_j von p_n können zu großen Änderungen der Eigenwerte führen. Diese mögliche hohe Empfindlichkeit der Eigenwerte nimmt in der Regel mit wachsendem Polynomgrad zu. Daher ist das charakteristische Polynom als Hilfsmittel zur numerischen Berechnung von Eigenwerten nicht geeignet.

Der Satz von Abel (Abel, Satz von) impliziert, daß jedes Verfahren zum Lösen des Eigenwertproblems iterativ sein muß, d.h., Eigenwerte können nur als Grenzwert einer unendlichen Folge von Näherungslösungen berechnet werden.

Je nach Aufgabenstellung sind einmal Eigenwerte und Eigenvektor, ein anderes Mal nur Eigenwerte oder Eigenvektoren der Matrix Λ gesucht. Es gibt Problemstellungen, bei denen man alle Eigenwerte/Eigenvektoren benötigt, aber auch solche, bei denen man nur an einigen wenigen interessiert ist, z. B. an dem betragsgrößten Eigenwert und dem zugehörigen Eigenvektor. Die Wahl des geeigneten Verfahrens zur Lösung hängt nicht nur von der jeweiligen Aufgabenstellung ab, sondern auch davon, welche speziellen Eigenschaften die Matrix A hat. Wichtige Kriterien sind hier, ob A symmetrisch oder nichtsymmetrisch ist, ob A klein- oder großdimensional ist, ob A dicht oder sparse besetzt ist.

Zur numerischen Lösung des Eigenwertproblems wurden zahlreiche Verfahren entwickelt, welche sich in zwei Verfahrensklassen einteilen lassen. Transformationsmethoden sind Methoden, bei denen die Matrix A selbst in jedem Iterationsschritt verändert wird. Typischerweise wird in jedem Iterationsschritt eine nichtsinguläre (idealerweise unitäre) Transformationsmatrix X_j berechnet und die Folge von Matrizen \begin{eqnarray} {{A}_{j}}=X_{j}^{-1}{{A}_{j-1}}{{X}_{j}} \end{eqnarray} gebildet, wobei A₀ = A gesetzt wird. Da alle Iterationsmatrizen A_j zueinander ähnlich sind, haben sie alle dieselben Eigenwerte. Ziel der Iteration ist es, A in eine Matrix zu transformieren, von der die Eigenwerte abgelesen werden können.

Man könnte nun versuchen, die X_j so zu wählen, daß die Folge der A_j gegen die Jordansche Normalform von A konvergiert. Dies führt aber auf einen numerisch nicht stabilen Algorithmus. Stattdessen wird meist das folgende Resultat verwendet:

Es existiert zu jeder reellen Matrix A eine unitäre Matrix Q so, daß\begin{eqnarray} {{Q}^{H}}AQ \end{eqnarray}eine obere Dreiecksmatrix ist (bzw. eine Diagonalmatrix, falls A symmetrisch, in diesem Fall kann Q orthogonal gewählt werden).

Von der Diagonalen dieser Dreiecksmatrix lassen sich dann die Eigenwerte ablesen, die Spalten von Q geben Informationen über die zugehörigen Eigenvektoren (bzw. invarianten Unterräume).

Da die Matrix Q nicht direkt in endlich vielen Schritten berechenbar ist, versucht man nun, durch eine geeignete Folge von Iterationsschritten Q als unendliches Produkt von unitären Matrizen X_j aufzubauen. Dabei wählt man X_j in jedem Schritt so, daß A_j in einem gewissen Sinne näher an einer oberen Dreiecksmatrix ist als A_j−1. Praktisch muß die Berechnung nach endlich vielen Schritten abgebrochen werden, z. B., wenn die Elemente im unteren Dreieck der Iterationsmatrix A_j klein genug sind. Die beiden bekanntesten Verfahren dieser Art sind das Jacobi-Verfahren zur Lösung von symmetrischen Eigenwertproblemen und der QR-Algorithmus zur Lösung eines symmetrischen oder eines allgemeinen Eigenwertproblems. Beide Verfahren sind hauptsächlich für Matrizen kleiner bis mittlerer Größe geeignet. Sie berechnen alle Eigenwerte einer Matrix A.

Die zweite Verfahrensklasse sind Iterationsverfahren, bei denen A erhalten bleibt. Es werden Folgen von Eigenwertnäherungen \begin{eqnarray} {{\beta}_{1}}\to {{\beta}_{2}}\to \cdots \to \cdots {{\lambda}_{j}} \end{eqnarray} und/oder Folgen von Eigenvektornäherungen \begin{eqnarray} {{z}_{1}}\to {{z}_{2}}\to \cdots \to {{x}_{j}} \end{eqnarray} erzeugt, die gegen gewisse Eigenwerte λ_j bzw. Eigenvektoren x_j konvergieren. Typische Vertreter dieser Verfahrensklasse sind Vektoriterationsverfahren, wie Potenzmethode oder inverse Iteration. Hierbei werden durch wiederholte einfache Matrix-Vektor-Multiplikationen oder wiederholtes Lösen eines Gleichungssystems Approximationen an einzelne Eigenvektoren einer Matrix A berechnet. Dabei fallen auch Approximationen an den zugehörigen Eigenwert ab.

Bei der Potenzmethode wird das Eigenwert-Eigenvektor-Paar zum betragsgrößten Eigenwert approximiert, bei der inversen Iteration ein beliebiges Eigenwert-Eigenvektor-Paar. Bei beiden Verfahren wird stets ein Eigenvektor approximiert.

Falls weitere Eigenvektoren benötigt werden, kann man wie folgt vorgehen: Die Matrix A wird mittels Deflation in eine Matrix B überführt, deren Spektrum alle Eigenwerte von A außer dem bereits berechneten enthält. Anschließend wendet man dann die Potenzmethode bzw. inverse Iteration auf B an.

Eine andere Möglichkeit ist, den Ansatz der Potenzmethode bzw. der inversen Iteration wie folgt zu verallgemeinern: Man startet die Iteration statt mit einem einzelnen Startvektor x₀ mit einer Start- matrix X₀ ∈ ℂ^n×p, um p Eigenwerte und den zugehörigen invarianten Unterraum gleichzeitig zu bestimmen. Dies führt auf Unterraum-Iterationsmethoden. Speziell für symmetrische Matrizen entwickelt wurden die Varianten Rayleigh-Quotienten-Verfahren und Rayleigh-Ritz-Verfahren.

Ist die Matrix A sehr groß und sparse, und ist man nur an einigen Eigenwerten und zugehörigen Eigenvektoren (bzw. invarianten Unterräumen) interessiert, so sind für symmetrisches A das Lanczos-Verfahren, bzw. für beliebige Matrizen das Arnoldi-Verfahren die geeigneten Verfahren. Das Lanczos-Verfahren war ursprünglich ein Verfahren zur Transformation einer symmetrischen Matrix auf Tridiagonalgestalt. Kombiniert mit einer Methode zur Bestimmung von Eigenwerten und Eigenvektoren symmetrischer Tridiagonalmatrizen (z.B. dem QR-Algorithmus oder der Sturm- schen Kette) ist es ein geeignetes Verfahren zur Lösung des symmetrischen Eigenwertproblems für große sparse Matrizen. Dabei wird die Eigenschaft, daß die Eigenwerte der bei der Konstruktion als Zwischenergebnisse auftretenden Tridiagonalmatrizen kleinerer Dimension häufig schon gute Approximationen an Eigenwerte von A sind, sehr erfolgreich dazu genutzt, einige Eigenwerte und Eigenvektoren von symmetrischen großen sparsen Matrizen zu berechnen. Die sparse Besetzung der Matrix kann hierbei gut genutzt werden, da in jedem Schritt nur die Multiplikation der gegebenen Matrix mit einem Vektor erforderlich ist. Das Arnoldi-Verfahren ist eine Verallgemeinerung des Lanczos-Verfahrens für nichtsymmetrische Matrizen. Es reduziert die gegebene Matrix auf obere Hessenberg-Form. Auch hier sind die Eigenwerte der bei der Konstruktion als Zwischenergebnisse auftretenden oberen Hessenberg-Matrizen kleinerer Dimension häufig schon gute Approximationen an Eigenwerte von A. Beide Verfahren lassen sich interpretieren als Berechnung einer orthogonalen Basis {q₁, q₂,…,q_n} für den Krylow-Raum \begin{eqnarray} \{{{q}_{1,}}A{{q}_{1,}}{{A}^{2}}{{q}_{1,\,.\,.\,.\,,}}{{A}^{n-1}}{{q}_{1}} \}, \end{eqnarray} bzw. als Berechnung einer QR-Zerlegung der Krylow-Matrix \begin{eqnarray} \begin{array}{*{35}{l}} k(A,{{q}_{1,}}n) & =({{q}_{1,}}{{Aq}_{1,}}{{A}^{2}}{{q}_{1,\,.\,.\,.\,,}}{{A}^{n-1}}{{q}_{1}}) \\ {} & =({{q}_{1,}}{{q}_{2,\,.\,.\,.,}}{{q}_{n}})R=QR. \\ \end{array} \end{eqnarray} Dabei werden im allgemeinen nicht alle Vektoren q₁,…,q_n berechnet, sondern man stoppt nach wenigen Schritten mit Q_k = (q₁, q₂, …,q_k), wobei k ≪ n. Dann ist \( Q_{k}^{T}A{{Q}_{k}} \) eine obere Hessenberg-Matrix (bzw. eine Tridiagonalmatrix, wenn A symmetrisch ist), deren Eigenwerte gute Approximationen an die Eigenwerte von A sind.

In den Anwendungen treten häufig allgemeinere Eigenwertprobleme als das hier betrachtete Standardeigenwertproblem auf. Ein stabiles Verfahren zur Lösung des verallgemeinerten Eigenwertproblems \begin{eqnarray} Ax=\lambda Bx \end{eqnarray} mit A, B ∈ ℝ^n×n ist das QZ-Verfahren.

Bezüglich weiterer Verfahren und Verfahren für andere Eigenwertprobleme muß auf die Spezialliteratur verwiesen werden.