Quelles méthodes simples existe-t-il pour échantillonner de manière adaptative une fonction 2D?

J'ai une fonction bidimensionnelle dont je voudrais échantillonner les valeurs. La fonction est très coûteuse à calculer et elle a une forme complexe, j'ai donc besoin de trouver un moyen d'obtenir le plus d'informations sur sa forme en utilisant le moins de points d'échantillonnage.$f(x,y)$

Quelles sont les bonnes méthodes pour y parvenir?

Ce que j'ai jusqu'ici

• Je pars d'un ensemble de points existant où j'ai déjà calculé la valeur de la fonction (cela pourrait être un réseau carré de points ou autre chose).

• Ensuite, je calcule une triangulation de Delaunay de ces points.

• Si deux points voisins dans la triangulation de Delaunay sont suffisamment éloignés ( ) et que la valeur de la fonction diffère suffisamment en eux ( > Δ f ), alors j'insère un nouveau point à mi-chemin entre eux. Je fais cela pour chaque paire de points voisins.$> \Delta X$$> \Delta f$

Quel est le problème avec cette méthode?

Eh bien, cela fonctionne relativement bien, mais sur des fonctions similaires à celle-ci, ce n'est pas idéal car les points d'échantillonnage ont tendance à "sauter" la crête et ne remarquent même pas qu'elle est là.

Il produit des résultats comme celui-ci (si la résolution de la grille de points initiale est suffisamment approximative):

Ce graphique ci-dessus montre les points où la valeur de la fonction est calculée (en fait des cellules de Voronoi autour d'eux).

Ce graphique ci-dessus montre l'interpolation linéaire générée à partir des mêmes points et la compare à la méthode d'échantillonnage intégrée de Mathematica (pour environ la même résolution de départ).

Comment l'améliorer?

Je pense que le principal problème ici est que ma méthode décide d'ajouter ou non un point de raffinement en fonction du gradient.

Il serait préférable de prendre en compte la courbure ou au moins la dérivée seconde lors de l'ajout de points de raffinement.

Question

Quelle est une manière très simple à mettre en œuvre pour prendre en compte la dérivée seconde ou la courbure lorsque les emplacements de mes points ne sont pas du tout contraints? (Je n'ai pas nécessairement un réseau carré de points de départ, cela devrait idéalement être général.)

Ou quels autres moyens simples existe-t-il pour calculer la position des points de raffinement de manière optimale?

Je vais implémenter cela dans Mathematica, mais cette question concerne principalement la méthode. Pour le bit "facile à implémenter", cela compte cependant que j'utilise Mathematica (c'est-à-dire que c'était facile à faire jusqu'à présent car il a un paquet pour faire la triangulation de Delaunay)

À quel problème pratique j'applique ceci

Je calcule un diagramme de phases. Il a une forme complexe. Dans une région, sa valeur est 0, dans une autre région, elle est comprise entre 0 et 1. Il y a un net saut entre les deux régions (c'est discontinu). Dans la région où la fonction est supérieure à zéro, il y a à la fois une certaine variation en douceur et quelques discontinuités.

La valeur de la fonction est calculée sur la base d'une simulation de Monte Carlo, donc parfois une valeur de fonction ou un bruit incorrect est à prévoir (cela est très rare, mais pour un grand nombre de points, cela se produit, par exemple lorsque l'état stationnaire n'est pas atteint en raison un facteur aléatoire)

_{J'ai déjà posé cette question sur Mathematica.SE mais je ne peux pas le lier car il est toujours en version bêta privée. Cette question concerne la méthode, pas la mise en œuvre.}

Répondre à @suki

Est-ce le type de division que vous proposez, c'est-à-dire mettre un nouveau point au milieu des triangles?

Ma préoccupation ici est qu'elle semble nécessiter une manipulation spéciale aux bords de la région, sinon elle donnera des triangles très longs et très fins, comme indiqué ci-dessus. Avez-vous corrigé cela?

MISE À JOUR

Un problème qui apparaît à la fois avec la méthode que je décris et avec la suggestion de @ suki de mettre la subdivision basée sur des triangles et de placer les points de subdivision à l'intérieur du triangle est que lorsqu'il y a des discontinuités (comme dans mon problème), recalculer la triangulation de Delaunay après une étape peut faire changer les triangles et peut-être apparaître de gros triangles qui ont des valeurs de fonction différentes dans les trois sommets.

Voici deux exemples:

Le premier montre le résultat final lors de l'échantillonnage autour d'une discontinuité droite. La seconde montre la distribution des points d'échantillonnage pour un cas similaire.

Quels sont les moyens simples d'éviter cela? Actuellement, je subdivise simplement les egdes qui disparaissent après une retriangulation, mais cela ressemble à un hack et doit être fait avec soin car dans le cas des maillages symétriques (comme une grille carrée), il existe plusieurs triangulations Delaunay valides, donc les bords peuvent changer au hasard après retriangulation.

J'ai travaillé sur un problème similaire à celui-ci il y a quelque temps.

Je pense que la principale différence entre nos implémentations est que je choisissais où ajouter des points en fonction des triangles, pas des bords. Je choisis également de nouveaux points à l'intérieur des triangles plutôt que sur les bords.

J'ai le sentiment que l'ajout de points à l'intérieur des triangles le rendrait plus efficace en donnant une petite augmentation de la distance moyenne entre les anciens points et les nouveaux.

Quoi qu'il en soit, une autre bonne chose à propos de l'utilisation de triangles au lieu d'arêtes est qu'elle donne une estimation du vecteur de gradient, au lieu de la pente le long de cette arête particulière.

Dans mon code matlab, j'ai utilisé une classe de base pour prendre en charge la plupart des machines, avec quelques méthodes abstraites:

• weight(self) pour décider de la priorité pour laquelle les triangles à subdiviser ensuite.
• choosePoints(self,npoints = "auto") pour décider de nouveaux points à évaluer en fonction du poids de chaque triangle.

J'ai trouvé cette configuration très flexible:

• définition d'une sous-classe weight() sur la zone du triangle produit une densité de maillage constante.
• réglage weight() pour calculer la valeur moyenne de la fonction multipliée par l'aire du triangle donne une sorte d'échantillonnage de probabilité quasi aléatoire.
• en utilisant var(triangle.zs) pourrait faire, pour les fonctions qui ont une sortie binaire, ce que je ressens est une généralisation d'une recherche de bissection à plus d'une dimension.
• l'utilisation area + var(triangle.zs)était assez efficace pour mettre une densité constante partout et une densité accrue le long de n'importe quelle pente (presque ce que vous avez maintenant).

J'ai utilisé la variance des valeurs z, pour approximer l'importance des effets de premier ordre (pente) car la variance n'atteindra jamais l'infini comme le peut la pente.

Pour le dernier exemple, la densité de fond était bonne parce que je cherchais des taches discontinues de haute valeur dans un espace de faible valeur. Donc, il remplirait lentement tout le maillage et lorsqu'il trouverait un blob, il se concentrerait sur le suivi du bord du blob tout autour en raison du poids élevé que j'ai mis sur le dégradé (et qu'il ne remplissait que le hautn triangles à chaque itération). À la fin, je pouvais savoir qu'il n'y avait pas de taches (de forme raisonnable) (ou trous dans mes gouttes) d'une taille supérieure à la densité de maillage d'arrière-plan résultante.

Comme vous, j'ai obtenu quelques mauvais points dans mes résultats, ils ne m'ont pas posé de problème car l'erreur était telle que si vous réexécutiez des points proches, ils donneraient probablement la bonne réponse. Je finirais juste avec des taches de densité de maillage accrue autour de mes mauvais points.

Quoi que vous fassiez, je recommande toujours de faire les poids liés à la taille du triangle afin que, toutes choses étant égales par ailleurs, les grands triangles soient d'abord brisés.

Peut-être qu'une solution pour vous est d'aller plus loin dans mon approche et au lieu d'évaluer les triangles en fonction du contenu de cette cellule triangulaire, d'évaluer en fonction de celui-ci et des trois triangles adjacents.

Cela contiendra suffisamment d'informations pour obtenir une estimation de la matrice hessoise complète. vous pouvez l'obtenir en faisant un ajustement des moindres carrés z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2sur tous les sommets des triangles d'intérêt (centrez d'abord le système de coordonnées sur le triangle).

Je n'utiliserais pas le gradient ou la Hesse (ces constantes) directement car ils iront à l'infini à une discontinuité.

Peut-être que l'erreur de somme des carrés des valeurs z par rapport à une approximation planaire de ces points serait une mesure utile de la façon dont les effets du second ordre sont intéressants.

Mis à jour:

Cela me semble raisonnable.

Je n'ai jamais réussi à mettre un boîtier spécial sur les bords. Cela me dérangeait un peu mais pour ce que je faisais, c'était suffisant pour commencer avec beaucoup de points sur les bords.

plus élégant serait de combiner nos deux approches, la pondération des bords et des triangles. Ensuite, si un bord est trop long, coupez-le en deux ... J'aime la façon dont ce concept se généralise à des dimensions plus élevées (mais les chiffres grossissent rapidement) ...

Mais puisque vous ne vous attendez pas à ce que le corps principal du maillage ait des triangles à rapport d'aspect élevé, vous pouvez donc utiliser une fonction comme la fonction freeboundary de Matlab pour trouver la frontière, puis exécuter le même algorithme dans une dimension de moins sur la frontière. Si c'est bien fait, sur un cube par exemple, vous pourriez obtenir la même densité de maillage sur les bords, sur les faces et à l'intérieur du cube. Intéressant.

Une chose pour laquelle je n'ai jamais trouvé de bonne solution était le fait que ma version n'explorerait jamais en dehors de la coque convexe de l'ensemble de points initial.

Je pense que le principal problème de votre heuristique est que vous ne considérez le gradient que dans une seule dimension et donc, dans les régions où dfdx est petit mais dfdy est grand (comme cela se produit au milieu de votre exemple), vous manquerez des points en regardant dans la "mauvaise" dimension.

Une solution rapide consisterait à considérer des ensembles de quatre points, en prenant leur centre de gravité et en rapprochant | dfdx | + | dfdy | en utilisant ces quatre points. Une autre alternative consiste à prendre trois points (c'est-à-dire un triangle) et à prendre le gradient maximum de la surface sur eux.