Technische Artikel

“Magische” Rekonstruktion: Compressed Sensing

Von Cleve Moler, MathWorks


Als ich zum ersten Mal von Compressed Sensing hörte, war ich skeptisch. Es gab Behauptungen, dass es die für die Darstellung von Signalen und Bildern erforderliche Datenmenge um enorme Größenordnungen verringert und anschließend die Originale exakt wiederherstellt. Aus dem Nyquist-Shannon-Abtasttheorem wusste ich, dass das unmöglich ist. Nachdem ich aber mehr über Compressed Sensing erfahren hatte, erkannte ich, dass sowohl die Behauptungen als auch das Theorem völlig korrekt sind, sofern bestimmte Bedingungen gelten.

Das Abtasttheorem nach Nyquist und Shannon besagt, dass ein Signal zu seiner exakten und eindeutigen Wiederherstellung mit mindestens der doppelten Frequenz abgetastet werden muss. Natürlich gilt dieses Theorem auch weiterhin; wenn Sie ein Byte in einem aus weißem Rauschen bestehenden Signal oder Bild auslassen, können Sie das Original nicht wiederherstellen. Die meisten interessanten Signale und Bilder bestehen aber nicht aus weißem Rauschen. Stellt man sie durch geeignete Basisfunktionen wie etwa trigonometrische Funktionen oder Wavelets dar, dann haben viele Signale nur relativ wenige von Null verschiedene Koeffizienten. In der Terminologie des Compressed (oder Compressive) Sensing sind sie dünn besetzt.

Das zugrunde liegende Matrizen-Problem

Natürlich war der Aspekt des Compressed Sensing, der zuerst meine Aufmerksamkeit erregt hat, die zugrunde liegende Matrizen-Mathematik. Ein Rohsignal oder -bild kann man als einen Vektor f mit Millionen von Komponenten auffassen. Wir nehmen an, dass f als Linearkombination bestimmter Basisfunktionen aufgefasst werden kann:

\[f = \psi c\]

Diese Basisfunktionen müssen für eine konkrete Anwendung geeignet sein. In unserem Beispiel ist \(\psi\) die diskrete Cosinus-Transformation. Wir nehmen außerdem an, dass die meisten Koeffizienten c praktisch gleich Null sind und c damit dünn besetzt ist. Zur Abtastung des Signals benötigen wir noch einen anderen linearen Operator,

\[b=\phi f\]

In unserem Beispiel besteht \(b\) aus einigen zufälligen Abtastwerten von \(f\), wodurch \(\phi\) eine Teilmenge der Zeilen des Identitäts-Operators darstellt. Es sind aber auch kompliziertere Abtastoperatoren möglich.

Zur Rekonstruktion des Signals müssen wir versuchen, die Koeffizienten wiederherzustellen, indem wir die Gleichung lösen:

\[Ax = b, \text{mit} A = \phi\psi\]

Haben wir erst die Koeffizienten, dann können wir das Signal selbst rekonstruieren durch Berechnung von:

\[f \approx \psi x\]

Da dies eine Kompression ist, ist \(A\) rechteckig mit viel mehr Spalten als Zeilen. Zur Berechnung der Koeffizienten \(x\) muss ein unterbestimmtes lineares Gleichungssystem, \(Ax = b\), gelöst werden. In dieser Situation gibt es viel mehr Unbekannte als Gleichungen. Der Schlüssel des fast schon magischen Rekonstruktionsprozesses ist die nichtlineare Regularisierung mit der \(\ell_1\)-Norm.

Eine sehr kleine Version dieser Situation habe ich 1990 in einem Cleve’s Corner unter dem Titel “The World’s Simplest Impossible Problem” beschrieben. Ich denke mir zwei Zahlen, deren Mittelwert gleich 3 ist. Wie lauten die Zahlen? Nachdem Sie sich beschwert haben, dass ich Ihnen nicht genug Informationen gegeben habe, antworten Sie vielleicht 2 und 4. Falls Sie das tun, dann wenden Sie unbewusst eine Art der Regularisierung an, die voraussetzt, dass das Ergebnis aus zwei konkreten Ganzzahlen besteht.

MATLAB® kann zwei verschiedene Antworten geben. In algebraischer Hinsicht ist das Problem ein 1-mal-2 lineares Gleichungssystem mit der Matrix

A = [1/2 1/2]

und der rechten Seite

b = 3

Wir wollen einen aus zwei Elementen bestehenden Vektor \(y\) finden, der die Gleichung \(Ay = b\) löst. Die Kleinste-Quadrate-Lösung wird durch die Pseudoinverse berechnet:

y = pinv(A)*b y = 3 3

Aber der Backslash-Operator erzeugt eine andere Lösung:

x = A\b x = 6 0

Beide Lösungen sind gültig, werden aber trotzdem von menschlichen Rätsellösern nur selten angeführt. Beachten Sie, dass die vom Backslash berechnete Lösung dünn besetzt ist; eine ihrer Komponenten ist gleich Null.

Eine größere Variante derselben Aufgabe

Die Wiederherstellung von Signalen oder Bildern ist eine umfangreichere Anwendung derselben Aufgabe. Vorgegeben sind tausende gewichtete Mittelwerte von mehreren Millionen Signal- oder Bildwerten. Unsere Aufgabe ist, das ursprüngliche Signal oder Bild zu rekonstruieren. Wir haben die komprimierte Probe \(b\) und müssen \(Ax = b\) lösen. Die Anzahl möglicher Lösungen ist gewaltig. Um die richtige herauspicken zu können, benötigt man Vektornormen. Die bekannte euklidische Abstandsnorm, \(\ell_2\), lautet

norm(x,2) = sqrt(sum(x.^2))

Die “Manhattan”-Norm, \(\ell_1\), ist nach der Reisezeit entlang eines von Straßen gebildeten quadratischen Gitters benannt.

norm(x,1) = sum(abs(x))

Die Notation l0 wird für eine Funktion verwendet, die nur die Zahl der von Null verschiedenen Komponenten zählt. Sie ist eigentlich keine Norm.

norm0(x) = sum(x~=0)

Diese Normen erlegen unserem unterbestimmten linearen Gleichungssystem drei verschiedene nichtlineare Regularisierungen oder Beschränkungen auf.

Finde unter allen \(x\), die \(Ax = b\) erfüllen, ein an \(x\), das \(\ell_{p}(x)\) minimiert.

Für unser 20 Jahre altes Beispiel ist \(\ell_{2}(y)\) kleiner als \(\ell_{2}(x)\), und \(\ell_{1}(y)\) und \(\ell_1(x)\) sind zufällig gleich. \(\ell_{0}(y)\) ist aber größer als \(\ell_{0}(x)\).

Dünne Besetzung und die \(\ell_1\)-Norm

Die Schlüssel zum Compressed Sensing sind dünne Besetzung und die \(\ell_1\)-Norm. Wenn die Darstellung des Ursprungssignals oder -bildes als Linearkombination der gewählten Basisfunktionen viele Koeffizienten enthält, die gleich Null sind, dann kann man das Signal oft exakt rekonstruieren. Im Prinzip würde man zur Berechnung dieser Rekonstruktion die von Null verschiedenen Koeffizienten mit \(\ell_0\) zählen. Die rechnerische Komplexität dieses kombinatorischen Problems lässt das aber in der Praxis gar nicht zu (es ist NP-schwer). Glücklicherweise kann, wie David Donoho, Emmanuel Candès und Terence Tao gezeigt haben, \(\ell_0\) durch \(\ell_1\) ersetzt werden. Sie erklären, dass die beiden Probleme „mit überwältigender Wahrscheinlichkeit“ dieselbe Lösung haben. Die Berechnung von \(\ell_1\) ist durchführbar, weil sie als lineares Programmierungsproblem formuliert und mit dem gebräuchlichen Simplex-Algorithmus oder modernen Innere-Punkt-Methoden gelöst werden kann.

Ein Beispiel

Die Basisfunktion für unser Beispiel ist die Diskrete Cosinus-Transformation (DCT). Das von der Taste „A“ auf einem Tonwahl-Telefon erzeugte Signal ist die Summe zweier sinusförmiger Töne mit unvereinbaren Frequenzen,

\[f(t) = \sin(1394 \pi t) + \sin(3266 \pi t)\]

Tasten wir diesen Ton für eine Achtelsekunde mit 40000 Hz ab, resultiert ein Vektor der Länge n = 5000. Das obere Diagramm in Abbildung 1 zeigt einen Ausschnitt dieses Signals sowie einige der m = 500 zufällig herausgegriffenen Abtastwerte. Das untere Diagramm zeigt die Koeffizienten \(c\), die durch Berechnung der inversen Diskreten Cosinus-Transformation von \(f\) erhalten wurden, mit zwei Spikes bei den korrekten Frequenzen. Weil die beiden Frequenzen unvereinbar sind, fällt dieses Signal nicht genau in den von den DCT-Basisfunktionen abgedeckten Bereich, weshalb es einige Dutzend signifikante von Null verschiedene Koeffizienten gibt.

nn10_cc_fig1_w.jpg
Abb. 1. Oben: Zufällige Stichproben des Ursprungssignals, das die „A“-Taste auf einem Tonwahltelefon erzeugt. Unten: Die inverse Diskrete Cosinus-Transformation des Signals.

Für unser Beispiel ist das komprimierte Signal ein Vektor b aus m zufälligen Stichproben des Ursprungssignals. Wir konstruieren eine Matrix \(A\), indem wir m Zeilen aus der n-mal-n-DCT-Matrix extrahieren

D = dct(eye(n,n)); A = D(k,:)

wobei k der Vektor der für die Stichprobe \(b\). Das resultierende lineare System, \(Ax = b\), hat die Dimension m-mal-n, also 500-mal-5000. Es gibt 10 Mal so viele Unbekannte wie Gleichungen.

Zur Rekonstruktion des Signals müssen wir diejenige Lösung für \(Ax = b\) finden, die die \(\ell_1\)-Norm von x minimiert. Dies ist ein nichtlineares Optimierungsproblem und es gibt eine ganze Reihe von MATLAB-basierten Programmen zu seiner Lösung. Ich habe mich für \(\ell_1\)-magic entschieden, das von Justin Romberg und Emmanuel Candès während ihrer Zeit am Caltech geschrieben wurde. Das obere Diagramm in Abbildung 2 zeigt die resultierende Lösung, \(x\). Wir sehen, dass sie relativ wenige große Komponenten besitzt und dass sie der DCT des Ursprungssignals sehr nahe kommt. Darüber hinaus sieht die im unteren Diagramm dargestellte Diskrete Cosinus-Transformation von \(x\) dem Ursprungssignal ebenfalls sehr ähnlich. Könnten wir hier Töne wiedergeben, dann würden Sie hören, dass die Audioausgaben der beiden Befehle sound(f) und sound(dct(x)) fast gleich klingen.

nn10_cc_fig2_w.jpg
Abb. 2. Die l1-Lösung für Ax = b erzeugt dct(x), ein Signal, das fast identisch mit dem Original ist.

Zum Vergleich zeigt Abbildung 3 die übliche l2- oder Kleinste-Quadrate-Lösung für \(Ay = b\), die berechnet wird durch

y = pinv(A)*b 
nn10_cc_fig3_w.jpg
Abb. 3. Die l2-Lösung für Ay = b erzeugt dct(y), ein Signal, das mit dem Original nur wenig Ähnlichkeit hat.

In den Diagrammen erkennt man einen schwachen Anflug des Ursprungssignals, aber sound(dct(y)) ist lediglich ein verrauschtes Brummen.

Es ist noch zu früh um vorauszusagen, wann – oder ob – wir Compressed Sensing in unseren Mobiltelefonen, Digitalkameras und Magnetresonanz-Scannern antreffen werden, aber ich finde die zugrunde liegende Mathematik und Software faszinierend.

Veröffentlicht 2010 - 91850v00

Literatur

Eingesetzte Produkte