Forum: Offtopic Überbestimmtes lineares GLS in C++ lösen


von Vorname N. (zoltan)


Lesenswert?

Hallo,

ich habe ein Gleichungssystem deiser Form mit mehreren 1000 Zielen:

4 = 3.1 - 4 + 6.7 + 3 + 5 - 7.7 + 2 + 5 + 1
1 = 4.2 + 2 + 2.7 - 1 + 2 + 7.8 + 3 + 2 + 1
.
.
3 = 5.5 + 3 + 2 - 3 + 2 + 5 + 2.1 + 7.1 + 1


Matlab löst dieses überbestimmtes Gleichungssystem mit dem einfachen
Befehl x = H\z; wobei "H" zum Beispiel ein Matrix mit der Grösse
5945x9 und "z" ein Vektor 5945x1 ist.

Ich würde gern diese Rechnung in C++ Programmieren, um von Matlab
unabhängig zu sein. Wie kann ich dies am besten machen?

Gruss Zoltan

von Chris (Gast)


Lesenswert?

Ein klassischer Fall für do_what_i_mean().



Im Ernst: Wo liegt das Problem? Überleg dir einen Algorithmus, mit dem
man überbestimmte Gleichungsssysteme lösen kann. Um diesen Schritt
kommst du bei keiner Programmiersprache (außer Matlab) herum. Diesen
Algorithmus programmierst du dann in C++, was das geringste Problem
darstellen sollte.

Bei konkreten Problemen wird dir hier gerne geholfen, etwas
Eigeninitiative solltest du aber schon zeigen.

von Unbekannter (Gast)


Lesenswert?

Wenn Gauß dies hätte lesen müssen, hätte er sicher gekotzt!

von Thomas Winkel (Gast)


Lesenswert?

Hallo Zoltan,
deine Frage ist absolut berechtigt und konkret gestellt!
Warum sich hier immer mal wieder Leute "auskotzen" müssen, die nicht
in der Lage sind etwas zum Thema beizutragen ist mir schleierhaft.
Solche Probleme tauchen in der Praxis auf, wenn man Gleichungssysteme
lösen will, bei denen jede Gleichung mit einer gewissen Unsicherheit
belastet ist, wie es beispielsweise bei physikalischen Messdaten der
Fall ist, die mit Rauschen und Messfehlern überlagert sind. Deswegen
macht es da oft Sinn, deutlich mehr Gleichungen als Unbekannte
aufzustellen. Das hat den Vorteil, dass man auch dann noch brauchbare
Ergebnisse bekommt, wenn einige Gleichungen völlig daneben liegen.
Und nun die konkrete Lösung:
x=(H'*H)^(-1)*H'*z
wobei H' die Transponierte ist. Diese Lösung ist optimal, weil x so
bestimmt wird, dass der Fehler minimal ist.
Den C++ Code solltest du damit wol selber hinbekommen.

Viel Erfolg,
Thomas

von Vorname N. (zoltan)


Lesenswert?

x=(H'*H)^(-1)*H'*z wäre so etwas wie die Pseudionverse, nicht wahr?
Also
x=H_PsInv*z, wobei in meinem Fall H_PsInv eine 9x9 Matrix wäre.

Zum Lösen dieser quadratischen Gleichung würde ich das "Gauss
Verfahren" oder die "Carmensche Regel" anwenden. Dazu hätte ich noch
eine Frage: Bietet ein iterativer Verfahren wie Jakobi oder
Gauss-Seiderverfahren irgend welche Vorteile beim Numerischen Lösen von
Gleichungssystemen?

Gruss Zoltan

von Vorname N. (zoltan)


Lesenswert?

Kleine Korrektur: Das sollte natürlich "quadratischen Matrix" und
nicht "quadratischen Gleichung" heißen.

von Thomas Winkel (Gast)


Lesenswert?

Ja, (H'*H)^(-1)*H' ist die Pseudoinverse, eine richtige Inverse lässt
sich ja nur von quadratischen Matrizen bilden. Diese hat allerdings die
Dimension H', ist also nicht quadratisch.
Von den ganzen Verfahren und Regeln, von denen du da sprichst hab ich
keine Ahnung, einige lassen sich anscheinend eh nur auf quadratische
Matrizen anwenden. Ich würde das jetzt einfach so programmieren, wie es
da steht. Die Transponierte und die Multiplikation dürften sich einfach
mit je 2 ineinander verschachtelten for Schleifen reallisieren lassen.
Bei der Inversen müsste ich jetzt auch nachschauen, wie das nochmal
geht, dazu muss man erst die Determinante bilden. Warscheinlich gibts
aber auch freie Matrix Algebra Bibliotheken, kannst ja mal in
PC-Programmierung nachfragen, würde mich auch interessieren.

Gruß,
Thomas

von Vorname N. (zoltan)


Lesenswert?

Ich habe einen Fehler gemacht:
x = (H'*H)^(-1)*H'*z kann man umforem nach
H'*H*x = H'*z
H'*H ist meine 9x9 Matrix und H'*z ist ein 9x1 Vektor.
Ich habe dazu eine super Seite gefunden:
http://ceee.rice.edu/Books/CS/chapter3/data34.html

Danke für die Hilfe Thomas,
Gruss Zoltan

von Detlef A (Gast)


Lesenswert?

Gute Quelle ist auch 'Numerical Recipes in C', googeln nach Titel
bringt nen link der Cornell Uni., der einzelne Kapitel sehen läßt.

Cheers
Detlef

von Marco S (Gast)


Lesenswert?

Hallo Zoltan.

Hatte auch mal x = H \ z zu lösen. Als ich spitz gekriegt habe, dass
MATLAB auch nur LAPACK-Routinen & co aufruft, habe ich von LAPACK die
Routine DGELS verwendet. Diese ist in Fortran programmiert, lässt sich
aber auch von C aus verwenden. Vorteil: Es funktioniert (mature-code).

Schon mein Informatik-Prof sagte damals: "Bevor Sie sich in C++
bezüglich der Lösung mathematischer Probleme verstricken, kaufen Sie
sich einfach eine fertige und erprobte Lösung; z.B. Numerical
Receipes"

Gruß
Marco
-

von Vorname N. (zoltan)


Angehängte Dateien:

Lesenswert?

Im Anhang befindet sich eine Datei mit zwei Versionen (c++ / Matlab) der
Householdertransformation zur Lösung von überbestimmten linearen
Systemen.

Gruss
Zoltan

von Daniel Hoberg (Gast)


Lesenswert?

Noch ein Hinweis zu der ganzen Diskussion:

Mit dem "Template Numerical Toolkit (TNT)" gibt es ein (freies) auf
Templates basierendes, einfach zu benutzendes Numerik-Paket für C++.
Lässt insgesamt ein wenig MATLAB-Flair aufkommen... ;)

Zu finden unter: http://math.nist.gov/tnt/

Ich hab's benutzt und find es wesentlich angenehmer als LAPACK. Zudem
ist es (für meine Implmentierungen) um ca. 25% schneller.

Grüße, Daniel

P.S.: die Überführung eines überbestimmten GLS der Form   H*x=z   in
die Form   H'*H*x = H'*z   nennt sich 'Gauss-Transformation'. (Wie
Thomas richtig sagte, liefert das eine exakte Lösung, im Sinne von
least squares.)
Die Koeffizientenmatrix H'*H ist dabei quadratisch. Das GLS lässt sich
dann mit einer LU- oder QR-Zerlegung, wie sie im TNT zu finden ist, ganz
'einfach' lösen...

von Detlef A (Gast)


Lesenswert?

H'*H ist quadratisch und symmetrisch und läßt sich Cholesky-zerlegen,
das ist numerisch günstiger.

von Michael (Gast)


Lesenswert?

Mal ne dumme Frage: Was gibt es denn hier zu lösen? Ich sehe gar keine
Unbekannten ?!?

von Jürgen Schuhmacher (Gast)


Lesenswert?

Aber Michael, die Nichterkennbarkeit von Variablen ist doch DAS
Wensensmerkmal von Unbekannten. Deshalb heissen sie ja so.

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.