Algorithmus von Samuelson-Berkowitz

Algorithmus von Samuelson-Berkowitz

Der Algorithmus von Samuelson-Berkowitz (nach Paul A. Samuelson und S. Berkowitz) ist ein Verfahren, das für beliebige quadratische Eingabematrizen  A\in R^{n\times n} die Koeffizienten des charakteristischen Polynoms von A ermittelt, d.h. insbesondere auch die Determinante von A. Im Gegensatz etwa zum Algorithmus von Faddejew-Leverrier sind die Voraussetzungen weniger restriktiv: Als Eingabe sind auch Matrizen A zulässig, deren Einträge Elemente eines beliebigen kommutativen Rings R mit Einselement sind, so dass das Verfahren völlig ohne Divisionen auskommt.

Inhaltsverzeichnis

Notation und Idee des Verfahrens

Wir bezeichnen mit

  •  I_r \in R^{r\times r} die r\times r-Einheitsmatrix
  •  A_r \in R^{r\times r} \; die r\times r-Submatrix von A bestehend aus den ersten r Zeilen und Spalten
  •  p_r \in R[\lambda] \; das charakteristische Polynom von  A_r\; , wobei p_r(\lambda)=\sum\limits_{k=0}^r c_{r-k} \lambda^k
  •  R_r \in R^{r} \; der Zeilenvektor mit den Komponenten  a_{r+1,j} \; mit  1 \leq j \leq r
  •  S_r \in R^{r} \; der Spaltenvektor mit den Komponenten  a_{i,r+1} \; mit  1 \leq i \leq r

und betrachten folgende Partitionierung von Ar + 1:

 A_{r+1} = \left[ \begin{array}{c|c}
                    A_r & S_r \\ \hline
                    R_r & a_{r+1,r+1}
                    \end{array}  \right]

Die grundlegende Idee des Verfahrens besteht darin, das charakteristische Polynom von  A_{r+1}\; rekursiv zu berechnen. Mit der obigen Notation gilt zunächst

 \det(A_{r+1}) = a_{r+1,r+1} \det(A_r) - R_r \; \textrm{Adj}(A_r) \; S_r \;

wobei Adj(Ar) die Adjunkte von  A_r\; bezeichnet (Begründung: Entwicklung nach letzter Zeile mittels Entwicklungssatz von Laplace, vgl. [1]).

Speziell für das charakteristische Polynom von  A_{r+1}\; bedeutet das also:

 \begin{align}
         p_{r+1}(\lambda) &= \det(\lambda I_{r+1} - A_{r+1}) \\
                 &= ( \lambda - a_{r+1,r+1}  ) \; p_r(\lambda) - R_r \; \textrm{Adj}(\lambda I_r - A_r ) \; S_r \qquad (*) \end{align}

Außerdem kann man leicht zeigen, dass sich die Adjunkte in (*) folgendermaßen schreiben lässt (siehe z.B. [1]):

 \begin{align}
\textrm{Adj}(\lambda I_r - A_r) &= \sum\limits_{k=1}^r \left( \sum\limits_{j=1}^k A_r^{j-1} c_{k-j} \right) \lambda^{r-k}\\
                                &= \sum\limits_{k=1}^r \left( A_r^{k-1}c_0 + \ldots + A_r^1 c_{k-2} + I_r c_{k-1} \right) \; \lambda^{r-k} \qquad (**) 
\end{align}

Hierin sind  c_0, \ldots, c_{r-1} die Koeffizienten des charakteristischen Polynoms von Ar.

Formel von Samuelson

Wir erhalten nun die gewünschte rekursive Darstellung für das charakteristische Polynom  p_{r+1}(\lambda) \; von  A_{r+1} \; (in der Literatur Formel von Samuelson genannt), indem wir die beiden obigen Beziehungen (*) und (**) zusammenfügen:

 \begin{align}
p_{r+1}(\lambda) &= (\lambda - a_{r+1,r+1}) \; p_r(\lambda) - \sum_{k=1}^r \left[ \sum_{j=1}^k (R_r A_r^{j-1} S_r) c_{k-j} \right] \lambda^{r-k} \\
                 &= (\lambda - a_{r+1,r+1}) \; p_r(\lambda) - \sum_{k=1}^r \left[ (R_r A_r^{k-1} S_r) c_0 + \ldots + R_r A_r^1 S_r c_{k-2} + (R_rS_r)c_{k-1}  \right] \lambda^{r-k}
\end{align}

Verfahren von Samuelson-Berkowitz

Matrix-Vektor Schreibweise

Um einen effektiven und leichter lesbaren Algorithmus formulieren zu können, transferieren wir nun die Formel von Samuelson in Matrix-Schreibweise. Dazu ordnen wir einem Polynom  \omega\; vom Grad d

 \omega(\lambda) = \sum_{k=0}^d \alpha_k \lambda^{d-k}

den Koeffizientenvektor

 \overrightarrow{\omega} = \left[ \begin{array}{c} \alpha_0 \\ \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_d \end{array} \right] \in R^{d+1}

sowie die folgende Toeplitz-Matrix (die zugleich eine untere Dreiecksmatrix ist) zu:

 \textrm{Toep}(\omega) = \left[ \begin{array}{ccccc}
                           \alpha_0 & 0   & 0 & \ldots & 0 \\
                           \alpha_1 & \alpha_0 & 0 & \ldots & 0 \\
                           \alpha_2 & \alpha_1 & \alpha_0 & \ldots & 0 \\
                           \vdots & \vdots & \vdots & \ddots & \vdots \\
                           \cdot  & \cdot  & \cdot  & \ldots & a_0 \\
                           \alpha_d    & \alpha_{d-1}& \alpha_{d-2} & \ldots & a_1 
                           \end{array} \right] \in R^{(d+1)\times d}

Genauer ist also der Eintrag an der Position (i,j) von Toep(ω) gegeben durch

 \left[ \textrm{Toep}(\omega) \right]_{i,j} = \left\{ \begin{array}{cl}
                                                     \alpha_k & \textrm{falls}\;i\ge j \; \textrm{und} \; i-j=k \\
                                                     0        & \textrm{falls}\;i<j    
                                                     \end{array}\right.

Definiert man nun noch das Polynom  q_{r+1} \; durch

 \begin{align}
q_{r+1}(\lambda) &= \lambda^{r+1} - a_{r+1,r+1} \lambda^r - \sum\limits_{k=1}^r (R_r A^{k-1} S_r) \lambda^{r-k} \\
                 &= \lambda^{r+1} - a_{r+1,r+1} \lambda^r - (R_rS_r) \lambda^{r-1} - \ldots - (R_r A_r^{r-1}S_r) \end{align}

dann lässt sich die Formel von Samuelson in der folgenden kompakten Form darstellen (vgl. [2] und [3]):

 \overrightarrow{p_{r+1}} = \textrm{Toep}(q_{r+1}) \overrightarrow{p_r}

Algebraische Top-Level Formulierung

Durch sukzessives Anwenden dieses Prinzips erhält man folgende zentrale Aussage (siehe [2] und [3]):

Mit  p_1(\lambda)=\lambda - a_{11} \; , also

 \overrightarrow{p_1} = \left[ \begin{array}{c} 1 \\ -a_{11} \end{array} \right] = \textrm{Toep}(q_1)

gilt:

  • Die Koeffizienten des charakteristischen Polynoms  p_r(\lambda) \; von  A_r \; für  1\leq r \leq n sind gegeben durch:
 \overrightarrow{p_r} = \textrm{Toep}(q_{r}) \cdot \textrm{Toep}(q_{r-1}) \cdots \textrm{Toep}(q_{1})
  • Insbesondere erhalten wir, falls r = n, die Koeffizienten des charakteristischen Polynoms  p_n(\lambda) \; von  A \; durch:
 \overrightarrow{p_n} = \textrm{Toep}(q_{n}) \cdot \textrm{Toep}(q_{n-1}) \cdots \textrm{Toep}(q_{1})

Algorithmus

Damit kann man nun folgenden Algorithmus formulieren (vgl. [4]):

 \textrm{Eingabe:} \;\; (n \times n) - \textrm{Matrix} \;  A\in R^{n\times n} 
 \overrightarrow{\textrm{Vect}} := \left[ \begin{array}{c} 1 \\ -a_{11} \end{array} \right] 
 \textbf{For} \; r=1,\ldots, n-1 \;  \textrm{berechne} 
   
  *  \left\{-R_r A_r^{k-1} S_r \right\}_{k=1}^r \; \textrm{(Koeffizienten}\;\textrm{von}\; q_{r+1}\;\textrm{)} 
  *  \overrightarrow{\textrm{Vect}} := \textrm{Toep}(q_{r+1}) \cdot \overrightarrow{\textrm{Vect}} \;
 \textbf{End}\;\textbf{For} 
 \textrm{Ausgabe:} \;\;  \overrightarrow{p_n} = \overrightarrow{\textrm{Vect}} \;\;\textrm{(Vektor}\;\textrm{der}\;\textrm{Koeffizienten}\;\textrm{des}\;\textrm{charakteristischen}\;\textrm{Polynoms)}

Der Algorithmus berechnet nicht nur die Koeffizienten des charakteristischen Polynoms von A, sondern darüber hinaus auch in jedem Schleifendurchlauf für die Untermatrix Ar.

Historisches

Die grundlegende Idee des Verfahrens wurde zuerst 1942 von Paul A. Samuelson beschrieben und publiziert [5]. Der Algorithmus in der oben präsentierten und heute gebräuchlichen Form geht auf Berkowitz [3] (parallele Version) und Abedeljaoued [2] (Beschreibung als serielles Verfahren) zurück, weswegen man manchmal auch die Bezeichnung Samuelson-Berkowitz-Abdeljaoued-Algorithmus (SBA-Algorithmus) in der Literatur findet [4].

Korrektheit des Algorithmus

Da im oben formulierten Verfahren nur endliche Schleifen auftreten, ist klar dass der Algorithmus terminiert. Die partielle Korrektheit folgt aus der Formel von Samuelson und der daraus abgeleiteten algebraischen Top-Level-Formulierung in Matrix-Vektor-Form (s.o., vgl. z.B. [1]). Genauer gesprochen beruht die Korrektheit auf folgender Schleifeninvariante: Am Ende des r ten Schleifendurchlaufs enthält der Vektor \overrightarrow{\textrm{Vect}} die Koeffizienten des charakteristischen Polynoms von Ar + 1(Formulierung als Nachbedingung).

Aufwand, Effizienz und Parallelisierbarkeit

Man kann zeigen (siehe [2], dass der Aufwand (Zeitkomplexität) des Algorithmus die Größenordnung  \mathcal{O}(n^4) hat. Eine genauere Schranke ist gegeben durch die Anzahl der arithmetischen Operationen  \frac 1 2 n^4 - n^3 + \frac 5 2 n^2 . Bei der Implementierung des Verfahrens kann man zudem ausnutzen, dass es für die Multiplikation von Toeplitz-Matrizen effektive Methoden gibt. Der Algorithmus lässt sich auch sehr gut parallelisieren, genaueres dazu findet man speziell in [3].

Numerisches Beispiel

Wir betrachten die Matrix


 A=\left[ \begin{array}{rrrr}
          5 &  5 & -3 & -7 \\
          2 &  1 &  9 &  6 \\
          4 &  2 & -6 & -5 \\
          5 & -8 & -9 &  2
          \end{array} \right]
Wir starten die Rekursion mit den charakteristischen Polynom der Matrix A1 = [5], für das  p_1 = \lambda - 5 \; gilt, d.h.
 \overrightarrow{\textrm{Vect}} = \left[ \begin{array}{r} 1 \\ -5 \end{array} \right]
  • r = 1:
Wir berechnen nun p2. Hierzu benötigen wir zunächst die Koeffizienten von q2:
  • k = 1:
 -R_1S_1 = -2*5=-10 \;
Also
 \begin{align}
         q_2 &= \lambda^2 - a_{22}\lambda - R_1S_1 \\
             &= \lambda^2 - 1 \lambda - 10\; 
        \end{align}
Hieraus resultiert nun die Toeplitz-Matrix
 \textrm{Toep}(q_2)=
 \left[ \begin{array}{rr}
           1 &  0 \\
         - 1 &  1 \\
         -10 & -1 
        \end{array} \right]
und damit
 \begin{align} 
        \textrm{Toep}(q_2) * \overrightarrow{p_1} &= \overrightarrow{\textrm{Vect}} \\
        \left[ \begin{array}{rr}
           1 &  0 \\
         - 1 &  1 \\
         -10 & -1 
        \end{array} \right] *
        \left[ \begin{array}{r}
        1 \\ -5
        \end{array} \right]
         &=
        \left[ \begin{array}{r}
           1  \\
         - 6  \\
         - 5  
        \end{array} \right] \end{align}
Das charakteristische Polynom von A2 lautet also  p_2=\lambda^2 -6 \lambda - 5 \;
  • r = 2:
Wir ermitteln die Koeffizienten von q3:
  • k = 1:
 - R_2 A_2^0 S_2 = - \left[4 \;\; 2\right] \cdot \left[ \begin{array}{r} -3 \\ 9 \end{array} \right] = -6
  • k = 2:
 -R_2 A_2^1 S_2 =  - \left[4 \;\; 2\right] \cdot \left[ \begin{array}{rr} 5 & 5 \\ 2 & 1 \end{array} \right] \cdot \left[ \begin{array}{r} -3 \\ 9 \end{array} \right] = -126
Also
 \begin{align}
        q_{3}(\lambda) &= \lambda^{3} - a_{3,3} \lambda^2 - (R_2S_2) \lambda^{1} - (R_2 A_2^{1}S_2) \\
                       &= \lambda^{3} + 6\lambda^2 - 6\lambda - 126
        \end{align}
und

\textrm{Toep}(q_3)=
 \left[ \begin{array}{rrr}
        1    &  0  & 0 \\
        6    &  1  & 0 \\
       -6    &  6  & 1 \\
       -126  & -6  & 6 
        \end{array} \right]
Damit erhalten wir

\begin{align} 
        \textrm{Toep}(q_3) * \overrightarrow{p_2} &= \overrightarrow{\textrm{Vect}} \\
\left[ \begin{array}{rrr}
        1    &  0  & 0 \\
        6    &  1  & 0 \\
       -6    &  6  & 1 \\
       -126  & -6  & 6 
        \end{array} \right] \cdot
        \left[ \begin{array}{r}
           1  \\
         - 6  \\
         - 5  
        \end{array} \right] &=
        \left[ \begin{array}{r}
         1 \\
         0 \\
        -47 \\
        -120  
        \end{array} \right]
\end{align}
Das charakteristische Polynom von A3 lautet daher  p_3=\lambda^3 -47 \lambda -120 \;
  • r = 3:
Wir ermitteln die Koeffizienten von q4:
  • k = 1:
 - R_3 A_3^0 S_3 = - \left[5 \;\; -8 \;\; -9\right] \cdot \left[ \begin{array}{r} -7 \\ 6 \\ -5 \end{array} \right] = 38
  • k = 2:
 - R_3 A_3^1 S_3 = - \left[5 \;\; -8 \;\; -9\right] \cdot 
           \left[\begin{array}{rrr}
           5 &  5 & -3 \\
           2 &  1 &  9 \\
           4 &  2 & -6
           \end{array}\right] \cdot
           \left[ \begin{array}{r} -7 \\ 6 \\ -5 \end{array} \right] = -348
  • k = 3:
 - R_3 A_3^2 S_3 = - \left[5 \;\; -8 \;\; -9\right] \cdot 
           \left[\begin{array}{rrr}
           5 &  5 & -3 \\
           2 &  1 &  9 \\
           4 &  2 & -6
           \end{array}\right]^2 \cdot
           \left[ \begin{array}{r} -7 \\ 6 \\ -5 \end{array} \right] = 679
Also
 \begin{align}
        q_{4}(\lambda) &= \lambda^{4} - a_{4,4} \lambda^3 - (R_3S_3) \lambda^{2} - (R_3 A_2^{1}S_3) \lambda - (R_3 A_2^{2}S_3) \\
                       &= \lambda^{4} -2 \lambda^3 +38 \lambda^2 - 348 \lambda + 679
        \end{align}
und

\textrm{Toep}(q_4)=
 \left[ \begin{array}{rrrr}
        1    &    0 &  0 &  0 \\
       -2    &    1 &  0 &  0 \\
        38   &   -2 &  1 &  0 \\
       -348  &   38 & -2 &  1 \\
        679  & -348 & 38 & -2
        \end{array} \right]
Die finale Matrix-Vektor-Multiplikation liefert nun die Koeffizienten des gesuchten charakteristischen Polynoms der gesamten Matrix A = A4:

\begin{align} 
        \textrm{Toep}(q_4) * \overrightarrow{p_3} &= \overrightarrow{\textrm{Vect}} \\
 \left[ \begin{array}{rrrr}
        1    &    0 &  0 &  0 \\
       -2    &    1 &  0 &  0 \\
        38   &   -2 &  1 &  0 \\
       -348  &   38 & -2 &  1 \\
        679  & -348 & 38 & -2
        \end{array} \right] \cdot
        \left[ \begin{array}{r}
         1 \\
         0 \\
        -47 \\
        -120  
        \end{array} \right] &=
       \left[ \begin{array}{r}
        1 \\ -2 \\ -9 \\ -374 \\ -867        
        \end{array} \right] \end{align}
Hieraus liest man das gesuchte Endergebnis ab:
 p_4 = \lambda^4 -2\lambda^3 -9 \lambda^2 - 374 \lambda - 867 \;
Insbesondere erhält man also für die Determinante von A
 p_4(0) = \det(-A) = (-1)^4\cdot \det(A) = -867  \; \Rightarrow \det (A) = -867

Literatur

  • J. Abdeljaoued, The Berkowitz algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring, MapleTech Vol. 4, No. 3, pp. 21-32, Birkhäuser Boston Basel Berlin, 1997
  • Stuart J. Berkowitz: On computing the determinant in small parallel time using a small number of processors, Information Processing Letters, 18, pp. 147-150, 1985, doi:10.1016/0020-0190(84)90018-8
  • G. Nakos and Robert M. Williams: A Fast Computation of the Characteristic polynomial, Mathematica in Education and Research, Vol. 9, No. 1, 2000
  • Paul A. Samuelson: A method for determining explicitly the characteristic equation, Annals of Mathematical Statistics, 13, pp. 424-429, 1942, doi:10.1214/aoms/1177731540
  • Günter Rote: Division-free algorithms for the determinant and the Pfaffian: algebraic and combinatorial approaches , in: Computational Discrete Mathematics, Editor: Helmut Alt, Lecture Notes in Computer Science 2122, Springer-Verlag, 2001, pp. 119-135, Online-Version
  • Michael Soltys: Berkowitz's Algorithm and Clow Sequences, The Electronic Journal of Linear Algebra (ELA), ISSN 1081-3810, Volume 9, pp. 42-54, April 2002, Online-Version

Einzelnachweise

  1. a b c Michael Soltys: Berkowitz's Algorithm and Clow Sequences, The Electronic Journal of Linear Algebra (ELA), ISSN 1081-3810, Volume 9, pp. 42-54, April 2002, Online-Version
  2. a b c d J. Abdeljaoued, The Berkowitz algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring, MapleTech Vol. 4, No. 3, pp. 21-32, Birkhäuser Boston Basel Berlin, 1997
  3. a b c d Stuart J. Berkowitz: On computing the determinant in small parallel time using a small number of processors, Information Processing Letters, 18, pp. 147-150, 1985, doi:10.1016/0020-0190(84)90018-8
  4. a b G. Nakos and Robert M. Williams: A Fast Computation of the Characteristic polynomial, Mathematica in Education and Research, Vol. 9, No. 1, 2000
  5. Paul A. Samuelson: A method for determining explicitly the characteristic equation, Annals of Mathematical Statistics, 13, pp. 424-429, 1942, doi:10.1214/aoms/1177731540

Wikimedia Foundation.

Schlagen Sie auch in anderen Wörterbüchern nach:

  • Algorithmus von Faddejew-Leverrier — Der Algorithmus von Faddejew Leverrier (nach Dmitri Konstantinowitsch Faddejew und Urbain Le Verrier) ist ein Verfahren, das für beliebige quadratische Matrizen die Koeffizienten des durch p(λ) = det(λI − A) definierten charakteristischen… …   Deutsch Wikipedia

  • Charakteristisches Polynom — Das charakteristische Polynom (Abk.: CP) ist ein Begriff aus dem mathematischen Teilgebiet der linearen Algebra. Dieses Polynom, das für quadratische Matrizen und Endomorphismen von endlichdimensionalen Vektorräumen definiert ist, gibt Auskunft… …   Deutsch Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”