Databases
StackOverflow2013
Postsdbo
POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
 

Note that there are some explanatory texts on larger screens.

plurals
  1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
    primarykey
    data
    text
    <p>I started a similar question on <a href="https://stackoverflow.com/questions/12479663/cvmat-cv-8u-product-error-and-slow-cv-32f-product/12527400#12527400">another thread</a>, but then I was focusing on how to use OpenCV. Having failed to achieve what I originally wanted, I will ask here exactly what I want.</p> <p>I have two matrices. Matrix a is 2782x128 and Matrix b is 4000x128, both unsigned char values. The values are stored in a single array. For each vector in a, I need the index of the vector in b with the closest euclidean distance.</p> <p>Ok, now my code to achieve this:</p> <pre><code>#include &lt;windows.h&gt; #include &lt;stdlib.h&gt; #include &lt;stdio.h&gt; #include &lt;cstdio&gt; #include &lt;math.h&gt; #include &lt;time.h&gt; #include &lt;sys/timeb.h&gt; #include &lt;iostream&gt; #include &lt;fstream&gt; #include "main.h" using namespace std; void main(int argc, char* argv[]) { int a_size; unsigned char* a = NULL; read_matrix(&amp;a, a_size,"matrixa"); int b_size; unsigned char* b = NULL; read_matrix(&amp;b, b_size,"matrixb"); LARGE_INTEGER liStart; LARGE_INTEGER liEnd; LARGE_INTEGER liPerfFreq; QueryPerformanceFrequency( &amp;liPerfFreq ); QueryPerformanceCounter( &amp;liStart ); int* indexes = NULL; min_distance_loop(&amp;indexes, b, b_size, a, a_size); QueryPerformanceCounter( &amp;liEnd ); cout &lt;&lt; "loop time: " &lt;&lt; (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) &lt;&lt; "s." &lt;&lt; endl; if (a) delete[]a; if (b) delete[]b; if (indexes) delete[]indexes; return; } void read_matrix(unsigned char** matrix, int&amp; matrix_size, char* matrixPath) { ofstream myfile; float f; FILE * pFile; pFile = fopen (matrixPath,"r"); fscanf (pFile, "%d", &amp;matrix_size); *matrix = new unsigned char[matrix_size*128]; for (int i=0; i&lt;matrix_size*128; ++i) { unsigned int matPtr; fscanf (pFile, "%u", &amp;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&lt;a_size; ++i) { min_distance = LONG_MAX; for (int j=0; j&lt;b_size; ++j) { distance=0; dataPtr = &amp;a[dataIndex]; vocPtr = &amp;b[vocIndex]; for (int k=0; k&lt;descrSize; ++k) { multiply = *dataPtr++-*vocPtr++; distance += multiply*multiply; // If the distance is greater than the previously calculated, exit if (distance&gt;min_distance) break; } // if distance smaller if (distance&lt;min_distance) { min_distance = distance; (*indexes)[i] = j; } vocIndex+=descrSize; } dataIndex+=descrSize; vocIndex=0; } } </code></pre> <p>And attached are the files with sample matrices.</p> <p><a href="https://dl.dropbox.com/u/1474325/matrixa" rel="nofollow noreferrer">matrixa</a> <a href="https://dl.dropbox.com/u/1474325/matrixb" rel="nofollow noreferrer">matrixb</a></p> <p>I am using windows.h just to calculate the consuming time, so if you want to test the code in another platform than windows, just change windows.h header and change the way of calculating the consuming time.</p> <p>This code in my computer is about 0.5 seconds. The problem is that I have another code in Matlab that makes this same thing in 0.05 seconds. In my experiments, I am receiving several matrices like matrix a every second, so 0.5 seconds is too much.</p> <p>Now the matlab code to calculate this:</p> <pre><code>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); </code></pre> <p>Ok. Matlab code is using that (x-a)^2 = x^2 + a^2 - 2ab.</p> <p>So my next attempt was to do the same thing. I deleted my own code to make the same calculations, but It was 1.2 seconds approx.</p> <p>Then, I tried to use different external libraries. The first attempt was Eigen:</p> <pre><code>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&lt;nframes; ++i) { for (int j=0; j&lt;descrSize; ++j) { a(i,j)=(int)*dataPtr++; } } unsigned char* vocPtr = matrixb; for (int i=0; i&lt;vocabulary_size; ++i) { for (int j=0; j&lt;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&lt;nframes; ++i) { d.row(i).minCoeff(&amp;index[i]); } </code></pre> <p>This Eigen code costs 1.2 approx for just the line that says: ab = a*b.transpose();</p> <p>A similar code using opencv was used also, and the cost of the ab = a*b.transpose(); was 0.65 seconds.</p> <p>So, It is real annoying that matlab is able to do this same thing so quickly and I am not able in C++! Of course being able to run my experiment would be great, but I think the lack of knowledge is what really is annoying me. How can I achieve at least the same performance than in Matlab? Any kind of soluting is welcome. I mean, any external library (free if possible), loop unrolling things, template things, SSE intructions (I know they exist), cache things. As I said, my main purpose is increase my knowledge for being able to code thinks like this with a faster performance.</p> <p>Thanks in advance</p> <p>EDIT: more code suggested by David Hammen. I casted the arrays to int before making any calculations. Here is the code:</p> <pre><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( &amp;liPerfFreq ); QueryPerformanceCounter( &amp;liStart ); a_int = (int*)malloc(a_size*descrSize*sizeof(int)); b_int = (int*)malloc(b_size*descrSize*sizeof(int)); for(int i=0; i&lt;descrSize*a_size; ++i) a_int[i]=(int)a[i]; for(int i=0; i&lt;descrSize*b_size; ++i) b_int[i]=(int)b[i]; QueryPerformanceCounter( &amp;liEnd ); cout &lt;&lt; "Casting time: " &lt;&lt; (liEnd.QuadPart - liStart.QuadPart) / long double(liPerfFreq.QuadPart) &lt;&lt; "s." &lt;&lt; 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&lt;a_size; ++i) { min_distance = LONG_MAX; for (int j=0; j&lt;b_size; ++j) { distance=0; dataPtr = &amp;a_int[dataIndex]; vocPtr = &amp;b_int[vocIndex]; for (int k=0; k&lt;descrSize; ++k) { multiply = *dataPtr++-*vocPtr++; distance += multiply*multiply; // If the distance is greater than the previously calculated, exit if (distance&gt;min_distance) break; } // if distance smaller if (distance&lt;min_distance) { min_distance = distance; (*indexes)[i] = j; } vocIndex+=descrSize; } dataIndex+=descrSize; vocIndex=0; } } </code></pre> <p>The entire process is now 0.6, and the casting loops at the beginning are 0.001 seconds. Maybe I did something wrong?</p> <p>EDIT2: Anything about Eigen? When I look for external libs they always talk about Eigen and their speed. I made something wrong? Here a simple code using Eigen that shows it is not so fast. Maybe I am missing some config or some flag, or ...</p> <pre><code>MatrixXd A = MatrixXd::Random(1000, 1000); MatrixXd B = MatrixXd::Random(1000, 500); MatrixXd X; </code></pre> <p>This code is about 0.9 seconds.</p>
    singulars
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. PTQuestion
    1. USCœur
    1. USmin.yong.yoon
    plurals
    1. PL
      singulars
      1. LTLinked
    2. PL
      singulars
      1. LTLinked
    1. PL
      singulars
      1. LTLinked
    1. This table or related slice is empty.
    1. PO
      singulars
      1. PTAnswer
    2. PO
      singulars
      1. PTAnswer
    3. PO
      singulars
      1. PTAnswer
    1. VO
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. USDaV
      1. VTFavorite
    3. VO
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. This table or related slice is empty.
      1. VTUpMod
    1. COYou compiled all tests in release mode?
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. USDenis Ermolin
    2. COYou may find it annoying that Matlab outperforms your code, but I, who make heavy use of Matlab, find it very satisfying. I don't have much in the way of concrete advice for you but the key to improving the performance of this type of code is often to make optimal (or at least very good) use of the memory hierarchy on modern CPUs. Another factor to consider is that much of Matlab's core functionality is now multi-threaded for execution on multi-core CPUs, it's not clear to me that any of your own code is multi-threaded; that would probably have some impact on performance.
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. USHigh Performance Mark
    3. COI don't know how to help you make your C/C++ code faster (your code looks more like C than C++. Evidence: `malloc`) yet, but I can see how you can make your Matlab code faster: Eliminate the `sqrt`. Given two non-negative numbers a and b, sqrt(a)>sqrt(b) ⇔ a>b.
      singulars
      1. POFastest way to calculate minimum euclidean distance between two matrices containing high dimensional vectors
      1. USDavid Hammen
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload