2010-12-12 5 views
2

Ich berechne eine LDA-Transformation (lineare Diskriminanzanalyse) für eine Anwendung, an der ich arbeite, und ich habe diese notes (ab Seite 36, insbesondere Folie 47 in grün) verfolgt).gute numerische Lösung für LDA-Transformation

Ich tue dies in Python (mit numpy und scipy), und das ist, was ich mit gekommen sind:

import numpy as np 
from scipy.linalg import sqrtm 
... 
sw_inv_sqrt = np.linalg.inv(sqrtm(self.sigma_within)) 
self.d, self.v = np.linalg.eig(
    np.dot(
     np.dot(sw_inv_sqrt, self.sigma_between), 
     sw_inv_sqrt 
     )) 
self.v = np.dot(sw_inv_sqrt, self.v) 

Ich weiß, dass diese Implementierung korrekt ist, wie ich es im Vergleich zu anderen haben. Meine Sorge ist, ob dies gut Lösung in numerischen Sinn ist. Wenn ich meine Lösung mit anderen vergleiche, stimmen sie nur mit etwa 6 Dezimalstellen überein. Gibt es einen besseren Weg, dies numerisch zu tun?

+0

Dies ist im Wesentlichen eine mathematische Frage - Sie fragen nach einer numerischen * Methode *, nicht nach ihrer * Implementierung *. Sie könnten besser auf http://math.stackexchange.com/ fragen. –

+0

Benötigen Sie tatsächlich mehr als 6 Dezimalstellen? (Ich würde das eigentlich ziemlich gut finden). Sind die anderen Lösungen, auf die Sie verweisen, auch in Python implementiert? Warum vermuten Sie, dass die Lösungen der anderen numerisch besser sind als Ihre (soweit ich mich erinnere, ist dies die Standardmethode zur Berechnung des Projektionsvektors von Fishers linearer Diskriminanzanalyse (http://en.wikipedia.org/wiki/Linear_discriminant_analysis) "Sie könnten auch auf stats.stackexchange.com fragen. –

+0

@Sven: Ich stimme nicht unbedingt zu. Implementierungsdetails können * Roundoff-Fehler * einführen, der sich von * Trunkierungsfehler * unterscheidet, der im Wesentlichen mathematisch ist. Darüber hinaus mit allen Den Respekt, den ich den Leuten von math.stackexchange schulde, denke ich, dass es hier bessere Zahlen gibt als dort. –

Antwort

2

Versuchen Sie eigh anstelle von eig. Da Sigma^{- 1/2} B_0 Sigma^{1/2} symmetrisch ist, verwende eine angepasste Routine.

Achten Sie auch auf die Verwendung des richtigen Algorithmus bei der Berechnung von B_0. Für einen einfacheren Fall (den Sie hier anpassen können) siehe http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass_algorithm.

+0

Danke, das ist ein guter Tipp, aber es hat mein Ergebnis nicht viel verändert ... Könnte sein, dass Andre Recht hat , und dass verschiedene Plattformen die Unterschiede verursachen – Anonymous

+0

In LAPACK (mindestens Intel MKL-Version) gibt es eine Möglichkeit, die Eigenvektoren einer symmetrischen Matrix mit hoher Genauigkeit zu berechnen.Die Methode ist weit entfernt von t obwohl rivial. Wenn Sie Ihren Code in C oder Fortran neu schreiben möchten, können Sie davon Gebrauch machen. –