Error message when trying to compute EOF with covariance matrix.

2 views (last 30 days)
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

Accepted Answer

Christine Tobler
Christine Tobler on 6 May 2022
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.

More Answers (0)

Categories

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by