hello,

i have a routine to calculate the 3x3 covariance matrix of a nx3 matrix, however it's not very efficient as it uses 3 nested for loops, and with large data sets it really clogs up. i was wondering if anyone could offer any advice to speed it up at all, the code is below, cheers,

Code:

double** getcovariancemat(Object* o, Object* p, int ndims)
{
int i, j, k;
//malloc space for cov-matrix here
Vertex* omean = (Vertex*) malloc(sizeof(Vertex));
Vertex* pmean = mean(p);
//number of vertices in the data object (the subsampled object)
int pverts = p->nvertices;
Vertex* ann[pverts];
void* kdtree = buildkdtree(o);
//create an array of pverts size filled with nearest neighbors of
//vertices of p contained in o
for (i = 0; i < pverts; ++i) {
ann[i] = nearestneighbour(&(p->vertices[i]), (struct kdtree*)kdtree);
omean->x += (ann[i]->x);
omean->y += (ann[i]->y);
omean->z += (ann[i]->z);
//printf("x: %g, y: %g, z: %g, \n",ann[i]->x, ann[i]->y, ann[i]->z);
}
omean->x /= pverts;
omean->y /= pverts;
omean->z /= pverts;
//arrays for calculating the variances from the data sets
double modelvar[pverts][3];
double datavar[pverts][3];
for (i = 0; i < pverts; ++i) {
modelvar[i][0] = (p->vertices[i].x - pmean->x);
modelvar[i][1] = (p->vertices[i].y - pmean->y);
modelvar[i][2] = (p->vertices[i].z - pmean->z);
datavar[i][0] = (ann[i]->x - omean->x);
datavar[i][1] = (ann[i]->y - omean->y);
datavar[i][2] = (ann[i]->z - omean->z);
}
//calculating the covariance matrix - not very efficient O(n^3)
for (i=0; i< ndims; ++i) {
for (j = 0; j < ndims; ++j) {
for (k=0; k < pverts; ++k) {
double val = (modelvar[k][i]*datavar[k][j]);
m[i][j] += val;
}
m[i][j] /= (pverts);
}
}
for (i = 0; i < pverts; ++i) {
free(ann[i]);
}
free(pmean);
free(omean);
return m;
}