DE4313627C1 - Verfahren zum Ermitteln der Eigenwerte einer Matrix aus Signalabtastwerten - Google Patents

Description

Die vorliegende Erfindung betrifft ein Verfahren zum Ermitteln der Eigenwerte einer Matrix, deren Elemente Abtastwerte eines Signals sind, wobei die Matrix der das Signal eindeutig charakterisierenden Eigenwerte durch Transformation der Matrix der Signalabtastwerte auf Hauptdiagonalform mit Hilfe von zyklischen Vektorrotationen nach dem Jacobi Verfahren ermittelt wird und die einzelnen Vektorrotationen mit Cordic Operationen realisiert werden.

Eigenwertberechnungen werden in der Signalverarbeitung angewendet, wenn es z. B. darum geht, das Spectrum eines Signals zu bestimmen oder Datenreduktionen bei Bildsignalen durchzuführen oder aus dem von einem zu lokalisierenden Objekt abgegebenen Signal und von mehreren Sensoren eines Antennenarrays empfangenen Signalproben die Richtung bzw. den Ort des Objekts zu schätzen. Vereinfacht gesagt, werden mit den Eigenwerten aus Abtastwerten eines Signals wenige charakteristische Werte abgeleitet, welche das Signal eindeutig beschreiben.

Die Eigenwerte Λ einer reellen symmetrischen n × n Matrix A, deren Elemente Abtastwerte (Proben) eines Signals sind, werden durch eine orthogonale Transformation

Λ = Q^TAQ, wobei Q^TQ = I (Einheitsmatrix)

bestimmt.

Um diese Transformation sehr schnell mit einem parallelarbeitenden Rechner ausführen zu können, wird die Jacobi Methode angewendet, wie z. B. bekannt aus R. P. Brent and F. T. Luk. The solution of singular value and symmetric eigenvalue problems on multiprocessor arrays, Society of Industrial and Applied Mathematics Journal, Scientific and Statistic Computing 6 (1985), 69-84.

Die Matrix Q besteht aus 2 × 2 Ebenen Rotationen, womit die Matrix A schrittweise auf Hauptdiagonalform reduziert wird.

Bei der klassischen Jacobi Methode werden folgende zwei Maßnahmen angewendet.

• a) das Verarbeitungsschema der einzelnen Elemente a_pq der Matrix A basiert darauf, daß die Zeilen und Spaltenindices p und q so gewählt werden, daß |a_pq| = max (|a_ÿ|), 1 i, j n, i = j
• b) Der Rotationswinkel wird iterativ bei jedem Transformationsschritt so gewählt, daß gilt: a_pq ^(k+1) = 0,wobei k der Iterationsindex ist.

Die klassische Jacobi Methode, welche gemäß a) die Indices p und q der Matrixelemente bestimmt, erschwert eine parallele Bearbeitungsprozedur. Außerdem sind für die Bestimmung der Rotationswinkel gemäß b) Operationen (Wurzelberechnungen, Divisionen) erforderlich, welche einen großen Nachteil für die Integration einer integrierten Schaltung darstellen.

Alle diese Nachteile werden vermieden, wenn - wie in der eingangs zitierten Veröffentlichung von R. P. Brent und F. T. Luk beschrieben - die Elemente der Nebendiagonalen der Matrix A zyklisch (d. h. in mehreren Iterationsschritten) Reihe für Reihe oder Spalte für Spalte einem Rotationsprozeß unterzogen werden. Wenn jedes Nebendiagonalelement einmal bearbeitet worden ist, ist ein Iterationsschritt beendet.

Für die Berechnung der Rotationswinkel ist es in Bezug auf die Integrierbarkeit der Schaltung vorteilhaft, das sogenannte Cordic (Coordinate Rotation Digital Computer)-Verfahren anzuwenden (vgl. z. B. J.-M. Delosme, Cordic algorithms: theory and extensions, Proc. SPIE Advanced Algorithms and Architectures for Signal Processing IV (1989), 131-145 oder EP 04 53 641 A2).

Mit dem Cordic-Verfahren werden üblicherweise die Rotationswinkel so bestimmt, daß die zweite Maßnahme b) des klassischen Jacobi-Verfahrens, nämlich jedes Nebendiagonalenelement a_pq ^(k+1) exakt zu Null zu machen, erfüllt wird. Doch dies ist nicht sehr sinnvoll, da ein in einem Arbeitsschritt zu Null transformiertes Element in einem nächsten Arbeitsschritt wieder einen von Null verschiedenen Wert annimmt. Aus diesem Grund wird in der Veröffentlichung von J.-P. Charlier, M. Vanbegin and P. van Dooren, On efficient implementations of Kogbetliantz's algorithm for computing the singular value decomposition. Numer. Math. 52 (1988), 279-300 vorgeschlagen, eine approximierte Rotation durchzuführen, die nicht a_pg ^(k+1) = 0 sondern nur

|a_pq ^(k+1)| < |a_pg ^(k)|

verlangt.

Bei dieser approximierten Rotationsmethode steigt zwar die Zahl der erforderlichen Iterationsschritte an aber der Rechenaufwand pron Iterationsschritt ist erheblich geringer als bei der exakten Rotationsmethode.

Daß approximierte Rotationen besonders vorteilhaft mit Cordic Prozessoren durchführbar sind, ist bereits von J.-M. Delosme, Bit-level systolic algorithms for real symmetric and hermitian eigenvalue problems. Journal of VLSI Signal Processing, 4(1): 69-88, 1992 dargelegt worden.

J.-M. Delosme schlägt vor, in jedem Rotationsprozeß nicht wie bisher eine fest vorgegebene Folge von Rotationswinkeln abzuarbeiten, sondern nur einige (3 oder 4) Rotationswinkel aus dieser Folge anzuwenden. Allerdings ist in dieser Literaturstelle nicht offenbart, nach welchen Kriterien diese Winkel ausgewählt werden.

Der Erfindung liegt nun die Aufgabe zugrunde, ein Verfahren der eingangs genannten Art anzugeben, welches für die Bestimmung der Eigenwerte mit einer möglichst geringen Zahl an Operationen (Additionen, Schiebeoperationen) auskommt, so daß eine diese Operationen ausführende Schaltung mit möglichst geringem Aufwand realisierbar ist.

Erfindungsgemäß wird diese Aufgabe durch die Merkmale des Hauptanspruchs gelöst. Eine vorteilhafte Weiterbildung der Erfindung geht aus dem Unteranspruch hervor.

Nachfolgend soll die Erfindung näher erläutert werden.

Zunächst wird die in der Beschreibungseinleitung schon erwähnte zyklische Jacobi Methode vorgestellt.

Nach der Jacobi Methode wird auf eine n × n Matrix A, die Abtastwerte eines Signals enthält, eine Folge von orthogonalen Transformationen

A^(k+1) = J^T _pq A^(k) J_pg (1)
(k = 0, 1, 2 . . . Iterationsindex)

durchgeführt, wobei J_pq eine Rotation in der p, q) Ebene ist, definiert durch die Parameter (c, s, -s, c) in den (pp, pq, qp, pp) Elementen einer n × n Einheitsmatrix.

c = cos Φ und s = sin Φ, Φ ist der Rotationswinkel.

Die Parameter (p, q) werden zyklisch ausgewählt, z. B. (1, 2), (1, 3) . . . (1, n) (2, 3) . . . (2, n) . . . (n-1, n).

Die Jacobi Rotationen J_pq sind so zu bestimmen, daß die Größe der Nebendiagonalelemente

mit jeder Transformation gemäß (1) abnimmt.

Es läßt sich zeigen, daß

[S^(k+1)]² = [S^(k)]² - [(a_pq)^(k)² - (a_pq)^(k+1)²] (3)

Die größte Reduktion von S^(k) wird erreicht bei a_pq ^(k+1): = 0.

Dies kann erreicht werden, wenn gilt:

mit

und daher

ist. Die vorausgehende Methode wird exakte Jacobirotation genannt wegen der Tatsache, daß a_pq ^(k) exakt zu Null gemacht wird.

Für die Reduktion von S^(k) ist es nicht notwendig, die exakte Jacobi Rotation, also (c, s) nach (5), anzuwenden.

Es reicht aus, (c, s) näherungsweise so zu bestimmen, daß die Orthogonalität erhalten bleibt und

|a_pq ^(k+1)| d · a_pq ^(k)| mit 0 |d| < 1 (6)

ist.

ln diesem Fall spricht man von einer approximierten Jacobi Rotatioin _pq.

Wegen

a_pq ^(k+1) = d · a_pq ^(k) mit d = c² - s² - 2τcs

gilt:

Der maximale Wert von |d| ist ein Maß für die Qualität der Approximation. Es läßt sich zeigen, daß für die exakte Jakobi Rotation, für die

d = 0 ist.

Es läßt sich auch zeigen, daß A^(k) zu einer Diagonalmatrix konvergiert, welche die Eigenwerte von A enthält, falls die verwendete Approximation |d(τ) | < 1 garantiert. Ferner wird quadratische Konvergenz erreicht, wenn |d(τ)| während des Verfahrens abnimmt, so daß gilt |d(τ)|_{τ →∞} → 0.

Wie bereits oben ausgeführt, wird die Jacobi Rotation nach dem Cordic-Verfahren ausgeführt. Deshalb sollen zunächst einige grundsätzliche Erläuterung zum Cordic-Verfahren folgen.

Eine Rotation eines Vektors (x, y)^T um einen Winkel Φ kann so ausgedrückt werden:

Diese Matrix-Vektor Multiplikation kann ausgeführt werden durch eine Folge von elementaren Rotationen Φ_k, welche durch einfache Schiebe- und Additionsoperationen realisierbar ist, falls Φ_k so gewählt wird, daß

Φ_k = ± arctan 2^-k.

Mit cos arctan 2^-k = 1/ gilt:

Diese Matrix-Vektor Multiplikation erfordert nur zwei Additionen und zwei Schiebeoperationen von k Bits. Eine Elementar Cordicoperation ist also die Matrix-Vektor- Multiplikation:

Das Ergebnis dieser Elementar-Operation wird mit dem Faktor K_k = skaliert.

Die Cordic Prozedur benutzt den Grundwinkel ± arctan 2^-k, um daraus den gewünschten Rotationswinkel Φ zu bilden.

Dies kann interpretiert werden als eine Darstellung des Rotationswinkels Φ in einem "arctan 2^-k" Zahlensystem mit den Ziffern f_kε {-1,1}. Damit ist nach b+1 Iterationsschritten der Eingangsvektor (x,y)^T um den Winkel Φ mit einer Genauigkeit von b Bits rotiert. Die Iterationsgleichungen für die Cordic Operation in zirkularen Koordinaten sind mit k = 0, 1, . . ., b:

X_k+1 = x_k + sign (τ_k) y_k2^-k
Y_k+1 = Y_k - sign (τ_k) x_k2^-k
τ_k+1 = τ_k - sign (τ_k) arctan 2^-k (11)

Während der Cordic Interationen geht der Parameter τ gegen Null und der Vektor (x, y)^T wird gedreht und gestreckt zu K_p (x′, y′)^T.

Der Skalierungsfaktor K_k eines Iterationsschrittes k hängt nicht von sign (τ_k) sondern nur von der Zahl der Iterationsschritte b ab, weil der Cosinus eine gerade Funktion ist. Der Vektor nach b+1 Iterationsschritten ist mit dem Faktor

zu skalieren, unabhängig vom Rotationwinkel.

In Matrixschreibweise läßt sich das Cordic Verfahren nun so darstellen.

Würde man diesen Cordic Prozeß mit b Iterationsschritten durchführen, wären insgesamt sehr viele Schiebe- und Additionsoperationen durchzuführen. Die Zahl der erforderlichen Operationen läßt sich sehr stark verringern, wenn in jedem Cordic Rotationsprozeß nur ein einziger Winkel (1 Iterationsschritt) ermittelt wird, der dem optimalen Rotationswinkel möglichst nahe kommt und nicht wie beim Stand der Technik in einem mehrschrittigen Iterationsprozeß man sich dem gewünschten Rotationswinkel über eine Folge von Winkel schrittweise annähert.

Man setzt für t aus Gleichung (4)

= sign (τ) 2^-l ≈ t (13)

wobei l = 0, 1, 2, . . ., b ist. l wird nun so gewählt, daß der Winkel = arctan 2^-l möglichst wenig vom optimalen Rotationswinkel Φ abweicht.

Die approximierte Jacobi Rotation ist:

Um _pq (l) zu bestimmen soll, wie gesagt, nur noch ein Iterationsschritt angewendet werden.

Für die Bestimmung der Jacobi Rotation _pq (l) reicht es aus, nur den Wert l zu bestimmen. Für die Bestimmung von l und für die Minimierung des maximalen Wertes von |d| ist d (|τ|, l) gemäß Gleichung (7) mit t: = sign (τ 2^-l:

Falls l = i für |τ|∈[τ_i, τ_i+a] gesetzt wird, läßt sich zeigen, daß |d(|τ|,l)| 1/3 ist. Die Grenzen des Intervalls für τ_i ergeben sich aus (15) mit |d(|τ_i|, l)| = 1/3 und l = i

|τ_i| = 1/3 (2ⁱ - 2^-i+1) (16)

Der Vergleich |τ_i| < |τ| < |τ_i+1| für die Bestimmung von 1 kann mit Schiebe- und Additionsoperationen durchgeführt werden. Mit a_D = |a_qq - a_pp| erhält man den folgenden Vergleich für die Bestimmung von l.

Falls

(2ⁱ - 2^-i+1)|a_pq| a_D+2^-1 a_D < (2ⁱ⁺¹ - 2^-i)|a_pq|, (17)

dann ist l = i.

Um eine große Zahl solcher Vergleiche zu vermeiden, wird das Intervall [τi, τi+1], d. h. l geschätzt.

Dazu wird ein Schätzwert von |τ| berechnet. Mit exp (a), der den Exponenten von a angibt, und mit |τ| = a_D/2|a_pq| erhält man

τ∈[2^k, 2^k+2 [mit k = exp (a_D) - exp (a_pq). (18)

Der Vergleich dieser Intervalle [2^k, 2^k+2], die einen Schätzwert für |τ| bedeuten, mit den Intervallen [τⁱ τⁱ⁺¹] (τi aus (16)), welche die Werte für l(l = i) bestimmen, zeigt, daß maximal drei Werte für l für ein gegebenes k möglich sind.

Deshalb, falls Gleichung (18) ein k < 0 ergibt, sind für die Bestimmung von l maximal drei Vergleiche gemäß (17) mit i = k-1, k, k+1 auszuführen.

Für die Approximation des Skalierungsfaktors für die Jacobi Rotation wird der nächst kleinere Rotationswinkel genommen, d. h. es wird 1: = 1+1 gesetzt, und die Rotatio _pq (1) gemäß (14) wird zweimal angewendet.

Der Skalierungsfaktor ergibt sich dann mit z. B. b = 32 aus

für jede Rotation gemäß (20).

Wenn das zuvor beschriebene Prinzip der Rotation mit einem Winkel angewendet wird, ist zwar die Konvergenz bei der Iteration nur linear (nichtquadratisch), d. h. die Anzahl der Iterationsschnitte steigt an. Der Aufwand an Schiebe- und Additionsoperationen ist aber trotzdem geringer als wenn jede Rotation in mehreren Winkelschritten durchgeführt würde und dadurch weniger Jacobi Iterationsschritte erforderlich wären (wegen quadratischer Konvergenz). Das Verhältnis von zusätzlichen Jacobi Iterationsschritten zu eingesparten Operationen ist kleiner 1.

Die quadratische Konvergenz kann aber trotz der Anwendung der Rotation mit nur einem Winkel erreicht werden, indem nämlich die Transformation mit der "Ein-Winkel-Rotation" ein oder mehrmals wiederholt wird. Die Genauigkeit der Rotation wird dabei von Mal zu Mal verbessert, wodurch schließlich quadratische Konvergent erreicht wird.

Nachfolgend ist an einem Beispiel (2 × 2 Matrix A) mit m = 1 . . . 5 wiederholten Transformationen, wobei in jeder die "Ein-Winkel-Rotation" angewendet worden ist, gezeigt.

Dieses Beispiel zeigt, daß mit nur 5 Cordic Rotationswinkeln Φ die exakte Transformation (Nebendiagonalelement von A sind Null) erreicht wird. Um dasselbe Ergebnis mit dem bekannten Transformationsverfahren zu erreichen, hätte man insgesamt b = 16 Cordic Rotationswinkel benötigt.

Claims (2)

1. Verfahren zum Ermitteln der Eigenwerte einer Matrix, deren Elemente Abtastwerte eines Signals sind, wobei die Matrix der das Signal eindeutig charakterisierenden Eigenwerte durch Transformation der Matrix der Signalabtastwerte auf Hauptdiagonalform mit Hilfe von zyklischen Vektorrotationen nach dem Jacobi Verfahren ermittelt wird und die einzelnen Vektorrotationen mit Cordic Operation realisiert werden, dadurch gekennzeichnet,
daß jede Vektorrotation mit einem einzigen Cordic Winkel vorgenommen wird und daß der Cordic Winkel jeweils so gewählt wird, daß durch die Transformation jedes Element der Matrix außerhalb der Hauptdiagonalen minimal wird.

2. Verfahren nach Anspruch 1, dadurch gekennzeichnet, daß die Transformation mindestens einmal wiederholt wird, wobei für jede Jacobi Rotation ein neuer Cordic Winkel verwendet wird, so daß mit jedem Transformationsprozeß die Elemente der Matrix außerhalb der Hauptdiagonalen weiter gegen Null konvergieren.