k-nearest neighbors using MATLAB with MEX

Question

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

Answer 1

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.

Answer 2

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.

Answer 3

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.