# k-nearest neighbors using MATLAB with MEX

I have implemented kNN (k-nearest neighbors) as follows, but it is very slow. I want to get an exact k-nearest-neighbor, not the approximate ones, so I didn't use the FLANN or ANN libraries.

mexFindNN.cpp

#include <iostream>
using namespace std;

#include "mex.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <string.h>
#include <vector>
#include <algorithm>

struct Pair{
int id;
double value;
Pair(int id, double value){
this->id=id;
this->value=value;
}
};

struct PairCompare {
bool operator()(Pair const &left, Pair const &right) {
return left.value < right.value;
}
};

template<typename T>
void FindNN(T *X, T *Y, int N, int d, int type, int inner_k, int outer_k, mxArray *innerM, mxArray *outerM){

if(!type)//just inner_k
{
vector<size_t> ir;
vector<size_t> jc;jc.push_back(0);
vector<double> pr;
size_t num_ele=0;

for(int i=0;i<N;i++){//X[j*N+i]
vector<Pair> inner;
for(int j=0;j<N;j++){
double temp=0.0;
for(int k=0;k<d;k++){
temp+=(X[k*N+i]-X[k*N+j])*(X[k*N+i]-X[k*N+j]);
}
if(Y[i]==Y[j]){
inner.push_back(Pair(j,sqrt(temp)));
}
}

std::sort(inner.begin(),inner.end(),PairCompare());
for(int j=1;j<=inner_k && j<inner.size();j++){
Pair x=inner[j];
ir.push_back(x.id);
pr.push_back(x.value);
num_ele++;
}
jc.push_back(num_ele);
}

size_t *pIr=(size_t *)mxGetIr(innerM);
size_t *pJc=(size_t *)mxGetJc(innerM);
double *pPr=(double *)mxGetPr(innerM);
memcpy(pIr,&ir[0],ir.size()*sizeof(size_t));
memcpy(pJc,&jc[0],jc.size()*sizeof(size_t));
memcpy(pPr,&pr[0],pr.size()*sizeof(double));
}
else
{
vector<size_t> ir,ir2;
vector<size_t> jc,jc2;jc.push_back(0);jc2.push_back(0);
vector<double> pr,pr2;
size_t num_ele=0;
size_t num_ele2=0;

for(int i=0;i<N;i++){//X[j*N+i]
vector<Pair> inner, outer;
for(int j=0;j<N;j++){
double temp=0.0;
for(int k=0;k<d;k++){
temp+=(X[k*N+i]-X[k*N+j])*(X[k*N+i]-X[k*N+j]);
}
if(Y[i]==Y[j]){
inner.push_back(Pair(j,sqrt(temp)));
}else{
outer.push_back(Pair(j,sqrt(temp)));
}
}

std::sort(inner.begin(),inner.end(),PairCompare());
std::sort(outer.begin(),outer.end(),PairCompare());
for(int j=1;j<=inner_k && j<inner.size();j++){
Pair x=inner[j];
ir.push_back(x.id);
pr.push_back(x.value);
num_ele++;
}
jc.push_back(num_ele);

for(int j=0;j<outer_k && j<outer.size();j++){
Pair x=outer[j];
ir2.push_back(x.id);
pr2.push_back(x.value);
num_ele2++;
}
jc2.push_back(num_ele2);
}

size_t *pIr=(size_t *)mxGetIr(innerM);
size_t *pJc=(size_t *)mxGetJc(innerM);
double *pPr=(double *)mxGetPr(innerM);
memcpy(pIr,&ir[0],ir.size()*sizeof(size_t));
memcpy(pJc,&jc[0],jc.size()*sizeof(size_t));
memcpy(pPr,&pr[0],pr.size()*sizeof(double));

size_t *pIr2=(size_t *)mxGetIr(outerM);
size_t *pJc2=(size_t *)mxGetJc(outerM);
double *pPr2=(double *)mxGetPr(outerM);
memcpy(pIr2,&ir2[0],ir2.size()*sizeof(size_t));
memcpy(pJc2,&jc2[0],jc2.size()*sizeof(size_t));
memcpy(pPr2,&pr2[0],pr2.size()*sizeof(double));
}
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//prhs[0]: X
//prhs[1]: Y
//prhs[2]: inner_k
//prhs[3]: outer_k

//plhs[0]: inner_kNN_Matrix
//plhs[1]: outer_kNN_Matrix
//mwSize dims_n=mxGetNumberOfDimensions(prhs[0]);
const mwSize *dims= mxGetDimensions(prhs[0]);
int type=(int)mxGetScalar(prhs[2]);
int inner_k=(int)mxGetScalar(prhs[3]);
int outer_k=(int)mxGetScalar(prhs[4]);

mxClassID clsID = mxGetClassID(prhs[0]);
if(clsID==mxSINGLE_CLASS){
int N=dims[0];
int d=dims[1];
float *X=(float *)mxGetPr(prhs[0]);
float *Y=(float *)mxGetPr(prhs[1]);
plhs[0]=mxCreateSparse(N,N,N*inner_k,mxREAL);
if(type)
{
plhs[1]=mxCreateSparse(N,N,N*outer_k,mxREAL);
FindNN<float>(X,Y,N,d,type,inner_k,outer_k,plhs[0],plhs[1]);
}
else
{
FindNN<float>(X,Y,N,d,type,inner_k,outer_k,plhs[0],NULL);
}
}else if(clsID==mxDOUBLE_CLASS){
int N=dims[0];
int d=dims[1];
double *X=(double *)mxGetPr(prhs[0]);
double *Y=(double *)mxGetPr(prhs[1]);
plhs[0]=mxCreateSparse(N,N,N*inner_k,mxREAL);
if(type)
{
plhs[1]=mxCreateSparse(N,N,N*outer_k,mxREAL);
FindNN<double>(X,Y,N,d,type,inner_k,outer_k,plhs[0],plhs[1]);
}
else
{
FindNN<double>(X,Y,N,d,type,inner_k,outer_k,plhs[0],NULL);
}
}
}


ConstructNNGraph2.m

function [innerG,outerG]=ConstructNNGraph2(X,Y,inner_k,outer_k)
[N,d]=size(X);
if isempty(Y)
Y=ones(N,1);
end
type=0;
if outer_k>0
type=1;
end
if(type)
[innerG,outerG]=mexFindNN(X,Y,1,inner_k,outer_k);
innerG = max(innerG, innerG');
outerG = max(outerG, outerG');
else
[innerG]=mexFindNN(X,Y,0,inner_k,0);
outerG=[];
end


The code above needs to be compiled in a MATLAB environment. The compile command is

mex -largeArrayDims mexFindNN.cpp


The sample input X and Y is as follows:

load fisheriris;
Y=zeros(150,1);
Y(1:50)=1;
Y(51:100)=2;
Y(101:end)=3;
X=meas;
[innerG,outerG]=ConstructNNGraph2(X,Y,3,5);


Use a run-time profiler (test and measure, instead of or as well as trying to guess what's slow). Who knows, just from looking at it: the most expensive line of code might be something innocuous-looking like your push_back method calls.

For a description of how to do this, see for example Profiling a mex-function.

Code review

The following apply to the small code fragment posted in the original version of this question:

• std::sort followed by for(int j=1;j<=k...) isn't the cheapest way to get the k smallest elements in a vector. Instead, std::nth_element has linear cost.

• It would be better to reserve a capacity for knn_samples, otherwise its doing (expensive) heap allocations and reallocations when you push_back. You could define it once outside your outer loop, and empty it (in order to reuse it) at the top of each loop.

The above is an inadequate/incomplete review (there was much more code added the new version of the question) but I don't have time to add to it now.

Algorithm review

• In the question you said, " I want to get an exact k-nearest-neighbor, not the approximate ones, so I didn't use the FLANN or ANN libraries"
• In his answer, @miniBill said, "Your algorithm is O(n^2), and as much as you can optimize, you can't do better with this."

Here's an idea for improving the algorithm (I don't know whether this idea might be helpful):

• Use one of the fast, "approximate" libraries to categorize your data set into zones
• Use your expensive, exact, O(n^2) algorithm on data already partitioned into much smaller zones.

For example, imagine that you need to run this algorithm against all stars in the universe. An O(n^2) algorithm would do that slowly. I think it would be faster if you used an inexact algorithm to partition the stars into galaxies, and then run your exact algorithm on the stars within each galaxy.

There's no need to get exact values between two stars in different galaxies: an approximate result is enough to tell you that this pair would not be nearest neighbours.

• Thanks, sir. I have updated the question. Maybe it is now easy to review. In the code, I want to generate k-nearest-neighbor matrix of a data matrix. – mining Jan 25 '14 at 1:51
• Welcome to CodeReview. I updated my answer to suggest that perhaps using a combination of algorithms might make it faster. – ChrisW Jan 25 '14 at 20:19
• thanks, sir. I think I still should learn the FLANN/ANN library, especially for large scale data set. – mining Jan 26 '14 at 0:11

Your code appears to be very C-like with some C++. I'll just give some feedback in regards to that:

• Try not to use using namespace std.

• <string.h> is a C library; use <string> with C++.

• In C++, prefer std::size_t over size_t from C.

• For testing algorithms such as these, it's good to provide your main() to show how you're doing your testing. Although the code there may already work, you cannot always determine how well an algorithm works if there is no test code.

• Instead of creating a new Pair structure, consider using std::pair from the STL. It is more idiomatic C++, and it already comes with a few functions and operator overloads.

Moreover, Pair is not a very descriptive name in relation to this program. All that's known is that it holds an int and a double.

Here's how you can change this with std::pair:

// this creates an alias for a new std::pair type
// this is just a generic type name for demonstration
typedef std::pair<int, double> SomePair;

// create a new std::pair
SomePair newPair;

// pass it to a function
void someFunc(SomePair pair /* ... */) {}

• For better readability, keep operators and operands separated by whitespace:

for (int i = 0; i < 10; ++i) {}

• In this function:

void FindNN(T *X, T *Y, int N, int d, int k)


It's not clear what these variables are for as they're single-character. An exception to this is loop counters, which can be single-character.

• The OP's question was how to make the code faster. – ChrisW Jan 24 '14 at 22:07
• @ChrisW: I know, but on Code Review, we could review any aspects of the code we'd like. I also prefer for there to be something else for other reviewers to address, as that means more answers. – Jamal Jan 24 '14 at 22:08
• @Jamal, Thanks. I'm not very familiar with the C++ syntax, so the code looks like C. – mining Jan 25 '14 at 1:45
• @NicolasZhong: That's understandable. You did use std::vector and std::sort, which is a start. When you get more familiar with the language, you'll see how useful these tools can be. – Jamal Jan 25 '14 at 1:47

Your algorithm is O(n^2), and as much as you can optimize, you can't do better with this.

I don't see any particularly slow path.

• Thanks, sir. I think it will be better to use paralleled code, but I'm not very familiar with that. – mining Jan 25 '14 at 1:52