Najszybszy sposób obliczenia minimalnej odległości euklidesowej między dwiema macierzami zawierającymi wektory wielowymiarowe

Zacząłem podobne pytanieinny wątek, ale potem skupiłem się na tym, jak korzystać z OpenCV. Nie udało mi się osiągnąć tego, co pierwotnie chciałem, zapytam tutaj dokładnie, czego chcę.

Mam dwie macierze. Matryca a to 2782x128, a Matrix b to 4000x128, obie niepodpisane wartości char. Wartości są przechowywane w pojedynczej tablicy. Dla każdego wektora w a potrzebuję indeksu wektora w b z najbliższą odległością euklidesową.

Ok, teraz mój kod, aby to osiągnąć:

#include <windows.h>
#include <stdlib.h>
#include <stdio.h>
#include <cstdio>
#include <math.h>
#include <time.h>
#include <sys/timeb.h>
#include <iostream>
#include <fstream>
#include "main.h"

using namespace std;

void main(int argc, char* argv[])
{
    int a_size;
    unsigned char* a = NULL;
    read_matrix(&a, a_size,"matrixa");
    int b_size;
    unsigned char* b = NULL;
    read_matrix(&b, b_size,"matrixb");

    LARGE_INTEGER liStart;
    LARGE_INTEGER liEnd;
    LARGE_INTEGER liPerfFreq;
    QueryPerformanceFrequency( &liPerfFreq );
    QueryPerformanceCounter( &liStart );

    int* indexes = NULL;
    min_distance_loop(&indexes, b, b_size, a, a_size);

    QueryPerformanceCounter( &liEnd );

    cout << "loop time: " << (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) << "s." << endl;

    if (a)
    delete[]a;
if (b)
    delete[]b;
if (indexes)
    delete[]indexes;
    return;
}

void read_matrix(unsigned char** matrix, int& matrix_size, char* matrixPath)
{
    ofstream myfile;
    float f;
    FILE * pFile;
    pFile = fopen (matrixPath,"r");
    fscanf (pFile, "%d", &matrix_size);
    *matrix = new unsigned char[matrix_size*128];

    for (int i=0; i<matrix_size*128; ++i)
    {
        unsigned int matPtr;
        fscanf (pFile, "%u", &matPtr);
        matrix[i]=(unsigned char)matPtr;
    }
    fclose (pFile);
}

void min_distance_loop(int** indexes, unsigned char* b, int b_size, unsigned char* a, int a_size)
{
    const int descrSize = 128;

    *indexes = (int*)malloc(a_size*sizeof(int));
    int dataIndex=0;
    int vocIndex=0;
    int min_distance;
    int distance;
    int multiply;

    unsigned char* dataPtr;
    unsigned char* vocPtr;
    for (int i=0; i<a_size; ++i)
    {
        min_distance = LONG_MAX;
        for (int j=0; j<b_size; ++j)
        {
            distance=0;
            dataPtr = &a[dataIndex];
            vocPtr = &b[vocIndex];

            for (int k=0; k<descrSize; ++k)
            {
                multiply = *dataPtr++-*vocPtr++;
                distance += multiply*multiply;
                // If the distance is greater than the previously calculated, exit
                if (distance>min_distance)
                    break;
            }

            // if distance smaller
            if (distance<min_distance)
            {
                min_distance = distance;
                (*indexes)[i] = j;
            }
            vocIndex+=descrSize;
        }
        dataIndex+=descrSize;
        vocIndex=0;
    }
}

I dołączone są pliki z przykładowymi macierzami.

matrixa matrixb

Używam windows.h tylko do obliczenia czasu zużycia, więc jeśli chcesz przetestować kod na innej platformie niż Windows, po prostu zmień nagłówek windows.h i zmień sposób obliczania czasu zużycia.

Ten kod w moim komputerze wynosi około 0,5 sekundy. Problem polega na tym, że mam inny kod w Matlab, który czyni to samo w 0,05 sekundy. W moich eksperymentach otrzymuję kilka macierzy, takich jak macierz, co sekundę, więc 0,5 sekundy to za dużo.

Teraz kod matlab do obliczenia tego:

aa=sum(a.*a,2); bb=sum(b.*b,2); ab=a*b'; 
d = sqrt(abs(repmat(aa,[1 size(bb,1)]) + repmat(bb',[size(aa,1) 1]) - 2*ab));
[minz index]=min(d,[],2);

Dobrze. Kod Matlaba używa tego (x-a) ^ 2 = x ^ 2 + a ^ 2 - 2ab.

Więc moją następną próbą było zrobić to samo. Usunąłem swój własny kod, aby wykonać te same obliczenia, ale to było około 1,2 sekundy.

Następnie próbowałem użyć różnych zewnętrznych bibliotek. Pierwsza próba to Eigen:

const int descrSize = 128;
MatrixXi a(a_size, descrSize);
MatrixXi b(b_size, descrSize);
MatrixXi ab(a_size, b_size);

unsigned char* dataPtr = matrixa;
for (int i=0; i<nframes; ++i)
{
    for (int j=0; j<descrSize; ++j)
    {
        a(i,j)=(int)*dataPtr++;
    }
}
unsigned char* vocPtr = matrixb;
for (int i=0; i<vocabulary_size; ++i)
{
    for (int j=0; j<descrSize; ++j)
    {
        b(i,j)=(int)*vocPtr ++;
    }
}
ab = a*b.transpose();
a.cwiseProduct(a);
b.cwiseProduct(b);
MatrixXi aa = a.rowwise().sum();
MatrixXi bb = b.rowwise().sum();
MatrixXi d = (aa.replicate(1,vocabulary_size) + bb.transpose().replicate(nframes,1) - 2*ab).cwiseAbs2();

int* index = NULL;
index = (int*)malloc(nframes*sizeof(int));
for (int i=0; i<nframes; ++i)
{
    d.row(i).minCoeff(&index[i]);
}

Ten kod Eigen kosztuje 1,2 ok. Tylko dla linii, która mówi: ab = a * b.transpose ();

Podobny kod wykorzystano także przy użyciu opencv, a koszt ab = a * b.transpose (); wynosił 0,65 sekundy.

To naprawdę denerwujące, że matlab jest w stanie zrobić to samo tak szybko i nie jestem w stanie w C ++! Oczywiście możliwość przeprowadzenia mojego eksperymentu byłaby świetna, ale myślę, że brak wiedzy jest tym, co naprawdę mnie denerwuje. Jak mogę osiągnąć co najmniej taką samą wydajność jak w Matlab? Każdy rodzaj rozwiązania jest mile widziany. Mam na myśli dowolną zewnętrzną bibliotekę (jeśli to możliwe bezpłatną), rozwijanie pętli, rzeczy szablonów, instrukcje SSE (wiem, że istnieją), cache rzeczy. Jak już powiedziałem, moim głównym celem jest zwiększenie wiedzy, ponieważ możliwość kodowania myśli w ten sposób z szybszą wydajnością.

Z góry dziękuję

EDYCJA: więcej kodu sugerowanego przez Davida Hammen. Przed wykonaniem jakichkolwiek obliczeń rzuciłem tablice na int. Oto kod:

void min_distance_loop(int** indexes, unsigned char* b, int b_size, unsigned char* a, int a_size)
{
    const int descrSize = 128;

    int* a_int;
    int* b_int;

    LARGE_INTEGER liStart;
    LARGE_INTEGER liEnd;
    LARGE_INTEGER liPerfFreq;
    QueryPerformanceFrequency( &liPerfFreq );
    QueryPerformanceCounter( &liStart );

    a_int = (int*)malloc(a_size*descrSize*sizeof(int));
    b_int = (int*)malloc(b_size*descrSize*sizeof(int));

    for(int i=0; i<descrSize*a_size; ++i)
        a_int[i]=(int)a[i];
    for(int i=0; i<descrSize*b_size; ++i)
        b_int[i]=(int)b[i];

    QueryPerformanceCounter( &liEnd );

    cout << "Casting time: " << (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) << "s." << endl;

    *indexes = (int*)malloc(a_size*sizeof(int));
    int dataIndex=0;
    int vocIndex=0;
    int min_distance;
    int distance;
    int multiply;

    /*unsigned char* dataPtr;
    unsigned char* vocPtr;*/
    int* dataPtr;
    int* vocPtr;
    for (int i=0; i<a_size; ++i)
    {
        min_distance = LONG_MAX;
        for (int j=0; j<b_size; ++j)
        {
            distance=0;
            dataPtr = &a_int[dataIndex];
            vocPtr = &b_int[vocIndex];

            for (int k=0; k<descrSize; ++k)
            {
                multiply = *dataPtr++-*vocPtr++;
                distance += multiply*multiply;
                // If the distance is greater than the previously calculated, exit
                if (distance>min_distance)
                    break;
            }

            // if distance smaller
            if (distance<min_distance)
            {
                min_distance = distance;
                (*indexes)[i] = j;
            }
            vocIndex+=descrSize;
        }
        dataIndex+=descrSize;
        vocIndex=0;
    }
}

Cały proces trwa teraz 0,6, a pętle rzutowe na początku to 0,001 sekundy. Może zrobiłem coś złego?

EDIT2: Coś o Eigen? Kiedy szukam zewnętrznych bibliotek, zawsze mówią o Eigenu i jego szybkości. Zrobiłem coś nie tak? Tutaj prosty kod wykorzystujący Eigen, który pokazuje, że nie jest tak szybki. Może brakuje mi konfiguracji lub jakiejś flagi albo ...

MatrixXd A = MatrixXd::Random(1000, 1000);
MatrixXd B = MatrixXd::Random(1000, 500);
MatrixXd X;

Ten kod to około 0,9 sekundy.

questionAnswers(3)

yourAnswerToTheQuestion