MF_SVdecompose MD_SVdecompose ME_SVdecompose
FunktionSingulärwert-Zerlegung (engl. Singular Value Decomposition, SVD)
Syntax C/C++#include <MFstd.h>
int MF_SVdecompose( fMatrix U, fMatrix V, fVector W, fMatrix MA, ui htA, ui lenA );
C++ MatObj#include <OptiVec.h>
int matrix<T>::SVdecompose( matrix<T> MV, vector<T> W, const matrix<T>& MA );
int matrix<T>::SVdecompose( matrix<T>* MV, vector<T>* W, const matrix<T>& MA );
Pascal/Delphiuses MFstd;
function MF_SVdecompose( MU, MV:fMatrix; W:fVector; MA:fMatrix; htA, lenA:UIntSize ): IntBool;
BeschreibungDie Matrix A mit den Dimensionen [ht, len] wird in ein Produkt
MA = MU * MW * MVT,
faktorisiert, wobei MU die Dimensionen [max(ht,len), len] hat. MW ist eine Diagonalmatrix, deren Elemente alle positiv oder 0 sind. Tatsächlich wird nur die Diagonale dieser Matrix in dem Vektor W der Größe [len] gespeichert. MV schließlich ist eine Quadratmatrix [len, len]. Sowohl MU als auch MV sind orthogonal:
MUT * MU = MVT * MV = (1).
Aufgrund dieser Orthogonalitäts-Beziehungen ist die Lösung eines linearen Gleichungssystem MA * X = B sehr einfach, hat man die beschriebene Faktorisierung erst einmal durchgeführt, da lediglich MW explizit invertiert werden muß:
X = MV * W-1 * MUT
(Diese Gleichung ist von rechts nach links zu berechnen.) Die Invertierung von MW ist auch wiederum einfach, da MW eine Diagonalmatrix ist und ihre Inverse aus den Kehrwerten der Diagonalelemente besteht. Diese Eigenschaft ist es, die SVD so nützlich macht, da sie sowohl eine Diagnose als auch eine Möglichkeit zur Behebung von Singularitäten bereitstellt. Falls eine Element von W sehr klein im Vergleich zu dem größten Element ist, entspricht dies einer Singularität (zumindest im Sinne numerischer Berechnungen, wo die Division durch extrem kleine Zahlen praktisch genauso schlecht ist wie eine Division durch 0). Eine Singularität bedeutet, daß das zugrundeliegende lineare Gleichungssystem unterdeterminiert ist. Dieses wiederum heißt, daß man willkürlich eine bestimmte Lösung aus dem unendlichen Lösungsraum herausgreifen kann. Insbesondere kann man W-1ii = 0 für sehr kleine Wii wählen - eine der seltenen Gelegenheiten, wo das Ergebnis einer Division durch (fast) 0 am besten gleich 0 anstelle von Unendlich zu setzen ist. Durch diese Editierung der isolierte Singulärwerte wird eine Dimensions-Reduktion des zugrundeliegenden Problems erreicht. Sie steht an zentraler Stelle bei allen auf SVD basierenden Routinen, wie z.B. MF_SVsolve,   MF_safeSolve) und VF_linfit. Die Schwelle für die Dimensions-Reduktion wird durch MF_SVDsetEdit gesetzt.

Die Singulärwerte sind in W in willkürlicher Reihenfolge gespeichert. Für manche Anwendungen (manuelle Inspektion der Singulärwerte, Dimensions-Reduktion usw.) ist es wünschenswert, die Singulärwerte in absteigender Reihenfolge zu sortieren. Hierbei müssen die Matrizen MU und MV entsprechend mit um-geordnet werden. Diese Aufgabe übernimmt die Funktion MF_SVsort.

Wenn Sie mittels MF_SVdecompose berechnete Ergebnisse mit SVD-Routinen anderer Anbieter vergleichen, beachten Sie bitte folgende Punkte:

  1. MF_SVdecompose nimmt selbst keine Dimensions-Reduktion (also keinen Ersatz sehr kleiner Singulärwerte durch 0) vor, sondern überläßt diesen Schritt nachfolgenden Funktionen oder ggf. manueller Inspektion durch den Anwender.
  2. Singulärwerte können vertauscht werden, wenn gleichzeitig auch die entsprechenden Spalten sowohl von MU als auch von MV vertauscht werden. Gelegentlich wird diese Eigenschaft genutzt, um die Singulärwerte in eine bestimmte (meist fallende) Reihenfolge zu bringen. Da der Vorzug geordneter Singulärwerte aber durch erhöhten Rechenaufwand erkauft wird, verzichtet MF_SVdecompose selbst auf diese Ordnung. Nur, wenn wirklich benötigt, sollte daher nach MF_SVdecompose noch die eben erwähnte Funktion MF_SVsort aufgerufen werden.
  3. Das Vorzeichen eines Singulärwertes kann vertauscht werden, wenn gleichzeitig die entsprechende Spalte von MV mit -1 multipliziert wird. MF_SVdecompose nutzt diese Eigenschaft, um alle nicht-verschwindenden Singulärwerte positiv zu machen.
SVD überprüfenDie folgende Routine demonstriert SVD und Überprüfung ihres Ergebnisses.

#include <VDstd.h>
#include <MDstd.h>
#include <conio.h>
void SVDtest( void )
{
    unsigned ht=386, len=52,            // beliebige Werte
             maxhtlen, sizeA, sizeV;
    dMatrix MA, MU, MV, MOne, MA2, MUV2, MDiffA, MDiffOne;
    dVector W;
    double err;
 
    maxhtlen = (ht >= len ? ht : len);
    sizeA = ht*len;
    sizeV = len*len;
 
    MA = MD_matrix( ht, len );   // Eingabe-Matrix
    MU = MD_matrix( maxhtlen, len );
    MV = MD_matrix( len, len );
    W = VD_vector( len );
    MOne = MD_matrix( len, len );   // Einheits-Matrix zum Vergleich
    MD_equ1( MOne, len );
 
    MA2 = MD_matrix( maxhtlen, len );
      // wird MU*W*MVT erhalten zum Vergleich mit der Eingabe-Matrix
      // führende Dimension maxhtlen ist für unter-determinierte Matrizen nötig
    MUV2 = MD_matrix( len, len );
      // wird MUT*MU, MVT*MV, MV*MVT erhalten zum Vergleich mit MOne
    MDiffA = MD_matrix( ht, len );  // für Zwischenergebnisse
    MDiffOne = MD_matrix( len, len );  // für Zwischenergebnisse
 
    MD_random( MA, ht, len, 1, -1., +1. );
      // MA mit Zufallszahlen zwischen -1 und +1 füllen
      // alternativ auf irgendeinem anderen Weg eine Test-Matrix erzeugen oder einlesen
    MD_SVdecompose( MU, MV, W, MA, ht, len );
 
      /* prüfen, ob MA = MU*W*MVT; Auswertung von rechts nach links: */
    MDdia_mulMT( MUV2, W, MV, len, len ); // MUV2 speichert W * MVT
    MD_mulM( MA2, MU, MUV2, maxhtlen, len, len );
    MD_subM( MDiffA, MA, MA2, ht, len );
      // für unter-determinierte Matrizen werden die Zeilen von ht bis maxhtlen-1 ignoriert
    err = VD_absmax( MDiffA[0], sizeA );
    printf( "SVDtest: Test für die SVD-Routine von OptiVec\n\n"
            "In den folgenden Tests in Double-Genauigkeit\n"
            "sind Fehler bis zu einigen Malen 1.e-14 in Ordnung\n"
            "bei Matrizen der Größenordnung 100x100 bis 1000x1000 Elemente:\n\n"
            "max. Fehler der Matrix-Rekonstruktion MA = MU*W*MVT: %lg", err );
 
    /* prüfe, ob MU orthogonal, also MUT*MU = (1): */
    MD_TmulM( MUV2, MU, MU, maxhtlen, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. Fehler von MUT*MU = (1): %lg", err );
 
    /* prüfe, ob MV orthogonal, also MVT*MV = (1): */
    MD_TmulM( MUV2, MV, MV, len, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. Fehler von MVT*MV = (1): %lg", err );
 
    /* prüfe, ob MV auch orthonormal, also MV*MVT = (1): */
    MD_mulMT( MUV2, MV, MV, len, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. Fehler von MV*MVT = (1): %lg", err );
    printf( "\n\nZum Beenden des Tests beliebige Taste drücken!" ); _getch();
 
    M_nfree( 8, MA, MU, MV, MOne, MA2, MUV2, MDiffA, MDiffOne );
    V_free( W );
}

QuerverweisMF_SVsort,   MF_SVsolve,   MF_safeSolve,   MF_solveBySVD,   Kap. 10

MatrixLib Inhaltsverzeichnis  OptiVec Home