Efficient Normal Correlation between 2 vectors

[sgsawant]

float crocor(float a[], float b[]){
	float numerator=0;
	float denominator=0;
	float t1=0,t2=0,t3=0,t4=0,t5=0;
	for(int i=0; i<=SCANRAD; i++){
		t1+=a[i]*b[i];
		t2+=a[i];
		t3+=b[i];
		//flsq simply squares the floating point number
		t4+=flsq(a[i]);	
		t5+=flsq(b[i]);
	}
	numerator=flsq((SCANRAD+1)*t1-t2*t3);
	denominator=((SCANRAD+1)*t4-flsq(t2))*((SCANRAD+1)*t5-flsq(t3));
	if(denominator==0) return(0);
	else	return(numerator/denominator);
}

The above code implements the formula shown in the attached image.

This code runs a number of times in my final code and the speed of the program hinges almost solely on this. Memory is not an issue. Please suggest any optimization if you can.

Regards.

P.S.: The (x 100) in the given formula is for calculating percent correlation and is not needed in the above case.

[Lindley]

That's probably about as good as you're going to get at that level. Any significant optimizations will more likely be related to higher-level code controlling how often that gets called.

Doing the computation entirely in integer math might be faster, but will probably be more difficult. You're liable to run into overflows, etc.

The only thing I see which you might consider changing is not computing the numerator until after you check whether the denominator is 0, but that's going to be of minimal help I suspect.

Make sure flsq() is an inline function.

[sgsawant]

Thanks for the tip to make flsq() inline.

Hadn't thought about it!!!

[OReubens]

Make sure that particular function is being optimised for speed.

It's generally not beneficial to have ALL your code optimised for speed since 90% of it isn't time critical, so optimising the bulk of your program for size and the critical parts for speed should give you the overall best speed.

If SCANRAD is high, you should be able to benefit from doing the loop calculation on multiple threads. Since you mention this code is called often you can create the worker threads once, and reuse them each time you get to this bit (thread creation and thread shutdown is rather expensive afterall).

This function itself is indeed about as good as it's going to get, however if this function is called often, then there may be any number of tweaks that may provide a positive benefit.
If the 2 tables don't change a lot in between calls then it may be advantagous to change this to a class that keeps a running summation, this kind of summation would benefit most for large values of SCANRAD and only a few elements of the tables changing in between each call to CroCor.

Alternatively, this code may be a good candidate to benefit from going to MMX/SSE/SSE2, but that would mean non portable code, and code that may not work as well on varying processors.

If SCANRAD is small, then unfolding the loop might be a good way to go.

If you only have 2 tables in your entire program where you're doing this calculation on. Make the tables global instead of passing them. Yes, it's not really that 'structured', but when performance matters...

[OReubens]

Without seeing more of your code, there's not much we can help with. But a typical rule of optimisation is that if you have a function eating a lot of your processing time... Finding clever ways to not having to call that function is one of the possible ways to speeding things up.

[sgsawant]

Thanks for the threading idea. Right now I am comfortable with classes. I also find the concept of "template" manageable. Keywords like "inline" etc. is where it becomes a bit difficult for me.

I have said this to give a general idea as to where I stand in C++. So can anyone recommend a good link to start reading about the concept of "threads"?

Regards.

P.S.: The SCANRAD usually lies between 20~70. Varies with the test case. I thought about ways to come up with different ways not to call the function, but I failed. Unfortunately, posting any more of my code would be a direct violation of my team guidelines.

[sgsawant]

Alternatively, this code may be a good candidate to benefit from going to MMX/SSE/SSE2, but that would mean non portable code, and code that may not work as well on varying processors.

It's all Greek & Latin to me! :P

[Paul McKenzie]

It's all Greek & Latin to me! :P

In other words, take advantage of your processor's instruction set to produce faster code.

Regards,

Paul McKenzie

[Lindley]

It's worth checking whether LAPACK's correlation routines will suit your needs. Something like that which is available in a form optimized for specific processors in some cases is probably going to be at least as fast as anything you write by hand.

[OReubens]

Again... if the tables are fairly stable, and only change a little in between calls to the correlation calculation. Keeping a runnign total might make drastic improvements.

basically.
Make a class. which has a running total member for t1 to t5.
Make the 2 tabels member of the class.
Make accessor functions to set a value in either table.
make an function that returns the total using t1-t5.

The accessor to set the tables would be something like:

void Set_A(int index, float NewA)
{
  t1 += (NewA - a[index]) * b[index];  // Remove old value, add new value.
  t2 += NewA - a[index];
  t4 += flsq(NewA) - flsq(a[index]);
  a[index] = NewA;
}

Similar type function for Set_B(),

the correlation calculation then can use the available t1 to t5 instead of havign to recalculate it all.

The above will only work if the table gets a low amount of changes for each call to correlation calculation. If the entire table changes each time it's not going to work out.