Schnellste Methode zur Berechnung des euklidischen Mindestabstands zwischen zwei Matrizen mit hochdimensionalen Vektoren

Ich fing eine ähnliche Frage anein anderer Thread, aber dann habe ich mich darauf konzentriert, wie man OpenCV benutzt. Nachdem ich nicht das erreicht habe, was ich ursprünglich wollte, werde ich hier genau fragen, was ich will.

Ich habe zwei Matrizen. Matrix a ist 2782x128 und Matrix b ist 4000x128, beides vorzeichenlose Zeichenwerte. Die Werte werden in einem einzelnen Array gespeichert. Für jeden Vektor in a benötige ich den Index des Vektors in b mit dem nächsten euklidischen Abstand.

Ok, jetzt mein Code, um dies zu erreichen:

#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;
    }
}

Und anbei die Dateien mit Beispielmatrizen.

matrixa Matrixb

Ich benutze windows.h nur, um den Zeitaufwand zu berechnen. Wenn Sie also den Code auf einer anderen Plattform als Windows testen möchten, ändern Sie einfach die Überschrift von windows.h und die Art und Weise, wie der Zeitaufwand berechnet wird.

Dieser Code in meinem Computer ist ungefähr 0,5 Sekunden. Das Problem ist, dass ich in Matlab einen anderen Code habe, der dasselbe in 0,05 Sekunden macht. In meinen Experimenten erhalte ich pro Sekunde mehrere Matrizen wie Matrix, also sind 0,5 Sekunden zu viel.

Nun der Matlab-Code, um dies zu berechnen:

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);

OK. Matlab-Code verwendet (x-a) ^ 2 = x ^ 2 + a ^ 2 - 2ab.

Also war mein nächster Versuch, dasselbe zu tun. Ich habe meinen eigenen Code gelöscht, um dieselben Berechnungen durchzuführen, aber es waren ca. 1,2 Sekunden.

Dann habe ich versucht, verschiedene externe Bibliotheken zu verwenden. Der erste Versuch war 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]);
}

Dieser Eigencode kostet nur für die Zeile, die besagt: ab = a * b.transpose ();

Ein ähnlicher Code mit opencv wurde ebenfalls verwendet, und die Kosten für die ab = a * b.transpose (); betrug 0,65 Sekunden.

Es ist wirklich ärgerlich, dass Matlab das auch so schnell kann und ich nicht in C ++ bin! Natürlich wäre es großartig, wenn ich mein Experiment durchführen könnte, aber ich denke, der Mangel an Wissen ist das, was mich wirklich nervt. Wie kann ich mindestens die gleiche Leistung erzielen wie in Matlab? Jede Art von Lösung ist willkommen. Ich meine, jede externe Bibliothek (kostenlos, wenn möglich), Loop-Unrolling-Dinge, Template-Dinge, SSE-Anweisungen (ich weiß, dass sie existieren), Cache-Dinge. Wie ich bereits sagte, besteht mein Hauptzweck darin, mein Wissen zu erweitern, damit ich in der Lage bin, solche Überlegungen mit einer schnelleren Leistung zu codieren.

Danke im Voraus

EDIT: mehr Code von David Hammen vorgeschlagen. Ich habe die Arrays auf int gesetzt, bevor ich irgendwelche Berechnungen anstellte. Hier ist der Code:

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;
    }
}

Der gesamte Prozess ist jetzt 0,6 und die Casting-Loops am Anfang sind 0,001 Sekunden. Vielleicht habe ich etwas falsch gemacht?

EDIT2: Irgendwas über Eigen? Wenn ich nach externen Bibliotheken suche, sprechen sie immer über Eigen und ihre Geschwindigkeit. Ich habe etwas falsch gemacht? Hier ein einfacher Code mit Eigen, der zeigt, dass er nicht so schnell ist. Vielleicht fehlt mir eine Konfiguration oder eine Flagge, oder ...

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

Dieser Code ist etwa 0,9 Sekunden.

Antworten auf die Frage(3)

Ihre Antwort auf die Frage