Error message when trying to compute EOF with covariance matrix.

I am trying to compute EOFs from a code I used before, but it does not seem to work here since my dataset has a too good spatial resolution (1/4 degree).

I want to compute the firsts 2 EOFs of mean June SLP northward of 70N for 43 years.

latitude has size 81 (4x20 = 80)

longitude has size 1440 (4x360 = 1440)

The data I want to compute the EOFs from has size 1440x81x43 (latitude x longitude x time)

When I want to compute the covariance matrix, I can't. I get the following error message:

Error using *

Requested 116640x116640 (101.4GB) array exceeds maximum array size preference (62.6GB). This might cause

MATLAB to become unresponsive.

Error in cov (line 155)

c = (xc' * xc) ./ denom;

Related documentation

%% Compute EOFs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%I CAN'T ATTACH MY DATA SINCE IT IS 31.1mb AND THE LIMIT IS 5 MB
load('mslpJune_dt.mat'); %load detrended data

mslpJune_dt_reshaped        = reshape(mslpJune_dt,[1440*81 43]);
mslpJune_dt_reshaped_trans  = transpose(mslpJune_dt_reshaped);

%I CAN'T COMPUTE COV, ERROR MESSAGE FROM MATLAB THAT IT IS TOO BIG
C                           = cov(mslpJune_dt_reshaped_trans);


[eigenvectors,eigenvalues]  = eigs(double(C),2);

%Principal components
pcs                         = mslpJune_dt_reshaped*eigenvectors;

%Scaling pcs
pcs1_norm(:,1)=pcs(:,1)/std(pcs(:,1));
pcs2_norm(:,2)=pcs(:,2)/std(pcs(:,2));

Does anyone know how to fix this issue and compute the EOFs?

Thank you

Answer 1

[Christine Tobler]

The problem is that the covariance matrix becomes very large here. Luckily, it's not necessary to compute this matrix explicitly: It's enough to compute a matrix xc such that C = xc' * xc. Then, instead of computing eig(C) it's equivalent to compute svd(xc). I explain a bit more in the code below:

mslpJune_dt_reshaped        = randn(1440*81, 43);
mslpJune_dt_reshaped_trans  = transpose(mslpJune_dt_reshaped);

% COV would return a matrix that's too large to store. Let's take the lines
% of code that COV does internally (see "edit cov") and see where that gets us:
x = mslpJune_dt_reshaped_trans;
[m,n] = size(x);
denom = m - 1;
xc = x - sum(x,1)./m;  % Remove mean
%c = (xc' * xc) ./ denom; % Still can't afford this last one
xc = xc ./ sqrt(denom); % now, C = xc' * xc

% Compute the SVD of xc:
[U, S, V] = svd(xc, 'econ'); % Could also use svds(xc, 2) to get only the two largest singular values

% We know that xc = U*S*V', therefore C = xc'*xc = V*S*U'*U*S*V'
% and because U is orthogonal, U'*U cancels out and we have
% C = V*S*S*V', and therefore we can treat V as the matrix of eigenvectors
% and the square of the singular values on the diagonal of S as the
% eigenvalues:
eigenvectors = V;
eigenvalues = S*S;

% Verify these are really the eigenvalues and eigenvectors of C = xc'*xc:
norm(xc'*(xc*eigenvectors) - eigenvectors*eigenvalues)

ans = 3.2419e-11

It seems like what you're computing may be similar to pca - perhaps that function could be useful.