Welche einfachen Methoden gibt es, um eine 2D-Funktion adaptiv abzutasten?

Ich habe eine zweidimensionale Funktion deren Werte ich abtasten möchte. Die Berechnung der Funktion ist sehr aufwendig und weist eine komplexe Form auf. Ich muss daher einen Weg finden, um mit der geringsten Anzahl von Abtastpunkten die meisten Informationen über ihre Form zu erhalten.$f(x,y)$

Welche guten Methoden gibt es dafür?

Was ich bisher habe

• Ich gehe von einer vorhandenen Menge von Punkten aus, bei denen ich den Funktionswert bereits berechnet habe (dies kann ein quadratisches Punktegitter oder etwas anderes sein).

• Dann berechne ich eine Delaunay-Triangulation dieser Punkte.

• Wenn zwei benachbarte Punkte in der Delaunay-Triangulation weit genug voneinander entfernt sind ( ) und sich der Funktionswert darin ausreichend unterscheidet ( ), füge ich einen neuen Punkt in der Mitte dazwischen ein. Ich mache das für jedes benachbarte Punktpaar.$> \Delta X$$> \Delta f$

Was ist los mit dieser Methode?

Nun, es funktioniert relativ gut, aber bei ähnlichen Funktionen ist es nicht ideal, da die Abtastpunkte dazu neigen, über den Grat zu springen und nicht einmal zu bemerken, dass sie dort sind.

Das Ergebnis ist wie folgt (wenn die Auflösung des Anfangspunktrasters ausreichend grob ist):

Dieses Diagramm oben zeigt die Punkte, an denen der Funktionswert berechnet wird (tatsächlich Voronoi-Zellen um sie herum).

Dieses Diagramm oben zeigt die lineare Interpolation, die aus denselben Punkten generiert wurde, und vergleicht sie mit der in Mathematica integrierten Abtastmethode (für ungefähr dieselbe Startauflösung).

Wie kann man es verbessern?

Ich denke, das Hauptproblem hierbei ist, dass meine Methode basierend auf dem Farbverlauf entscheidet, ob ein Verfeinerungspunkt hinzugefügt wird oder nicht.

Es ist besser, die Krümmung oder zumindest die zweite Ableitung beim Hinzufügen von Verfeinerungspunkten zu berücksichtigen.

Frage

Was ist eine sehr einfach zu implementierende Methode, um die zweite Ableitung oder Krümmung zu berücksichtigen, wenn die Positionen meiner Punkte überhaupt nicht eingeschränkt sind? (Ich habe nicht unbedingt ein quadratisches Gitter mit Startpunkten, dies sollte idealerweise allgemein sein.)

Oder welche anderen einfachen Möglichkeiten gibt es, um die Position von Verfeinerungspunkten optimal zu berechnen?

Ich werde dies in Mathematica implementieren, aber bei dieser Frage geht es hauptsächlich um die Methode. Für das "einfach zu implementierende" Bit zählt allerdings, dass ich Mathematica verwende (dh dies war bisher einfach, da es ein Paket für die Delaunay-Triangulation enthält).

Auf welches praktische Problem wende ich das an?

Ich berechne ein Phasendiagramm. Es hat eine komplexe Form. In einer Region ist der Wert 0, in einer anderen zwischen 0 und 1. Zwischen den beiden Regionen besteht ein starker Sprung (diskontinuierlich). In dem Bereich, in dem die Funktion größer als Null ist, gibt es sowohl eine sanfte Variation als auch ein paar Diskontinuitäten.

Der Funktionswert wird auf Basis einer Monte-Carlo-Simulation berechnet, so dass gelegentlich ein falscher Funktionswert oder Rauschen zu erwarten ist (dies ist sehr selten, aber bei einer großen Anzahl von Punkten tritt dies auf, z. B. wenn der eingeschwungene Zustand nicht erreicht wird) ein zufälliger Faktor)

_{Ich habe dies bereits auf Mathematica.SE gefragt , kann aber keinen Link dazu erstellen, da es sich noch in der privaten Beta befindet. Bei dieser Frage geht es um die Methode, nicht um die Implementierung.}

Antworte auf @suki

Ist dies die Art der Unterteilung, die Sie vorschlagen, dh einen neuen Punkt in die Mitte der Dreiecke zu setzen?

Meine Sorge ist hier, dass es an den Rändern der Region anscheinend eine besondere Behandlung erfordert, da es sonst sehr lange und sehr dünne Dreiecke gibt, wie oben gezeigt. Hast du das korrigiert?

AKTUALISIEREN

Ein Problem, das sowohl bei der von mir beschriebenen Methode als auch bei dem Vorschlag von @ suki auftritt, eine Unterteilung basierend auf Dreiecken vorzunehmen und die Unterteilungspunkte innerhalb des Dreiecks zu setzen, besteht darin, dass bei Diskontinuitäten (wie in meinem Problem) die Delaunay-Triangulation nach einem Schritt möglicherweise neu berechnet wird bewirken, dass sich Dreiecke ändern und möglicherweise große Dreiecke mit unterschiedlichen Funktionswerten in den drei Scheitelpunkten angezeigt werden.

Hier sind zwei Beispiele:

Die erste zeigt das Endergebnis beim Abtasten um eine gerade Diskontinuität. Die zweite zeigt die Verteilung der Abtastpunkte für einen ähnlichen Fall.

Welche einfachen Möglichkeiten gibt es, dies zu vermeiden? Gegenwärtig unterteile ich einfach jene Egdes, die nach einer Retriangulation verschwinden, aber dies fühlt sich wie ein Hack an und muss mit Sorgfalt durchgeführt werden, da es bei symmetrischen Netzen (wie einem quadratischen Gitter) mehrere gültige Delaunay-Triangulationen gibt, daher können sich Kanten ändern zufällig nach retriangulation.

Ich habe vor einiger Zeit an einem ähnlichen Problem gearbeitet.

Ich denke, der Hauptunterschied zwischen unseren Implementierungen besteht darin, dass ich anhand der Dreiecke und nicht anhand der Kanten ausgewählt habe, wo Punkte hinzugefügt werden sollen. Ich wähle auch neue Punkte innerhalb der Dreiecke anstatt an den Kanten.

Ich habe das Gefühl, dass das Hinzufügen von Punkten innerhalb der Dreiecke die Effizienz erhöhen würde, indem der durchschnittliche Abstand zwischen alten und neuen Punkten geringfügig erhöht wird.

Eine weitere nette Sache bei der Verwendung von Dreiecken anstelle von Kanten ist, dass sie eine Schätzung des Gradientenvektors anstelle der Steigung entlang dieser bestimmten Kante liefert.

In meinem Matlab-Code habe ich eine Basisklasse verwendet, um den größten Teil der Maschinerie mit ein paar abstrakten Methoden zu verwalten:

• weight(self) um die Priorität zu bestimmen, für die Dreiecke als nächstes unterteilt werden sollen.
• choosePoints(self,npoints = "auto") neue Punkte zu bestimmen, die basierend auf dem Gewicht jedes Dreiecks bewertet werden sollen.

Ich fand dieses Setup sehr flexibel:

• Wenn Sie die weight()Funktion einer Unterklasse auf die Fläche des Dreiecks einstellen, wird eine konstante Maschendichte erzeugt.
• Die Einstellung weight()zur Berechnung des durchschnittlichen Funktionswerts multipliziert mit der Fläche des Dreiecks ergibt eine Art quasi zufällige Wahrscheinlichkeitsabtastung.
• using var(triangle.zs)könnte für Funktionen mit Binärausgabe eine Verallgemeinerung einer Halbierungssuche auf mehr als eine Dimension bedeuten.
• Das Verwenden von area + var(triangle.zs)war ziemlich effektiv, um überall eine konstante Dichte und eine erhöhte Dichte entlang einer beliebigen Steigung zu erzielen (fast so, wie Sie es jetzt haben).

Ich habe die Varianz der z-Werte verwendet, um die Wichtigkeit von Effekten erster Ordnung (Steigung) zu approximieren, da die Varianz niemals unendlich wird, wie es die Steigung kann.

Für das letzte Beispiel war die Hintergrunddichte gut, weil ich auf der Suche nach diskontinuierlichen Blobs mit hohem Wert in einem Raum mit geringem Wert war. Also füllte es langsam das gesamte Netz aus und wenn es einen Blob fand, konzentrierte es sich darauf, dem Rand des Blobs rundherum zu folgen, da ich das Gefälle mit hohem Gewicht belastete (und nur die oberen nDreiecke ausfüllte) bei jeder Iteration). Am Ende konnte ich feststellen, dass es keine (recht geformten) Blobs (oder Löcher in meinen Blobs) gab, die größer waren als die resultierende Hintergrund-Maschendichte.

Wie Sie, ich habe einige schlechte Punkte in meinen Ergebnissen erhalten, sie waren kein Problem für mich, da der Fehler so war, dass sie wahrscheinlich die richtige Antwort gaben, wenn Sie nahegelegene Punkte erneut ausführen. Ich würde nur mit Streifen von erhöhter Maschendichte um meine schlechten Punkte enden.

Was auch immer Sie tun, ich empfehle immer, die Gewichte in Bezug auf die Dreiecksgröße so zu gestalten, dass bei allen anderen gleichen Werten zuerst große Dreiecke aufgebrochen werden.

Vielleicht besteht eine Lösung für Sie darin, meinen Ansatz einen Schritt weiter zu gehen und statt Dreiecke auf der Grundlage des Inhalts dieser dreieckigen Zelle zu bewerten, auf der Grundlage dieses einen und aller drei benachbarten Dreiecke zu bewerten.

Das wird genügend Informationen enthalten, um eine Schätzung der gesamten hessischen Matrix zu erhalten. Sie können dies erreichen, indem Sie eine Anpassung der kleinsten Quadrate z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2über alle Scheitelpunkte in den Dreiecken von Interesse durchführen (zentrieren Sie das Koordinatensystem zuerst auf dem Dreieck).

Ich würde den Gradienten oder das Hessische (diese Konstanten) nicht direkt verwenden, weil sie bei einer Diskontinuität ins Unendliche gehen.

Vielleicht wäre der Summenquadratfehler der z-Werte in Bezug auf eine planare Approximation dieser Punkte ein nützliches Maß dafür, wie interessant Effekte zweiter Ordnung sind.

Aktualisiert:

Das erscheint mir vernünftig.

Ich bin eigentlich nie dazu gekommen, die Kanten speziell zu umhüllen. Es hat mich ein wenig gestört, aber für das, was ich tat, war es genug, nur mit vielen Punkten an den Rändern zu beginnen.

eleganter wäre es, unsere beiden Ansätze zu kombinieren, Kanten und Dreiecke zu gewichten. Wenn eine Kante dann zu lang ist, halbieren Sie sie ... Ich mag die Art und Weise, wie das Konzept auf höhere Dimensionen verallgemeinert wird (aber die Zahlen werden schnell groß) ...

Da Sie jedoch nicht erwarten, dass der Hauptteil des Netzes Dreiecke mit hohem Seitenverhältnis aufweist, können Sie eine Funktion wie die Matlab- Funktion für die freie Grenze verwenden, um die Grenze zu finden, und dann denselben Algorithmus in einer Dimension weniger auf der Grenze ausführen. Wenn Sie dies richtig machen, z. B. bei einem Würfel, erhalten Sie an den Kanten, an den Flächen und im Inneren des Würfels die gleiche Maschendichte. Interessant.

Eine Sache, für die ich nie eine gute Lösung gefunden habe, war die Tatsache, dass meine Version niemals außerhalb der konvexen Hülle des Anfangspunktsatzes erkunden würde.

Ich denke, das Hauptproblem in Ihrer Heuristik besteht darin, dass Sie den Gradienten nur in einer Dimension betrachten und daher in Regionen, in denen dfdx klein, aber dfdy groß ist (wie in der Mitte Ihres Beispiels), Punkte beim Betrachten verpassen in der "falschen" Dimension.

Eine schnelle Lösung wäre, Mengen von vier Punkten zu betrachten, deren Schwerpunkt zu bestimmen und | dfdx | + | dfdy | anzunähern mit diesen vier Punkten. Eine andere Alternative ist, drei Punkte (dh ein Dreieck) zu nehmen und den maximalen Gradienten der Oberfläche über ihnen zu nehmen.