[ExternalEmail] Pearson Correlation Speed
On Tue, 16 Dec 2008, Nathan S. Watson-Haigh wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Charles C. Berry wrote:
On Mon, 15 Dec 2008, Nathan S. Watson-Haigh wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Nathan S. Watson-Haigh wrote:
I'm trying to calculate Pearson correlation coefficients for a large matrix of size 18563 x 18563. The following function takes about XX minutes to complete, and I'd like to do this calculation about 15 times and so speed is some what of an issue.
I think you are on the wrong track, Nathan. The matrix you are starting with is 18563 x 18563 and the result of finding the correlations amongst the columns of that matrix is also 18563 x 18563. It will require more than 5 Gigabytes of memory to store the result and the original matrix.
Yes the memory usage is somewhat large - luckily I have the use of a cluster with lots of shared memory! However, I'm interested to learn how you came about the calculation to determine the memory requirements.
The original object is
18563^2*8/1024^3
[1] 2.567358
Gigabytes, and so is the result. I added them together.
Likely the time needed to do the calc is inflated because of caching issues and if your machine has less than enough memory to store the result and all the intermediate pieces by swapping as well. You can finesse these by breaking your problem into smaller pieces, say computing the correlations between each pair of 19 blocks of columns (columns 1:977, 977+1:977, ... 18*977+1:977 ), then assembling the results.
This is possibly, however why is something like this not implemented internally in the cor() function if it poorly scales due to the large memory requirements?
Because nobody ever really needed it? Seriously, optimizing something like this is machine dependent, and R-core probably has higher priorities. cor() provides lots of options - it handles NAs, for example - and it is probably not worth the trouble to try to optimize over those options. The calculation sans NAs is a simple one and can be done using the built in BLAS (as crossprod() does), which BLAS can in turn be tuned to the machine used. So, if your environment has a tuned or multithreaded BLAS, you might be better off to use crossprod() and scale the result.
--- BTW, R already has the necessary machinery to calculate the crossproduct matrix (etc) needed to find the correlations. You can access the low level linear algebra that R uses. You can marry R to an optimized BLAS if you like. So pulling in some other code to do this will not save you anything. If you ever do decide to import C[++] code there is excellent documentation in the Writing R Extensions manual, which you should review before attempting to import C++ code into R.
Thanks, I have seen this and it seemed quite technical to use as a starting point for someone unfamiliar with both C++ and incorporating C++ code into R.
Well, in that case the path of least resistance is to start the process
when you leave for the night and pick up the results the next morning.
HTH,
Chuck
Charles C. Berry (858) 534-2098
Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901