Algorithmus gesucht

Freier Radikaler

Hallo,

ich suche einen Algorithmus für folgendes mathematische Problem:

Gegeben sind n Punkte in einem dreidimensionalen Raum deren Koordinaten bekannt sind. n >= 1 in der Regel aber >= 4.
Neben den Punkten sind Verhältniswerte zu einem Punkt Px zwischen allen diesen bekannten Punkten n bekannt. Zu berücksichtigen ist, dass diese Verhätniswerte nur näherungsweise genau sind.

Gesucht wird dieser Punkt P.

Noch ein konkretes Beispiels zum besseren Verständnis.

Man stelle sich 4 Punkte vor. P1 (1|1|1), P2 (2|2|2), P3 (3|3|3) und P4(4|4|4).
Die Verhältnisse zum gesuchten Punkt Px sind zb. 4 zu 3 im vergleich zu P1 und P2, 1 zu 2 im Vergleich zu P1 und P3 und 4 zu 8 im Vergleich zu P3 und P4 usw...
(Das ist jetzt ein ausgedachtes Beispiels, welches ggf. nicht besonders Sinn macht, aber das Problem beschreiben soll)

Der gesuchte Punkt Px ist nun jener Punkt, an dem der Fehler der Verhältniswerte zu allen ermittelten Referenzpunkten (P1-4) am geringsten ist.

Normalerweise müsste es sich dabei um ein iteratives Verfahren handeln. Die Frage ist nur: Welcher Algorithmus ist bei dem Problem am effizientesten mit Hinsicht auf die Komplexität.

Die entscheidende Frage ist, was genau bedeutet "Fehler am geringsten"?

Wenn sich Px im Verhältnis a:b zu Pi und Pj verhält, ist der erwartete Ort Pxij = a/(a+b)*Pi + b/(a+b)*Pj.
Mach das für alle gegebenen Verhältnisse und du erhältst eine Punkmenge von erwarteten Orten.
Davon kannst du z.B. das arithmetische Mittel nehmen (wenn du die Summe der Abstände minimieren möchtest).

Die Komplexität wäre hier ziemlich gering.

knivil

Ich verstehe die Beschreibung von Verhaeltnis zwischen Punkten nicht oder meinst du Strecken zwischen zwei Punkten. Druecke das Verhaeltnis doch mal etwas mathematischer aus und nicht nur als Beispiel! Dann beschreibst du, was Fehler genau bedeutet! Und dann waehle bitte ein Beispiel, be der nicht alle Punkte auf einer Geraden liegen mit Loesung fuer den gesuchten Punkt!

Freier Radikaler

Das Beispiel sollte nur die Dimensionalität verdeutlichen. Ich gebe zu, dass es nicht perfekt ist. In meinem Anwendungsfall sind die Referenzpunkte normalerweise linear unabhängig gewählt.

Mit "Fehler am geringsten" meine ich die Summe der Verhältnisabweichungen an dem ermittelten Punkt muss global betrachtet am geringsten sein. Beispiel:

Der aktuell untersuchte Punkt hat ein Verhältnis von 1,5:1 zu den Punkten P1 und P2. Gemessen wurde allerdings ein Verhältnis von 1,4:1. Dann ist der Fehler prozentual zur 1,0 gesehen = 0,1. Gleiches würde man dann noch mit allen anderen Verhältnissen zu den restlichen Punkten durchführen und die dort ermittelten Fehler zu der 0,1 hinzu addieren. Somit hätte man eine Fehlerwert am aktuellen Punkt. Per Brute force müsste man jetzt alle Punkte um die gegebenen Punkte herum abklappern um den Punkt mit dem so geringsten Fehler zu finden. Das Problem ist aber, dass Bruteforce im dreidimensionalen Raum O((n/g)³) hat. Wobei n die Kantenlänge des Bruteforce Würfels ist und g die Ganularität mit der gesucht wird. Beispiel:

Würfel von 3 Meter, Granularität von 0,1m: (3/0,1)³ = 27000 Brute-Force Schritte

Ich habe aber eben die starke Vermutung, dass das einfacher geht. Aber wie
Ich hab mir jetzt mal K-means angeschaut. Mein Problem ist glaube ich ein spezialfall dieses Algorithmuses bei dem Cluster = 1 und die Centroidfunktion im prinzip das ist was ich suche. Allerdings eben nicht mit konkreten Abständen, sondern lediglich mit Verhältnissen. Kann aber auch sein, dass ich da voll auf dem Holzweg mit bin...

Nochmal zum Verhältnis für knivil:

Mit Verhältnis zwischen zwei Punkten meine ich den Quotienten der Abstände. Das ist allerdings kein fester Abstand sondern wie gesagt ein Verhältnis. Konkretes Beispiel:

Px ist 3 mal so weit von P1 entfernt wie Px von P3 entfernt ist. Das ergibt ein Verhälntnis von 3:1. Mathematisch betrachtet ergibt sich dadurch ein Ellipsoid um den Punkt P3. Auf diesem Ellipsoid liegt der gesuchte Punkt Px. Da es sich dabei um einen Körper handelt, müssen zur Bestimmung von Px noch weitere Verhältnisse zu anderen bekannten Raumpunkte bekannt sein, die dann ihrerseits wieder einen Ellipsoid bilden. Betrachtete man nun alle diese so erzeugen Ellipsoiden, gibt es einen Punkt an dem die Entferung zu allen Ellipsoiden-Oberflächen am geringsten ist. Genau das ist der Punkt, den ich suche...

@Barycenter: Deine Formel kann nicht wirklich richtig sein, da es keinen konkreten Punkt sondern einen Elipsoiden bei der reinen Betrachtung von Verhältnissen gibt. Einen Punkt kann man nur berechnen, wenn man feste Entferungswerte zu Grunde legt. Das würde im zweidimensionalen mit zwei Referenzpunkten, was du aufgezeigt hast allerdings eine doppelte Lösung implizieren. Das ist nur ab drei linear unabhängigen Referenzpunkten eindeutig lösbar...
Bei Verhältnissen handelt es sich wie gesagt um eine Ellipse im Zweidimensionalen und einen Ellipsoiden im Dreidimensionalen. Wenn das Problem so einfach gelagert wäre, hätte ich kein Problem

krümelkacker

Freier Radikaler schrieb:

[...] Summe der Verhältnisabweichungen [minimieren] [...]

[...] einen Punkt an dem die Entferung zu allen Ellipsoiden-Oberflächen am geringsten ist. [...]

Diese Wünsche sind inkompatibel. Es sind zwei unterschiedliche Probleme mit jeweiligen Lösungen, die nicht übereinstimmen müssen. Und beide sind einigermaßen schwierig. Das erste ist auch sicher nicht das, was du willst.

Ohne eindeutig mathematisch formuliertes Problem gibt es auch keine Lösung. Dein Problem ist also (noch) nicht die Lösung eines Problems, sondern die Formulierung eines Problems. Welches Problem du lösen möchtest, hängt sicher auch davon ab, woher die Eingaben kommen und wie sie verfälscht sind.

Wenn du die Summe der Quadrate der Verhältnisabweichungen minimieren willst, dann ist das ein nicht-lineares Ausgleichsproblem. Da bietet sich das Levenberg-Marquardt-Verfahren an. Allerdings muss dir klar sein, dass wenn du Zähler und Nenner bei den Verhältnissen vertauschst, eine andere Lösung heraus kommt. Deswegen sprach ich auch vorhin die Herkunft der Eingaben an. Du willst doch sicher eine Lösung bestimmen, die möglichst viel mit der Realität zu tun hat, die möglichst wahrscheinlich ist. Stichwort "Maximum-Likelihood". Oft reduziert sich das (ggf mit Vereinfachungen) auf ein Ausgleichsproblem (least squares, also Minimierung der Summe von Quadraten von Abweichungen).

Klingt aber interessant. Wo tauchen denn solche Probleme auf?

Freier Radikaler

Stimmt. Gerade ist mir mein Modellfehler auch aufgefallen. Es handelt sich überhaupt nicht um eine Ellipse. Sondern um einen seltsamen Körper der aussieht wie ein Kegel mit abgerundetem boden und Spitze. Also nichts, was mann einfach so bestimmen könnte glaube ich.

Das was du meintest, schließt sich bei Ellipsen aus aber ich meinte eben keine Ellipsen s.o. Sondern eben diesen komischen Verhältniskörper dessen Oberfläche alle Punkte bilden, die exakt das gewünschte Verhältnis zu den beiden Punkten haben.

Mit der maximum likelyhood methode hast du Recht. Genau das will ich erreichen. Bei den Daten, die ich verarbeiten möchte handelt es sich um Signalfeldstärken von Funksignalen. Dabei habe ich ein Sensornetz von Empfangsstationen, deren genaue Position bekannt ist. Ein Sender innerhalb dieses Feldes sendet ein Signal aus, welches in Abhängigkeit der Entferung von allen Empfängern mit unterschiedlicher Empfangsstärke wahrgenommen wird. Ein beliebter Fehler ist, die dBM-Werte einfach 1:1 in Entferung umzurechnen und dann die Spheren, die um die Basisstationen so entstehen auf Schnittpunkte zu untersuchen. Allerdings vernachlässigt man dabei, dass das mapping dBm Wert <-> Entferung nicht konstant ist, sondern sich im Laufe der Zeit ändert. Aus diesem Grund muss man die Berechnung auf Basis der Signalstärkenquotienten vornehmen.
Der gesuchte Punkt ist dann der warscheinlichste Aufenthaltspunkt des Senders.

Ich hoffe das hat das Problem jetzt verständlich genug gemacht

Ich werd mir mal das Levenberg-Marquardt-Verfahren, welches du vorgeschlagen hast anschauen. Mal schaun was das taugt und ob ich es verstehen

Freier Radikaler

Hmm wenn ich mir das Verfahren bei Wikipedia so anschaue, verlangt das Teil eine Funktion. Die Funktion kenne ich aber gar nicht

knivil

Mit Verhältnis zwischen zwei Punkten meine ich den Quotienten der Abstände.

Na also ...

Mathematisch betrachtet ergibt sich dadurch ein Ellipsoid um den Punkt P3

Das halte ich fuer ein Geruecht, was macht den Punkt P3 im Vergeich zu P1 besonders? Warum nicht um P1?

dass das mapping dBm Wert <-> Entferung nicht konstant ist, sondern sich im Laufe der Zeit ändert.

Achso, du hast also eine Welle/Sender und moechtest anhand der Zeitreihe deiner Messstationen den Ursprung bestimmen. Jetzt ist natuerlich die Frage: Wie weit sind deine Empfaenger untereinader entfernt, d.h. spielt die Geschwindigkeit eine Rolle oder kann instantan angenommen werden? Wenn instantan, dann kannst du einfach eine quadratische Funktion in Abhaengigkeit von der Intensitaet zu einem fixen Zeitpunkt fitten. Funktioniert unter idealen Bedingungen ... im Vakuum also. Wie gut in anderen Medien kann durch Simulation mit dieser Modellanahme geprueft werden.

dBm Wert <-> Entferung nicht konstant ist

Konstant? Nein, die Intensitaet sollte quadratisch von der Entfernung abhaengen. Mit dBm als Einheit hatte ich aber noch nie was am Hut. Aber wenn sie mit groesserer Entfernung abnimmt reichts vielleicht trotzdem.

Die Funktion kenne ich aber gar nicht

Das ist deine Modellannahme.

PS: Dein Ansatz mit den Verhaeltnissen ist quark.

Freier Radikaler

Obala da hab ich mich verschrieben. Ich meinte nicht Quotienten der Abstände sondern Quotienten der RSSI-Werte (RSSI = Radio Signal Strength Indikation).

WIe du schon sagtest ist die quatratische abnahme der Ideale fall.

Das halte ich fuer ein Geruecht, was macht den Punkt P3 im Vergeich zu P1 besonders? Warum nicht um P1?

Ja Ellipsoid ist es nicht das habe ich oben schon geschrieben. Aber auf jeden Fall ein mathematisch definierbarer Körper, der sich um den Punkt herum aufzieht, dessen Verhältniszahl kleiner ist. Also beim verhältnis 1:2 um den Punkt mit Verhältniszahl 1.

Konstant? Nein, die Intensitaet sollte quadratisch von der Entfernung abhaengen. Mit dBm als Einheit hatte ich aber noch nie was am Hut. Aber wenn sie mit groesserer Entfernung abnimmt reichts vielleicht trotzdem.

Hast du das überhaupt richtig gelesen? Ich habe damit das Mapping von dBm wert auf einen bestimmte Entferungszahl gemeint. Mir ist auch klar, dass die dBm-Werte quatratisch mit der Entfernung abnehmen. Aber genau das gillt eben nur unter idealbedingungen, die allerdings so in der Wirklichkeit fast nicht auftreten. Dort gibt es massenweise Störeinflüsse.
Das Mapping an sich ist 1:1, da immer ein fester rssi wert einer bestimmten entfernung so zugeordnet wird. Das ist aber hochgradig Fehleranfällig. Deswegen der Umweg über die Verhältnisse!

Der Ansatz über die Verhältnisse ist nicht quark, da er in unterschiedlichsten Papers so vorgeschlagen und umgesetzt wurde.

krümelkacker

Freier Radikaler schrieb:

Mir ist auch klar, dass die dBm-Werte quatratisch mit der Entfernung abnehmen.

Ich denke, dass sie das nicht tun. Die dBm-Zahlen sind ja logarithmierte Leistungen.

Nichts destro trotz stellt sich die Frage, wie die Messfehler bei diesen dBm-Zahlen verteilt sind; denn daraus ergibt, wie du das Ausgleichsproblem aufstellen musst. In diesem Fall ist das natürlich schwer zu schätzen, zumal das ein sehr einfaches Modell ist, was wohl nicht viel mit der Realität zu tun haben wird (da hat man ja noch Reflektionen, Interferenzen, ...).

Aber ich würde mir schon ein paar Gedanken dazu machen, wie ich die "Gleichungen" des Systems aufschreibe und gewichte. Das macht bei einem Ausgleichsproblem schon einen Unterschied.

$\frac{\left|| P\_x - P\_1 \right||}{\left|| P\_x - P\_2 \right||} \approx \frac{4}{3}$

$\frac{\left|| P\_x - P\_2 \right||}{\left|| P\_x - P\_1 \right||} \approx \frac{3}{4}$

$\frac{\left|| P\_x - P\_1 \right||}{\left|| P\_x - P\_1 \right|| + \left|| P\_x - P\_2 \right||} \approx \frac{4}{3+4}$

werden dir unterschiedliche Ergebnisse liefern, wenn alle Gleichungen nicht exakt erfüllbar sind. Gewichten könnte man so eine Gleichung auch noch per pow(10.0,(max(dBm1,dBm2)/10.0), um nur mal ein Beispiel zu nennen. Wie sinnvoll das ist, weiß ich nicht.

Freier Radikaler schrieb:

[...] da immer ein fester rssi wert einer bestimmten entfernung so zugeordnet wird. Das ist aber hochgradig Fehleranfällig. Deswegen der Umweg über die Verhältnisse!

Das erschließt sich mir jetzt allerdings auch nicht auf Anhieb, warum der Umweg besser sein sollte.

Freier Radikaler schrieb:

Der Ansatz über die Verhältnisse ist nicht quark, da er in unterschiedlichsten Papers so vorgeschlagen und umgesetzt wurde.

Ok.

So, und wo ist nun das Problem? Du hast einen Haufen von m Gleichungen, die jeweils die gesuchten x/y/z-Koordinaten mit etwas verbinden, was sich aus Messergebnissen ergibt. Eine solche Gleichung könnte diese sein:

$\frac{\left|| P\_x - P\_1 \right||}{\left|| P\_x - P\_1 \right|| + \left|| P\_x - P\_2 \right||} \approx \frac{4}{3+4}$

Nimmst du alle Gleichungen zusammen, kannst du das so schreiben:

f : \mathbb{R}^3 \mapsto \mathbb{R}^m \; \mbox{ mit } \; f(p) = b

mit Vektor p als gesuchten Punkt und Vektor b, der die Messergebnisse irgendwie enthält. Der Levenberg-Marquardt-Algorithmus kann dir jetzt eine Lösung für p liefern, die

$\left|| f(p) - b \right||_2$

minimiert. Das einzig komplizierte dabei ist die Berechnung der Jacobi-Matrix J von f für ein bestimmtes p. Implementieren tut man Levenberg-Marquardt auch nicht unbedingt selbst, und wenn doch, dann bitteschön auf Basis schon vorhandener Linear-Algebra-Pakete (wie z.B. Eigen). In diesem Fall (mit nur drei Unbekannten) bekommt man das noch recht einfach hin mit den Gleichungssystemen. Irgendwo muss die suche mit einem p₀ beginnen. Und dann passiert bei den Iterationen über k "nur" folgendes:

$\left(J^T\_k J\_k + \lambda\_k I \right) s\_k = J^T\_k \left( b - f(p\_k) \right)$

nach s_k lösen (Korrekturschritt), wobei $\lambda_k I$ der Levenberg-Marqaudt-spezifische Regularisierungsparameter ist. Die Matrix ist positiv definit. Du kannst hier also mit einer Cholesky-Zerlegung rangehen. Und dann das update:

$p_{k+1} = p\_k + s\_k$

Ach, und

$J\_k = \left( \frac{ \partial f\_i(p\_k) }{ \partial (p\_k)\_j } \right)\_{i,j}$

Das schöne ist allerdings, dass Du nicht J und b als Matrix und Vektor komplett aufbauen musst, sondern direkt J^T J und J^T b berechnen kannst. Damit sparst du Platz.

Mit den Ableitungen kann dir "automatic differentiation" helfen. Aber das ist nochmal eine ganz andere Baustelle.

Freier Radikaler

Stimmt, die dBm-Zahlen nehmen nicht quadratisch ab sondern die zu grunde liegenden Sendeleistungspegel. Zu berechnen mit 10^(dBm-Zahl / 10).

Konkret verhalten sich die verhältnisse von P1Geschätzt / P2Geschätzt = Wurzel ( (10^(dBm-Zahl2 / 10)) / (10^(dBm-Zahl / 10)) ).
Zu grunde liegt da das Abstandsgesetz: http://de.wikipedia.org/wiki/Abstandsgesetz

Das was du da als Verfahren ansprichst nennt sich soweit ich das weis aber "Gradientenverfahren". Nach dem habe ich das jetzt implementiert. Die Jakobimatrix war mittels Wolfram alpha schnell aufgestellt. Man hätte das auch allgemeiner wie du es gesagt hast machen können, doch in meinem Fall wird es immer |R³ sein. Es sei denn uns schenkt jemand eine weitere Dimension. Das wird jedoch in absehbarer Zeit wohl nicht der Fall sein :D. Bei der Ableitung tauchen auch viele Konstanten wieder auf wie etwa Entferungen und deren Wurzeln. Das lässt sich sehr gut in jedem Iterationsschritt optimiert vorberechnen, sodass die Berechnung des Gradienten aus der Jakobimatrix sehr performant machbar ist.
Mit Hilfe einer optimalen Schrittweitenbestimmung ist man auch schnell mit der Iteration fertig.
Da man so allerdings nur lokale Minima findet, und das konkrete lokale Minima von der Wahl des Startpunktes abhängt, wäre es interessant zu wissen, wie viele Minima es global betrachtet gibt.
Den Fall zu erkennen in dem es nur ein lokales Minima gibt welches dann dadurch auch das globale Minima wäre, könnte man durch Bildung der Hesse Matrix und bestimmen aller Teildeterminanten > 0 von dieser in Erfahrung bringen. Das ist jedoch sehr umständlich und bringt im Programm selbst nur wenig Mehrwert.

Gibts da vllt. was besseres wie man die Anzahl lokaler Minimas berechnen kann?
Wenn man die Zahl wüsste, könnte man so oft von einem anderen Punkt aus anfangen, bis man alle lokalen Minimas gefunden hat, um dann das tiefste Minima herauszunehmen...

knivil

Gibts da vllt. was besseres wie man die Anzahl lokaler Minimas berechnen kann?

Nein. Die Anzahl der Minima wird durch die Modellfunktion und dessen zu schaetzenden Parametern festgelegt. Im Idealfall sollte es nur ein Minima geben.

Und wenn deine Groesse mit steigender Entfernung streng monoton faellt und du nur eine Quelle hast, dann gibt es auch nur ein Minima/Maxima bzw. einen Punkt wo sie sein kann.