Accelerazione della distanza L1 tra tutte le coppie in un set di terra

Question

Ho una matrice NxM (di solito 10k X 10k elementi) che descrive un set di terra. Ogni linea rappresenta un oggetto e ogni colonna è una caratteristica specifica. Ad esempio, nella matrice

f1 f2 f3 
x1 0 4 -1 
x2 1 0 5 
x3 4 0 0 
x4 0 1 0 

dell'oggetto x1 ha valore 0 in funzione 1, valore 4 nella caratteristica 1, e il valore 0 in funzione -1. I valori di questo sono numeri reali generali (double).

Devo calcolare diverse distanze/dissomiglianze personalizzate tra tutte le coppie di oggetti (tutte le coppie di linee). Per confrontare, voglio calcolare le distanze L1 (manhattan) e L2 (euclidee).

Ho utilizzato la libreria Eigen per eseguire la maggior parte dei miei calcoli. Per calcolare la L2 (euclideo), utilizzo la seguente osservazione: per due vettori un e b di dimensioni n, abbiamo:

 
||a - b||^2 = (a_1 - b_1)^2 + (a_2 - b_2)^2 + ... +(a_n - b_n)^2 
      = a_1^2 + b_1^2 - 2 a_1 b_1 + a_2^2 + b_2^2 - 2 a_2 b_2 + ... + a_n^2 + b_n^2 - 2 a_n b_n 
      = a . a + b . b - 2ab 

In altre parole, si riscrive la norma quadrata usando il prodotto punto del vettore da solo e sottrarre il doppio del prodotto punto tra di loro. Da quello, prendiamo il quadrato e abbiamo finito. Col tempo, ho trovato questo trucco molto tempo fa e sfortunatamente ho perso il riferimento dell'autore.

Comunque, questo permette di scrivere un codice di fantasia utilizzando Eigen (in C++):

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> matrix, XX, D; 

// Load matrix here, for example 
// matrix << 0, 4, -1, 
//   1, 0, 5, 
//   4, 0, 0, 
//   0, 1, 0; 

const auto N = matrix.rows(); 

XX.resize(N, 1); 
D.resize(N, N); 

XX = matrix.array().square().rowwise().sum(); 

D.noalias() = XX * Eigen::MatrixXd::Ones(1, N) + 
       Eigen::MatrixXd::Ones(N, 1) * XX.transpose(); 

D -= 2 * matrix * matrix.transpose(); 
D = D.cwiseSqrt(); 

Per matrici 10k X 10k, siamo in grado di calcolare la distanza L2 per tutte le coppie di oggetti/linee in meno di 1 min (2 core/4 thread), che personalmente considero una buona prestazione per i miei scopi. Eigen è in grado di combinare le operazioni e utilizzare diverse ottimizzazioni di livello basso/alto per eseguire questi calcoli. In questo caso, Eigen utilizza tutti i core disponibili (e, ovviamente, possiamo regolarli).

Tuttavia, ho ancora bisogno di calcolare la distanza L1, ma non sono riuscito a capire una buona forma algebrica da utilizzare con Eigen e ottenere prestazioni soddisfacenti. Fino ad ora ho il seguente:

const auto N = matrix.rows(); 
for(long i = 0; i < N - 1; ++i) { 
    const auto &row = matrix.row(i); 

    #ifdef _OPENMP 
    #pragma omp parallel for shared(row) 
    #endif 
    for(long j = i + 1; j < N; ++j) { 
     distance(i, j) = (row - matrix.row(j)).lpNorm<1>(); 
    } 
} 

Ma questo prendere tempo molto più lungo: per la stessa matrice 10k X 10k, questo codice utilizza 3,5 min, che è molto peggio se si considera che L1 e L2 sono molto vicini nel loro originale forma:

L1(a, b) = sum_i |a_i - b_i| 
L2(a, b) = sqrt(sum_i |a_i - b_i|^2) 

Qualsiasi idea di come il cambiamento L1 utilizzare calcoli veloci con Eigen? O una forma migliore per farlo e io non l'ho capito.

Grazie mille per il vostro aiuto!

Carlos

Answer 1

lascia solo calcolare entrambe le distanze allo stesso tempo. Condividono davvero solo la differenza di riga (mentre entrambi potrebbero essere la differenza assoluta, la distanza euclidea usa una casella che non è realmente equivalente). Quindi ora stiamo solo scorrendo le righe n^2.

const auto N = matrix.rows(); 
for(long i = 0; i < N - 1; ++i) { 
    const auto &row = matrix.row(i); 

    #ifdef _OPENMP 
    #pragma omp parallel for shared(row) 
    #endif 
    for(long j = i + 1; j < N; ++j) { 
     const auto &rowDiff = row - matrix.row(j); 
     distanceL1(i, j) = rowDiff.cwiseAbs().sum(); // or .lpNorm<1>(); if it's faster 
     distanceL2(i, j) = rowDiff.norm() 
    } 
} 

EDIT un'altra più memoria di metodo/testata potrebbe calcolare una riga intera distanza ogni iterazione (non so se questo sarebbe un miglioramento o meno)

const auto N = matrix.rows(); 
#ifdef _OPENMP 
#pragma omp parallel for shared(matrix) 
#endif 
for(long i = 0; i < N - 1; ++i) { 
    const auto &row = matrix.row(i); 
    // you could use matrix.block(i,j,k,l) to cut down on the number of unnecessary operations 
    const auto &mat = matrix.rowwise() - row; 

    distanceL1(i) = mat.cwiseAbs().sum().transpose(); 
    distanceL2(i) = mat.rowwise().norm().transpose(); 
} 

Answer 2

Sono due operazioni molto comuni nell'elaborazione delle immagini. Il primo è il Sum of Squared Differences (SSD), il secondo è il Sum of Absolute Differences (SAD).

Hai determinato correttamente che SSD ne richiede solo uno per calcolare lo cross-correlation tra le due serie come il calcolo principale. Tuttavia, si può prendere in considerazione l'utilizzo della FFT per calcolare questi termini a.b, ridurrà il numero di operazioni necessarie significativamente per il caso L2 (tuttavia quanto non lo so, dipende da quale algoritmo di moltiplicazione matrice-matrice Eigen usi.) Se hai bisogno di me per spiegare questo posso, but I figure you can also look it up as its a standard use of FFTs. OpenCV ha un'implementazione (piuttosto scadente/errata) di corrispondenza dei modelli, che è ciò che si desidera quando si utilizza CV_TM_SQDIFF.

L1 caso è più complicato. Il caso L1 non può essere scomposto in modo carino, ma è anche una delle operazioni più semplici che puoi fare (solo aggiunte & valori assoluti.) Come tale, a lot of computation architectures have parallelized implementations di questo come istruzioni o funzioni implementate hardware. Altre architetture hanno researchers experimenting with the best way to compute this.

Si potrebbe anche voler guardare in Intel Imaging Primitives per il cross-correlazione, e le librerie veloce FFT come FFTW & CUFFT. Se non puoi permetterti Intel Imaging Primitves, puoi usare lo SSE instructions per accelerare notevolmente l'elaborazione quasi allo stesso effetto.