c++ - Acelerando la distancia L1 entre todos los pares en un conjunto de suelo
algorithm matrix (2)
Tengo una matriz NxM (generalmente 10k X 10k elementos) que describe un conjunto de suelo. Cada línea representa un objeto y cada columna una característica específica. Por ejemplo, en la matriz
f1 f2 f3
x1 0 4 -1
x2 1 0 5
x3 4 0 0
x4 0 1 0
el objeto x1 tiene el valor 0 en la característica 1, el valor 4 en la característica 1 y el valor 0 en la característica -1. Los valores de esto son números reales generales (dobles).
Tengo que calcular varias distancias / disimilitudes personalizadas entre todos los pares de objetos (todo el par de líneas). Para comparar, quiero calcular las distancias L1 (Manhattan) y L2 (euclídea).
He usado la biblioteca Eigen para realizar la mayor parte de mis cálculos. Para calcular el L2 (euclidiano), utilizo la siguiente observación: para dos vectores ayb de tamaño n , tenemos:
||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
En otras palabras, reescribimos la norma al cuadrado usando el producto de puntos del vector por sí mismos y restamos dos veces el producto de puntos entre ellos. De eso, simplemente tomamos el cuadrado y terminamos. Con el tiempo, encontré este truco hace mucho tiempo y desafortunadamente perdí la referencia para el autor.
De todos modos, esta habilitación para escribir un código de fantasía usando Eigen (en 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();
Para la matriz 10k X 10k, podemos calcular la distancia L2 para todos los pares de objetos / líneas en menos de 1 min (2 núcleos / 4 hilos), lo que personalmente considero un buen rendimiento para mis propósitos. Eigen puede combinar las operaciones y usar varias optimizaciones de bajo / alto nivel para realizar estos cálculos. En este caso, Eigen usa todos los núcleos disponibles (y, por supuesto, podemos sintonizar eso).
Sin embargo, todavía necesito calcular la distancia L1, pero no pude encontrar una buena forma algebraica para usar con Eigen y obtener un buen rendimiento. Hasta ahora tengo lo siguiente:
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>();
}
}
Pero esto lleva mucho más tiempo: para la misma matriz de 10k X 10k, este código usa 3.5 min, lo cual es mucho peor teniendo en cuenta que L1 y L2 están muy cerca en su forma original:
L1(a, b) = sum_i |a_i - b_i|
L2(a, b) = sqrt(sum_i |a_i - b_i|^2)
¿Alguna idea de cómo cambiar L1 para usar cálculos rápidos con Eigen? O una forma mejor de hacer eso y simplemente no me di cuenta.
¡Muchas gracias por su ayuda!
Carlos
Permite calcular ambas distancias al mismo tiempo. En realidad, solo comparten la diferencia de fila (mientras que ambos podrían ser una diferencia absoluta, la distancia euclidiana utiliza un cuadrado que no es realmente equivalente). Entonces ahora solo estamos recorriendo las filas 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()
}
}
EDITAR otro método más intensivo en memoria / no probado podría ser calcular una fila de distancia completa en cada iteración (no sé si esto sería una mejora o no)
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();
}
Estas son dos operaciones muy comunes en el procesamiento de imágenes. El primero es la suma de las diferencias cuadradas (SSD) , el segundo es la suma de las diferencias absolutas (SAD) .
Ha determinado correctamente que SSD solo requiere uno para calcular la correlación cruzada entre las dos series como el cómputo principal. Sin embargo, es posible que desee considerar el uso de la FFT para calcular estos términos ab , reducirá el número de operaciones necesarias de forma significativa para el caso L2 (sin embargo, no sé cuánto depende del algoritmo de multiplicación matriz-matriz que Eigen utiliza). .) Si necesitas que explique esto, puedo, pero creo que también puedes buscarlo ya que es un uso estándar de las FFT .OpenCV tiene una implementación (bastante mala / defectuosa) de coincidencia de plantilla, que es lo que quiere cuando usa CV_TM_SQDIFF.
El caso L1 es más complicado. La carcasa L1 no puede descomponerse tan bien, pero también es una de las operaciones más simples que puede hacer (solo adiciones y valores absolutos). Como tal, muchas arquitecturas de computación tienen implementaciones paralelizadas de esto como instrucciones o funciones implementadas por hardware. . Otras arquitecturas tienen a los investigadores experimentando con la mejor manera de calcular esto.
También es posible que desee buscar en Intel Imaging Primitives para la correlación cruzada, y bibliotecas de FFT rápidas como FFTW y CUFFT . Si no puede pagar Intel Image Primitves, puede usar las instrucciones SSE para acelerar su procesamiento casi con el mismo efecto.