Wie wählt man eine gleichmäßig verteilte Teilmenge eines teilweise dichten Datensatzes?

Question

P ist eine n * d-Matrix, die n d-dimensionale Proben enthält. P in einigen Bereichen ist mehrere Male dichter als andere. Ich möchte eine Teilmenge von P auswählen, in der der Abstand zwischen irgendwelchen Paaren von Proben mehr als d0 ist, und ich brauche es, um über die Fläche verteilt zu werden. Alle Abtastwerte haben dieselbe Priorität, und es ist nicht notwendig, irgendetwas zu optimieren (z. B. die abgedeckte Fläche oder die Summe der paarweisen Abstände).

Hier ist ein Beispielcode, der das tut, aber es ist wirklich langsam. Ich brauche einen effizienteren Code, da ich ihn mehrmals aufrufen muss.

%% generating sample data 
n_4 = 1000; n_2 = n_4*2;n = n_4*4; 
x1=[ randn(n_4, 1)*10+30; randn(n_4, 1)*3 + 60]; 
y1=[ randn(n_4, 1)*5 + 35; randn(n_4, 1)*20 + 80 ]; 
x2 = rand(n_2, 1)*(max(x1)-min(x1)) + min(x1); 
y2 = rand(n_2, 1)*(max(y1)-min(y1)) + min(y1); 
P = [x1,y1;x2, y2]; 
%% eliminating close ones 
tic 
d0 = 1.5; 
D = pdist2(P, P);D(1:n+1:end) = inf; 
E = zeros(n, 1); % eliminated ones 
for i=1:n-1 
    if ~E(i) 
     CloseOnes = (D(i,:)<d0) & ((1:n)>i) & (~E'); 
     E(CloseOnes) = 1; 
    end 
end 
P2 = P(~E, :); 
toc 
%% plotting samples 
subplot(121); scatter(P(:, 1), P(:, 2)); axis equal; 
subplot(122); scatter(P2(:, 1), P2(:, 2)); axis equal; 

Edit: Wie gross ist die Teilmenge enthalten sein sollte?

Wie j_random_hacker in Kommentaren hingewiesen, kann man sagen, dass P(1, :) die schnellste Antwort ist, wenn wir keine Einschränkung für die Anzahl der ausgewählten Proben definieren. Es zeigt zierlich Inkohärenz des Titels! Aber ich denke, der aktuelle Titel beschreibt den Zweck besser. Lassen Sie uns eine Einschränkung definieren: "Versuchen Sie, m Beispiele auszuwählen, wenn es möglich ist". Mit der impliziten Annahme von m=n können wir nun die größtmögliche Teilmenge erhalten. Wie ich bereits erwähnt habe, übertrifft eine schnellere Methode diejenige, die die optimale Antwort findet.

Answer 1

Wenn Sie die nächsten Punkte immer wieder finden, wird eine andere Datenstruktur vorgeschlagen, die für die räumliche Suche optimiert ist. Ich schlage eine Delaunay-Triangulation vor.

Die folgende Lösung ist "ungefähr" in dem Sinne, dass wahrscheinlich mehr Punkte als unbedingt notwendig entfernt werden. Ich binde alle Berechnungen auf und entferne alle Punkte in jeder Iteration, die zu zu langen Abständen beitragen, und in vielen Fällen kann das Entfernen der Kante, die später in der gleichen Iteration erscheint, in vielen Fällen entfernt werden. Wenn dies von Bedeutung ist, kann die Kantenliste weiter verarbeitet werden, um Duplikate zu vermeiden, oder sogar um zu entfernende Punkte zu finden, die sich auf die größte Anzahl von Entfernungen auswirken.

Das ist schnell.

dt = delaunayTriangulation(P(:,1), P(:,2)); 
d0 = 1.5; 

while 1 
    edge = edges(dt); % vertex ids in pairs 

    % Lookup the actual locations of each point and reorganize 
    pwise = reshape(dt.Points(edge.', :), 2, size(edge,1), 2); 
    % Compute length of each edge 
    difference = pwise(1,:,:) - pwise(2,:,:); 
    edge_lengths = sqrt(difference(1,:,1).^2 + difference(1,:,2).^2); 

    % Find edges less than minimum length 
    idx = find(edge_lengths < d0); 
    if(isempty(idx)) 
     break; 
    end 

    % pick first vertex of each too-short edge for deletion 
    % This could be smarter to avoid overdeleting 
    points_to_delete = unique(edge(idx, 1)); 

    % remove them. triangulation auto-updates 
    dt.Points(points_to_delete, :) = []; 

    % repeat until no edge is too short 
end 

P2 = dt.Points; 

Answer 2

Sie geben nicht an, wie viele Punkte Sie auswählen möchten. Dies ist entscheidend für das Problem.

Ich sehe nicht leicht einen Weg, um Ihre Methode zu optimieren.

Unter der Annahme, dass die euklidische Entfernung als Abstandsmaß akzeptabel ist, ist die folgende Implementierung viel schneller bei der Auswahl nur einer kleinen Anzahl von Punkten und schneller sogar beim Versuch der Teilmenge mit 'allen' gültigen Punkten (beachten Sie das Maximum finden mögliche Anzahl von Punkten ist schwer).

%% 
figure; 
subplot(121); scatter(P(:, 1), P(:, 2)); axis equal; 

d0 = 1.5; 

m_range = linspace(1, 2000, 100); 
m_time = NaN(size(m_range)); 

for m_i = 1:length(m_range); 
    m = m_range(m_i) 

    a = tic; 
    % Test points in random order. 
    r = randperm(n); 
    r_i = 1; 

    S = false(n, 1); % selected ones 
    for i=1:m 
     found = false; 

     while ~found 
      j = r(r_i); 
      r_i = r_i + 1; 
      if r_i > n 
       % We have tried all points. Nothing else can be valid. 
       break; 
      end 
      if sum(S) == 0 
       % This is the first point. 
       found = true; 
      else 
       % Get the points already selected 
       P_selected = P(S, :); 
       % Exclude points >= d0 along either axis - they cannot have 
       % a Euclidean distance less than d0. 
       P_valid = (abs(P_selected(:, 1) - P(j, 1)) < d0) & (abs(P_selected(:, 2) - P(j, 2)) < d0); 
       if sum(P_valid) == 0 
        % There are no points that can be < d0. 
        found = true; 
       else 
        % Implement Euclidean distance explicitly rather than 
        % using pdist - this makes a large difference to 
        % timing. 
        found = min(sqrt(sum((P_selected(P_valid, :) - repmat(P(j, :), sum(P_valid), 1)) .^ 2, 2))) >= d0; 
       end 
      end 
     end 
     if found 
      % We found a valid point - select it. 
      S(j) = true; 
     else 
      % Nothing found, so we must have exhausted all points. 
      break; 
     end 
    end 
    P2 = P(S, :); 
    m_time(m_i) = toc(a); 
    subplot(122); scatter(P2(:, 1), P2(:, 2)); axis equal; 
    drawnow; 
end 
%% 
figure 
plot(m_range, m_time); 
hold on; 
plot(m_range([1 end]), ones(2, 1) * original_time); 
hold off; 

wobei original_time die Zeit ist, die Ihre Methode einnimmt. Dies ergibt die folgenden Zeiten, wobei die rote Linie Ihre Methode ist und das Blau meiner ist, mit der Anzahl der Punkte, die entlang der x-Achse ausgewählt sind. Beachten Sie, dass die Linie flacher wird, wenn "alle" Punkte ausgewählt wurden, die den Kriterien entsprechen.

Wie Sie in Ihrem Kommentar sagen, Leistung ist stark abhängig von dem Wert von d0. In der Tat, wie d0 reduziert wird, wird das Verfahren über eine noch größere Verbesserung der Leistung zu haben (dies ist für d0=0.1):

ist jedoch zu beachten, dass diese von anderen Faktoren auch abhängig ist wie die Verteilung von deine Daten. Diese Methode nutzt spezifische Eigenschaften Ihres Datensatzes und reduziert die Anzahl teurer Berechnungen, indem Punkte herausgefiltert werden, bei denen die Berechnung der euklidischen Entfernung sinnlos ist. Dies funktioniert besonders gut für die Auswahl weniger Punkte, und es ist tatsächlich schneller für kleinere d0, weil es weniger Punkte im Datensatz gibt, die den Kriterien entsprechen (so sind weniger Berechnungen der euklidischen Distanz erforderlich). Die optimale Lösung für ein Problem wie diese ist normalerweise spezifisch für den verwendeten Datensatz.

Beachten Sie auch, dass in meinem obigen Code die manuelle Berechnung der euklidischen Distanz viel schneller ist als das Aufrufen von pdist. Die Flexibilität und Allgemeingültigkeit der Matlab-Einbauten wirkt sich in einfachen Fällen oft nachteilig auf die Leistung aus.