Lapack++: Hat LaLinearSolve Probleme mit fast singulären LGS?

Kiddycat

Hallo Darkfate,

darkfate schrieb:

In solchen Szenarien wird wohl jeder Algorithmus Schwierigkeiten haben.

Ja, vermutlich hast du Recht.

darkfate schrieb:

Hast du versucht dein Gleichungssystem mit einem Vorkonditionierer zu behandeln?
Hast du es schonmal versucht mit MATLAB/Mathematica/Octave zu lösen?

Ich muss zugeben, dass ich auf die Idee einen Vorkonditionierer zu benutzen noch nicht gekommen bin. Aber ich werde mal suchen, ob ich was brauchbares finde. Das LGS wird von Matlab gelöst.

darkfate schrieb:

Immer das selbe LGS?

Nein. Wie gesagt, A und g werden in jedem Schleifendurchgang neu berechnet. Und A und g hängen ja auch von meinem Datensatz ab. Bei manchen Datensätzen wird das LGS gelöst, bei anderen Datensätzen rechnet das Programm eine Weile (das LGS wird gelöst) und irgendwann bricht es ab.

darkfate schrieb:

Poste doch mal den Code.

Der Code ist lang. Was brauchst du genau?
An sich läuft das ganze so ab:

INPUT: LaGenMatDouble* x
OUTPUT: LaGenMatDouble* n_new
BEGIN LevenbergMarquardt: 
	Initialisiere alle Variablen (Zwischenwerte, Schrittweitensteuerung, Dämpfer, etc.)
	double f = functionValue(*x)
	LaGenMatDouble* J = calcJacobi(*x)

	LaGenMatDouble* A = new LaGenMatDouble(J->rows(),J->rows());
	Blas_Mat_Mat_Mult(*J,*J,*A,false,true,1.0,0.0);
	LaGenMatDouble* g = new LaGenMatDouble(*J);
	*g *= -f;

	Überprüfe, ob Abbruchbedingung schon erfüllt ist (Ob schon nah genug an Nullstelle...)

	LaGenMatDouble* h = new LaGenMatDouble(A->rows(),g->cols());
	LaGenMatDouble* x_new = new LaGenMatDouble(x);
	LaGenMatDouble A_mu_I;

BEGIN WHILE (Abbruchbedingung nicht erfüllt und Höchstanzahl Iterationen nicht überschritten):

	A_mu_I = LaGenMatDouble::eye(5);
	A_mu_I *= mu;
	A_mu_I = *A - A_mu_I;

	LaLinearSolve(A_mu_I,*h,*g);

BEGIN IF Abbruchbedigung erfüllt: Return x
BEGIN ELSE

	x = x+*h;
	J = calcJacobi(&x);	
	g = J;
	f = functionValue(&x);
	*g *= -f;
	Blas_Mat_Mat_Mult(*J,*J,*A,false,true,1.0,0.0);

	Überprüfe Abbruchbedingung
	Aktualisiere mu

END ELSE
END WHILE
END LevenbergMarquardt

Kiddycat

Kiddycat schrieb:

f: IR^5->IR ist eine Funktion.
J_f die Jacobimatrix von f.
double mu

A = J'*J (J transponiert mal J, also eine symmetrische 5x5-Matrix)
g = -J'*f ( - J transponiert mal f, also eine 5x1-Matrix)
I = LaGenMatDouble::eye(5)

finde h, sodass
(A-mu*I)*h = g

Verrat mal mehr über Dein Problem. Dass f von 5->1 Dimensionen abbildet, finde ich komisch. Das entspricht 5 Freiheitsgraden aber nur einem Messpunkt. Dadurch wird A singulär. Das GLS, das Du aufstellst, ist auch falsch. Es muss heißen (A**+**mu*I)*h = g -- jedenfalls, wenn mu>0 gilt. Mit dem mu wird das System "regularisiert". Das heißt, dass A singulär sein darf. A+mu*I mit mu>0 hat vollen Rang (ist nicht singulär). Das ist ja der Witz an dem Verfahren. Es vereinigt die besten Eigenschaften des Gauss-Newton-Verfahrens und des Gradientenabstiegs.

Kiddycat

krümelkacker schrieb:

Das GLS, das Du aufstellst, ist auch falsch. Es muss heißen (A+mu*I)*h = g -- jedenfalls, wenn mu>0 gilt.

Wow, danke! Da hab ich was falsch gelesen!

Über die Dimensionen denk ich noch mal nach.
Das Programm läuft, wenn ich die diagonale Vorkonditionierung benutze und B=D^{-1} wähle. Das LGS ABB(B^{-1}B{-1}x)=g wird gelöst...
Das finde ich jetzt aber nicht so zufrieden stellend. Wer sagt mir denn, dass zweimal B dran multiplizieren immer ausreicht?

Ich mache mir jetzt mal Gedanken darüber, ob meine Funktion so stimmen kann.
Ich habe 5 Freiheitsgrade und mindestens 6 Messpunkte, es können aber auch wesentlich mehr werden. Ich will die Summe der quadratischen Abstände der Punkte zu einer Oberfläche, die durch die 5 Freiheitsgrade beschrieben wird, minimieren.

Die Messpunkte gehen in die Funktionsauswertung und die Jacobimatrix mit ein.

d(S(x_1,x_2,x_3,x_4,x_5),p_i) ist der Abstand vom Punkt p zur Oberfläche S.
Wenn P={p_1,...,p_n} die Menge der Punkte ist und x=(x_1,...,x_n) der Vektor, dann ist meine Funktion
f(x) = sum_{i=1}^n{d(S(x),p_i)}

Ich denke, Vorkonditionierung kannst Du hier ignorieren -- es sei denn, Du willst das LGS iterativ lösen (mit dem CG-Algorithmus zB). A+mu*I ist symmetrisch positiv definit, kann also bei 5 Freiheitsgraden elegant über eine Cholesky-Zerlegung gelöst werden.

Kiddycat schrieb:

dann ist meine Funktion
f(x) = sum_{i=1}^n{d(S(x),p_i)}

Das sieht jedenfalls ungeschickt aus. Es geht wahrscheinlich viel besser. Was ist S? Wie rechnest Du d aus?

Ich hätte jetzt mit einem f der Art
f : R^5 --> R^(6*3)
gerechnet. Fünf Freiheitsgrade werden auf 6 dreidimensionale Fehlervektoren abgebildet. Mit dem LM-Algorithmus kannst Du f(x)^T f(x) = ||f(x)||^2 minimieren, was wirklich die Summe der quadratischen Abstände ist. Und wenn J^T f(x) verschwindet (Fehlervektor ist 0 oder steht senkrecht auf den durch J aufgespannten Unterraum) bist Du fertig. So, wie Du es gemacht hast, verschenkst Du Informationen.

darkfate

Kiddycat schrieb:

Ich will die Summe der quadratischen Abstände der Punkte zu einer Oberfläche, die durch die 5 Freiheitsgrade beschrieben wird, minimieren.

Orthogonalregression? Meinst du die lotrechten Abstände?

Kiddycat schrieb:

d(S(x_1,x_2,x_3,x_4,x_5),p_i) ist der Abstand vom Punkt p zur Oberfläche S.
Wenn P={p_1,...,p_n} die Menge der Punkte ist und x=(x_1,...,x_n) der Vektor, dann ist meine Funktion
f(x) = sum_{i=1}^n{d(S(x),p_i)}

Wieso nutzt du nicht das strenge Gauß-Helmert Modell (GHM), den Allgemeinfall der Ausgleichungsrechung?
Es ist schlechthin Standardverfahren um Formen wie Quadriken auszugleichen, ist in der DIN normiert und liefert immer das exakte Minimum der Lagrange Funktion.

Bei vielen anderen Verfahren hat sich in der letzten Zeit herausgestellt, dass sie nicht die selbe Funktion minimieren wie das strenge Gauß-Helmert Modell und somit je nach Form nur in Ausnahmefällen ein qualitativ identisch gutes Ergebnis liefern.

Zudem kannst du jedes GHM in ein qualitativ gleichwertiges Gauß-Markov Modell überführen womit du ein weiteres Gütesiegel für deine Ergebnisse hättest.

Kiddycat

krümelkacker schrieb:

Ich hätte jetzt mit einem f der Art
f : R^5 --> R^(6*3)
gerechnet. Fünf Freiheitsgrade werden auf 6 dreidimensionale Fehlervektoren abgebildet. Mit dem LM-Algorithmus kannst Du f(x)^T f(x) = ||f(x)||^2 minimieren, was wirklich die Summe der quadratischen Abstände ist. Und wenn J^T f(x) verschwindet (Fehlervektor ist 0 oder steht senkrecht auf den durch J aufgespannten Unterraum) bist Du fertig. So, wie Du es gemacht hast, verschenkst Du Informationen.

Hm. Aber wieso IR^{6x3}?
Angenommen ich hätte f \in IR^{6x3} und angenommen, ich wollte f'f \in IR^{3x3} minimieren, dann hätte ich eine Jacobi"matrix" der Dimension 3x3x5. Das sieht dann aber gar nicht mehr so geschickt aus...

Aber ich kann f:IR^5->IRk (k=Anzahl der Punkte) wählen.
Nur, dann ist J kx5-dimensional, A=J'J 5x5 dimensional und g=J'f 5x1 dimensional. Was hat sich dann großartig geändert?
Ich kann auch f:IR^5->IR{(6*3)x1} wählen. Ich weiß nicht, ob du das nicht gemeint hast. Aber dann ist J'f wieder 5x1 und J'J wieder 5x5...

darkfate schrieb:

Wieso nutzt du nicht das strenge Gauß-Helmert Modell (GHM), den Allgemeinfall der Ausgleichungsrechung?
Es ist schlechthin Standardverfahren um Formen wie Quadriken auszugleichen, ist in der DIN normiert und liefert immer das exakte Minimum der Lagrange Funktion.

Hast du da mal einen Link zu? Ich finde leider keine Erklärung.

Ein paar Fragen, die sich mir stellen sind:
Wie verhält sich Gauß-Helmert bei großen Datensätzen? (Laufzeit/Speicherbedarf)
Du sagst, es sei für Quadriken Standardverfahren, was ist, wenn man andere Formen hat?
Gibt es einen Beweis, dass es immer das exakte Minimum liefert?
Brauche ich Startwerte?

darkfate

Kiddycat schrieb:

Hast du da mal einen Link zu? Ich finde leider keine Erklärung.

Das strenge Gauß-Helmert Modell ist der allgemeinste Fall der Ausgleichungsrechnung.

Wenn man es ganz allgemein Lesen will:
[1] http://www.archive.org/details/ausgleichungsrec00helmrich
Würde ich jetzt persönlich nicht empfehlen, da es bereits wesentlich kompakter und mit neueren mathematischen Notationen zusammengefasst wurde.

Als Referenz würde ich dir folgende Links empfehlen:
Lenzmann & Lenzmann: Allgemeine Notation + Beispiel (Parabel)
[2] http://www.herbert-wichmann-verlag.de/ai/resources/756dafe4425.pdf
Prof. Neitzel: Beispiel (Geraden).
[3] http://edoc.gfz-potsdam.de/gfz/get/11374/0/13a7fdd69e092576ee682dd36e29cd83/11374.pdf
Beispielimplementierung der Ebene in MATLAB (Hinweis auf Unterschiede zwischen genähertem und strengen Modell):
[4] http://forum.derletztekick.com/index.php?id=1005
[5] Meine Theisis (bei Bedarf PN)
Nutze ausschließlich diese Quellen, weil nahezu 90% der Bücher ohne Hinweis das genäherte und nicht das strenge GHM Modell benutzen. Es unterscheidet sich je nach Funktion und rauschen sehr stark in der Qualität der Ergebnisse. Solange man sich nicht intensiv mit diesem Thema befasst hat, wird man es sehr schwer haben das strenge und das genäherte Modell auseinanderzuhalten.

Kiddycat schrieb:

Wie verhält sich Gauß-Helmert bei großen Datensätzen? (Laufzeit/Speicherbedarf)

Ich habe das GHM anhand von ~40.000 (X,Y,Z) Tripeln (Laserscannerdaten) untersucht. Meine Zielfunktion war: n_xX+n_yY+n_zZ+d (Ebene in der HNF). Für die komplette Ausgleichung bis zur Ergebnisausgabe benötigte MATLAB 100ms. Ist es zu langsam, so gibt es mehrere Möglichkeiten dem Vorzubeugen. Man sollte allerdings das Gleichungssystem so aufstellen, dass alle Matrizenelemente (im allgemeinen) möglichst gleich groß werden und das Gleichungssystem nicht mit steigender Anzahl der Punkte wächst. Stichwort: Bedingungsgleichungen in die Matrix einführen.

Kiddycat schrieb:

Du sagst, es sei für Quadriken Standardverfahren, was ist, wenn man andere Formen hat?

Laut DIN werden lineare sowie nicht lineare Zielfunktionen akzeptiert. Eine Einschränkung bezüglich der Verwendung wird daher nicht vorgenommen.

Kiddycat schrieb:

Gibt es einen Beweis, dass es immer das exakte Minimum liefert?

Link [2] Seite 1.

Kiddycat schrieb:

Brauche ich Startwerte?

Wenn du linearisieren musst: Ja
Richtig aufgestellt konvergiert es sehr gut in meinem Fall sogar mit willkürlichen Näherungswerten.

Kiddycat schrieb:

Ich kann auch f:IR^5->IR{(6*3)x1} wählen. Ich weiß nicht, ob du das nicht gemeint hast. Aber dann ist J'f wieder 5x1 und J'J wieder 5x5...

Ich meinte jedenfalls mit R^(6*3) keine 6x3 Matrix, sondern einen 18-dim-Vektor, der sich ergibt, wenn man 6 3-dimensionale "Fehlervektoren" "untereinanderhängt". Wenn Du f so aufstellst, ist q(x) := 1/2 f(x)^T f(x) proportional zur Summe der quadrierten Abstände und die Jacobi-Matrix J von f ist eine 18x5 Matrix. Weiter ist J^T J die richtige Hessenmatrix von q, gibt also die Krümmung an von q in allen Richtungen an.

Bei Deinem Ansatz ist f ein Skalar, J eine 1x5 Matrix und J^T J eine singuläre 5x5 Matrix, welche eine sehr schlechte Näherung der Hessenmatrix der eigentich zu minimierenden Funktion ist. Das, was dabei rauskommt, ist nichts weiter als ein kompliziert implementierter Gradientenabstieg.

Kiddycat

Hab meinen Fehler verstanden. Danke für die Hilfe!

otze

Ehhh hab ich das richtig verstanden, das wir hier einfach über lineare Regression reden? Wozu dann so ein kompliziertes Verfahren?

http://www.neuroinformatik.ruhr-uni-bochum.de/PEOPLE/igel/LT/LinearRegressionAndDiscrimination.pdf Folie 13. Gilt natürlich nur, wenn sichd er Threadstarter dazu entscheidet anzunehmen, dass seine Fehler Normalverteilt sind und er deswegen den Mean-Squared-Error direkt ansetzen kann.

"...Oberfläche, die durch die 5 Freiheitsgrade..."

klingt nicht nach einer Geraden oder einer Ebene, sondern nach etwas potentiell gekrümmten. Und je nachdem, was der OP genau machen will, wird das resultierende Problem beliebig kompliziert/nichtlinear.

otze

ahh ich sehs. jemand anders hat die Ebene ins Spiel gebracht, und die blieb dann haften. Naja, wenns beliebig kompliziert werden kann, sind Verfahren die linear Annähern beliebig schlecht :).

darkfate

otze schrieb:

Ehhh hab ich das richtig verstanden, das wir hier einfach über lineare Regression reden? Wozu dann so ein kompliziertes Verfahren?

http://www.neuroinformatik.ruhr-uni-bochum.de/PEOPLE/igel/LT/LinearRegressionAndDiscrimination.pdf Folie 13. Gilt natürlich nur, wenn sichd er Threadstarter dazu entscheidet anzunehmen, dass seine Fehler Normalverteilt sind und er deswegen den Mean-Squared-Error direkt ansetzen kann.

Hast du dazu irgendwas implementiertes oder rechne mal bitte folgende Werte für eine Gerade durch:

X Y
0 0
1 1
2 4
3 9

Zielfunktion: y = a*x + b

otze

nö, hab da nichts implementiert, und für das Problem kann mans ja noch analytisch lösen (sofern man sich nicht verrechnet, wobei ich mir fast sicher bin...ich kenn mich!): 5*x-4

Das Verfahren gibt aber, sofern man es richtig implementiert (d.h. ohne Moore-Penrose Inverse als letzte "Vereinfachung"...) das selbe, numerisch stabile, Ergebnis.

darkfate

otze schrieb:

nö, hab da nichts implementiert, und für das Problem kann mans ja noch analytisch lösen (sofern man sich nicht verrechnet, wobei ich mir fast sicher bin...ich kenn mich!): 5*x-4

Das Verfahren gibt aber, sofern man es richtig implementiert (d.h. ohne Moore-Penrose Inverse als letzte "Vereinfachung"...) das selbe, numerisch stabile, Ergebnis.

Es geht mir nicht um irgend ein Ergebnis.
Folie 12 weist ein Minima nur in eine Achse (y_i) nach.
Die Achse x_i wird überhaupt nicht berücksichtigt.
Ist ein wenig verwirrend weil aus Folie 6 hervorgehen sollte
dass beide Achsen eine Verbesserung erfahren.

Das ist ein schlichtes Gauß-Markov Modell ohne
Berücksichtigung der Varianzen und Kovarianzen.
Stichwort: Satz von Gauß-Markov.

Damit kommt man jedoch nicht zur kleinsten Summe der Verbesserungen.

otze

darkfate schrieb:

Es geht mir nicht um irgend ein Ergebnis.
Folie 12 weist ein Minima nur in eine Achse (y_i) nach.
Die Achse x_i wird überhaupt nicht berücksichtigt.

natürlich nicht. Weil das der Eingaberaum ist, über den das Ergebnis generalisiert wird. In dem Sinne ist "x" keine Zufallsvariable mehr, sondern nur y.

Folie 6 ist übrigens eine andere Darstellung. da gibts keine x und y Achse sondern nur x1 und x2. Die abgebildete Funktion die du dort siehst ist nur die gerade für f(x)=0.

Siehe desweiterne Folie 9 und 10.

//edit und weil ich selber in der Vorlesung dring esessen hab: Folie 6 war eine einfache Erklärung, was eine lineare Funktion ist, weil das in späteren Kapiteln (es geht fast nur um lineare binäre Klassifikation) relevant ist. von daher hat Folie 6 keine tiefere Aussagekraft,