¿Qué métodos simples existen para muestrear adaptativamente una función 2D?

Tengo una función bidimensional cuyos valores me gustaría muestrear. La función es muy costosa de calcular y tiene una forma compleja, por lo que necesito encontrar una manera de obtener la mayor información sobre su forma utilizando el menor número de puntos de muestra.$f(x,y)$

¿Qué buenos métodos hay para hacer esto?

Lo que tengo hasta ahora

• Comienzo desde un conjunto de puntos existente donde ya he calculado el valor de la función (esto podría ser una red cuadrada de puntos u otra cosa).

• Luego calculo una triangulación de Delaunay de estos puntos.

• Si dos puntos vecinos en la triangulación de Delaunay están lo suficientemente lejos ( ) y el valor de la función difiere lo suficiente en ellos ( > Δ f ), entonces inserto un nuevo punto a medio camino entre ellos. Hago esto para cada par de puntos vecinos.$> \Delta X$$> \Delta f$

¿Qué tiene de malo este método?

Bueno, funciona relativamente bien, pero en funciones similares a esta no es ideal porque los puntos de muestra tienden a "saltar" sobre la cresta y ni siquiera se dan cuenta de que está allí.

Produce resultados como este (si la resolución de la cuadrícula de puntos inicial es suficientemente aproximada):

Este gráfico anterior muestra los puntos donde se calcula el valor de la función (en realidad, las celdas de Voronoi a su alrededor).

Este gráfico anterior muestra la interpolación lineal generada a partir de los mismos puntos y la compara con el método de muestreo incorporado de Mathematica (para aproximadamente la misma resolución inicial).

¿Cómo mejorarlo?

Creo que el problema principal aquí es que mi método decide si agregar un punto de refinamiento o no en función del gradiente.

Sería mejor tener en cuenta la curvatura o al menos la segunda derivada al agregar puntos de refinamiento.

Pregunta

¿Cuál es una forma muy sencilla de implementar para tener en cuenta la segunda derivada o curvatura cuando las ubicaciones de mis puntos no están limitadas en absoluto? (No necesariamente tengo una red cuadrada de puntos de partida, esto idealmente debería ser general).

¿O qué otras formas simples existen para calcular la posición de los puntos de refinamiento de manera óptima?

Voy a implementar esto en Mathematica, pero esta pregunta es principalmente sobre el método. Para el bit "fácil de implementar", sí cuenta que estoy usando Mathematica (es decir, esto fue fácil de hacer hasta ahora porque tiene un paquete para hacer la triangulación de Delaunay)

¿A qué problema práctico estoy aplicando esto?

Estoy calculando un diagrama de fase. Tiene una forma compleja. En una región su valor es 0, en otra región está entre 0 y 1. Hay un salto brusco entre las dos regiones (es discontinuo). En la región donde la función es mayor que cero, existe una variación suave y un par de discontinuidades.

El valor de la función se calcula en base a una simulación de Monte Carlo, por lo que ocasionalmente se espera un valor de función incorrecto o ruido (esto es muy raro, pero ocurre en un gran número de puntos, por ejemplo, cuando no se alcanza el estado estable debido algún factor aleatorio)

_{Ya he preguntado esto en Mathematica.SE, pero no puedo vincularlo porque todavía está en beta privada. Esta pregunta aquí es sobre el método, no la implementación.}

Responder a @suki

¿Es este el tipo de división que sugiere, es decir, poner un nuevo punto en el medio de los triángulos?

Mi preocupación aquí es que parece requerir un manejo especial en los bordes de la región, de lo contrario dará triángulos muy largos y muy delgados, como se muestra arriba. ¿Corrigiste esto?

ACTUALIZAR

Un problema que aparece tanto con el método que describo como con la sugerencia de @ suki de subdividir en base a triángulos y colocar los puntos de subdivisión dentro del triángulo es que cuando hay discontinuidades (como en mi problema), volver a calcular la triangulación de Delaunay después de un paso puede hacen que los triángulos cambien y tal vez aparezcan algunos triángulos grandes que tienen diferentes valores de función en los tres vértices.

Aquí hay dos ejemplos:

El primero muestra el resultado final al muestrear alrededor de una discontinuidad directa. El segundo muestra la distribución del punto de muestreo para un caso similar.

¿Qué formas simples hay para evitar esto? Actualmente, simplemente estoy subdividiendo esos egdes que desaparecen después de una retriangulación, pero esto se siente como un truco y debe hacerse con cuidado ya que en el caso de mallas simétricas (como una cuadrícula cuadrada) hay varias triangulaciones válidas de Delaunay, por lo tanto, los bordes pueden cambiar aleatoriamente después de la retriangulación.

Trabajé en un problema similar a esto hace un tiempo.

Creo que la principal diferencia entre nuestras implementaciones es que estaba eligiendo dónde agregar puntos basados ​​en los triángulos, no en los bordes. También elijo nuevos puntos dentro de los triángulos en lugar de en los bordes.

Tengo la sensación de que agregar puntos dentro de los triángulos lo haría más eficiente al dar un pequeño aumento en la distancia promedio de los puntos antiguos al nuevo.

De todos modos, otra cosa interesante sobre el uso de triángulos en lugar de bordes es que proporciona una estimación del vector de gradiente, en lugar de la pendiente a lo largo de este borde en particular.

En mi código matlab utilicé una clase base para cuidar la mayoría de la maquinaria, con algunos métodos abstractos:

• weight(self) para decidir la prioridad para qué triángulos subdividir a continuación.
• choosePoints(self,npoints = "auto") para decidir nuevos puntos para evaluar en función del peso de cada triángulo.

Encontré esta configuración muy flexible:

• establecer una función de subclase weight()en el área del triángulo produce una densidad de malla constante.
• la configuración weight()para calcular el valor de la función promedio por el área del triángulo da una especie de muestreo de probabilidad cuasialeatorio.
• usando var(triangle.zs)podría hacer, para funciones que tienen salida binaria, lo que siento es una generalización de una búsqueda de bisección a más de 1 dimensión.
• el uso area + var(triangle.zs)fue bastante efectivo para colocar una densidad constante en todas partes, y una mayor densidad a lo largo de cualquier pendiente (casi lo que tienes ahora).

Utilicé la varianza de los valores z para aproximar la importancia de los efectos de primer orden (pendiente) porque la varianza nunca llegará al infinito como puede hacerlo la pendiente.

Para el último ejemplo, la densidad de fondo era buena porque estaba buscando manchas discontinuas de alto valor en un espacio de bajo valor. Por lo tanto, llenaría lentamente toda la malla y, cuando encontrara una gota, se concentraría en seguir el borde de la gota todo el tiempo debido al alto peso que puse en el gradiente (y que solo llenaba los ntriángulos superiores en cada iteración). Al final pude saber que no había manchas (o agujeros de forma razonable) (o agujeros en mis manchas) de un tamaño mayor que la densidad de malla de fondo resultante.

Al igual que usted, obtuve algunos puntos negativos en mis resultados, no fueron un problema para mí porque el error fue tal que si volviera a ejecutar los puntos cercanos, probablemente darían la respuesta correcta. Simplemente terminaría con manchas de densidad de malla aumentada alrededor de mis puntos negativos.

Hagas lo que hagas, siempre recomiendo hacer los pesos relacionados con el tamaño del triángulo para que, al igual que todos los demás, los triángulos grandes se dividan primero.

Quizás una solución para usted es llevar mi enfoque un paso más allá y, en lugar de evaluar triángulos basados ​​en el contenido de esa celda triangular, evaluar en base a ese y los tres triángulos adyacentes.

Eso contendrá suficiente información para obtener una estimación de la matriz completa de Hesse. puede obtenerlo haciendo un ajuste de mínimos cuadrados z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2sobre todos los vértices en los triángulos de interés (centre primero el sistema de coordenadas en el triángulo).

No usaría el gradiente o Hessian (esas constantes) directamente porque irán al infinito en una discontinuidad.

Tal vez el error de suma cuadrada de los valores z en relación con una aproximación plana de esos puntos sería una medida útil de lo interesantes que son los efectos de segundo orden.

Actualizado:

Eso me parece razonable.

En realidad nunca llegué a revestir los bordes de manera especial. Me molestó un poco, pero por lo que estaba haciendo, fue suficiente para comenzar con muchos puntos alrededor de los bordes.

más elegante sería combinar nuestros dos enfoques, ponderar bordes y triángulos. Luego, si un borde es demasiado largo, córtelo por la mitad ... Me gusta la forma en que ese concepto se generaliza a dimensiones más altas (pero los números se hacen grandes rápidamente) ...

Pero dado que no espera que el cuerpo principal de la malla tenga triángulos de alta relación de aspecto, puede usar una función como la función de límite libre de Matlab para encontrar el límite, luego ejecutar el mismo algoritmo en una dimensión menos en el límite. Si se hace correctamente, en un cubo, por ejemplo, podría obtener la misma densidad de malla en los bordes, en las caras y dentro del cubo. Interesante.

Una cosa para la que nunca encontré una buena solución fue el hecho de que mi versión nunca exploraría fuera del casco convexo del conjunto de puntos inicial.

Creo que el principal problema en su heurística es que está considerando el gradiente solo en una dimensión y, por lo tanto, en regiones donde dfdx es pequeño pero dfdy es grande (como sucede en el medio de su ejemplo), perderá puntos al mirar en la dimensión "equivocada".

Una solución rápida sería considerar conjuntos de cuatro puntos, tomar su centro de gravedad y aproximar | dfdx | + | dfdy | usando esos cuatro puntos. Otra alternativa es tomar tres puntos (es decir, un triángulo) y tomar el gradiente máximo de la superficie sobre ellos.