Quais métodos simples existem para amostrar adaptativamente uma função 2D?

Eu tenho uma função bidimensional cujos valores eu gostaria de provar. A função é muito cara de calcular e tem uma forma complexa, portanto, preciso encontrar uma maneira de obter o máximo de informações sobre sua forma usando o menor número de pontos de amostra.$f(x,y)$

Que bons métodos existem para fazer isso?

O que eu tenho até agora

• Começo com um conjunto de pontos existente em que já calculei o valor da função (isso poderia ser uma estrutura quadrada de pontos ou outra coisa).

• Então eu calculo uma triangulação de Delaunay desses pontos.

• Se dois pontos vizinhos na triangulação de Delaunay são distantes o suficiente ( ) e o valor da função difere suficientemente neles ( > Δ f ), insiro um novo ponto no meio deles. Eu faço isso para cada par de pontos vizinho.$> \Delta X$$> \Delta f$

O que há de errado com esse método?

Bem, funciona relativamente bem, mas em funções semelhantes a essa não é ideal porque os pontos de amostra tendem a "pular" a cordilheira e nem percebem que está lá.

Produz resultados como este (se a resolução da grade de pontos inicial for suficientemente aproximada):

Este gráfico acima mostra os pontos onde o valor da função é calculado (na verdade, as células Voronoi ao seu redor).

Este gráfico acima mostra a interpolação linear gerada a partir dos mesmos pontos e a compara ao método de amostragem embutido do Mathematica (para aproximadamente a mesma resolução inicial).

Como melhorá-lo?

Penso que a questão principal aqui é que meu método decide se deve adicionar um ponto de refinamento ou não com base no gradiente.

Seria melhor levar em consideração a curvatura ou pelo menos a segunda derivada ao adicionar pontos de refinamento.

Questão

O que é uma maneira muito simples de implementar para levar em consideração a segunda derivada ou curvatura quando os locais dos meus pontos não estão restritos? (Eu não tenho necessariamente uma estrutura quadrada de pontos de partida, isso idealmente deve ser geral.)

Ou que outras maneiras simples existem para calcular a posição dos pontos de refinamento da maneira ideal?

Vou implementar isso no Mathematica, mas essa pergunta é principalmente sobre o método. Para o bit "fácil de implementar", é importante contar que eu estou usando o Mathematica (isto é, isso foi fácil de fazer até agora porque tem um pacote para fazer a triangulação de Delaunay)

Que problema prático eu estou aplicando isso a

Estou calculando um diagrama de fases. Tem uma forma complexa. Em uma região, seu valor é 0, em outra região, está entre 0 e 1. Há um salto acentuado entre as duas regiões (é descontínuo). Na região onde a função é maior que zero, existem variações suaves e algumas descontinuidades.

O valor da função é calculado com base em uma simulação de Monte Carlo; portanto, ocasionalmente, é esperado um valor ou ruído incorreto da função (isso é muito raro, mas ocorre um grande número de pontos, por exemplo, quando o estado estacionário não é atingido devido algum fator aleatório)

_{Já perguntei isso no Mathematica.SE, mas não consigo vinculá-lo porque ele ainda está na versão beta privada. Esta pergunta aqui é sobre o método, não a implementação.}

Responder para @suki

Esse é o tipo de divisão que você sugere, ou seja, colocando um novo ponto no meio dos triângulos?

Minha preocupação aqui é que parece exigir um tratamento especial nas bordas da região, caso contrário, resultará em triângulos muito longos e muito finos, como mostrado acima. Você corrigiu isso?

ATUALIZAR

Um problema que aparece tanto no método que descrevo quanto na sugestão de @ suki de subdividir com base em triângulos e colocar os pontos de subdivisão dentro do triângulo é que, quando há descontinuidades (como no meu problema), recalcular a triangulação de Delaunay após uma etapa pode faça com que os triângulos mudem e talvez alguns grandes triângulos apareçam com diferentes valores de função nos três vértices.

Aqui estão dois exemplos:

O primeiro mostra o resultado final ao amostrar em torno de uma descontinuidade direta. O segundo mostra a distribuição do ponto de amostragem para um caso semelhante.

Que maneiras simples existem para evitar isso? Atualmente, estou simplesmente subdividindo os egdes que desaparecem após uma retriangulação, mas isso parece um hack e precisa ser feito com cuidado, pois no caso de malhas simétricas (como uma grade quadrada) existem várias triangulações válidas de Delaunay, portanto as bordas podem mudar aleatoriamente após retriangulação.

Eu trabalhei em um problema semelhante a isso há algum tempo.

Eu acho que a principal diferença entre as nossas implementações é que eu estava escolhendo onde adicionar pontos com base nos triângulos, não nas arestas. Também escolho novos pontos dentro dos triângulos, em vez de nas bordas.

Tenho a sensação de que adicionar pontos dentro dos triângulos tornaria mais eficiente, dando um pequeno aumento na distância média dos pontos antigos ao novo.

De qualquer forma, outra coisa interessante sobre o uso de triângulos em vez de arestas é que ela fornece uma estimativa do vetor gradiente, em vez da inclinação ao longo dessa aresta em particular.

No meu código do matlab, usei uma classe base para cuidar da maioria das máquinas, com alguns métodos abstratos:

• weight(self) para decidir a prioridade para a qual os triângulos serão subdivididos a seguir.
• choosePoints(self,npoints = "auto") decidir novos pontos a serem avaliados com base no peso de cada triângulo.

Achei essa configuração muito flexível:

• definir a função de uma subclasse weight()na área do triângulo produz uma densidade de malha constante.
• A configuração weight()para calcular o valor médio da função vezes a área do triângulo fornece uma espécie de amostragem de probabilidade quase aleatória.
• using var(triangle.zs)poderia fazer, para funções com saída binária, o que eu sinto é uma generalização de uma pesquisa de bissecção para mais de uma dimensão.
• O uso area + var(triangle.zs)foi bastante eficaz para colocar uma densidade constante em todos os lugares e um aumento da densidade ao longo de qualquer declive (quase o que você tem agora).

Eu usei a variação dos valores z, para aproximar a importância dos efeitos de primeira ordem (inclinação), porque a variação nunca chegará ao infinito como a inclinação pode.

Para o último exemplo, a densidade do fundo foi boa porque eu estava procurando por blobs descontínuos de alto valor em um espaço de baixo valor. Então, ele preenchia lentamente toda a malha e, quando encontrava um blob, concentrava-se em seguir a borda do blob por toda a volta por causa do alto peso que colocava no gradiente (e que preenchia apenas os ntriângulos superiores em cada iteração). No final, eu sabia que não havia blobs (de forma razoável) (ou orifícios nos blobs) com um tamanho maior que a densidade de malha de fundo resultante.

Como você, obtive alguns pontos negativos nos meus resultados, eles não foram um problema para mim, porque o erro era tal que, se você repetisse os pontos próximos, provavelmente dariam a resposta correta. Acabaria com pontinhos de maior densidade de malha em torno dos meus pontos negativos.

O que quer que você faça, eu sempre recomendo fazer os pesos relacionados ao tamanho do triângulo, para que, sendo o resto igual, os triângulos grandes sejam divididos primeiro.

Talvez uma solução para você seja levar minha abordagem um passo adiante e, em vez de avaliar triângulos com base no conteúdo dessa célula triangular, avalie com base nesse e nos três triângulos adjacentes.

Isso conterá informações suficientes para obter uma estimativa da matriz hessiana completa. você pode obtê-lo fazendo um ajuste mínimo de quadrados z = c1*x + C2*y c11*x^2+c12*x*y+c22*y^2sobre todos os vértices nos triângulos de interesse (centralize o sistema de coordenadas no triângulo primeiro).

Eu não usaria o gradiente ou o Hessiano (aquelas constantes) diretamente, porque eles irão para o infinito com uma descontinuidade.

Talvez o erro da soma quadrada dos valores z em relação a uma aproximação planar desses pontos seja uma medida útil de quão interessantes são os efeitos de segunda ordem.

Atualizada:

Isso me parece razoável.

Eu nunca cheguei a um revestimento especial nas bordas. Isso me incomodou um pouco, mas pelo que eu estava fazendo, bastava começar com muitos pontos nas bordas.

mais elegante seria combinar nossas duas abordagens, ponderar arestas e triângulos. Então, se uma aresta for muito longa, corte-a ao meio ... Gosto da maneira como esse conceito se generaliza para dimensões mais altas (mas os números ficam grandes rapidamente) ...

Mas como você não espera que o corpo principal da malha tenha triângulos de alta proporção, então você pode usar uma função como a função de fronteira livre do Matlab para encontrar o limite e executar o mesmo algoritmo em uma dimensão a menos no limite. Se bem feito, em um cubo, por exemplo, você pode obter a mesma densidade de malha nas bordas, nas faces e dentro do cubo. Interessante.

Uma coisa para a qual nunca encontrei uma boa solução foi o fato de que minha versão nunca exploraria fora do casco convexo do conjunto de pontos inicial.

Acho que o principal problema em sua heurística é que você está considerando o gradiente em apenas uma dimensão e, portanto, em regiões onde o dfdx é pequeno, mas o dfdy é grande (como acontece no meio do seu exemplo), você perderá pontos ao olhar na dimensão "errada".

Uma solução rápida seria considerar conjuntos de quatro pontos, tomando seu centro de gravidade e aproximando | dfdx | + | dfdy | usando esses quatro pontos. Outra alternativa é pegar três pontos (ou seja, um triângulo) e obter o gradiente máximo da superfície sobre eles.