What is the correct way to calculate the explained variance of each EOF as calculated from a gappy data set?

Question

I am trying to determine the correct amount of variance explained by each mode of an Empirical Orthogonal Function (EOF) analysis (similar to "PCA") as applied to a gappy data set. (i.e., containing NaNs). The following question builds on an earlier one that I had regarding the differing results obtained from the decomposition of the data set's covariance matrix using either eigen or svd. In essence, the problem is that I have read that both decompositions can be used interchangeably for obtaining the EOFs from a square covariance matrix. This does seem to be the case when the data set is not gappy (as illustrated below):

###Make complete and gappy data set
set.seed(1)
x <- 1:100
y <- 1:100
grd <- expand.grid(x=x, y=y)

#complete data
z <- matrix(rnorm(dim(grd)[1]), length(x), length(y))
image(x,y,z, col=rainbow(100))

#gappy data
zg <- replace(z, sample(seq(z), length(z)*0.5), NaN)
image(x,y,zg, col=rainbow(100))


###Covariance matrix decomposition
#complete data
C <- cov(scale(z), use="pair")
E <- eigen(C)
S <- svd(C)

#sum of lambda
sum(E$values)
sum(S$d)
sum(diag(C))

The sum of lambda in both eigen and svd equals the sum of the diagonal of the covariance matrix. So far, so good - Both methods explain the correct amount of variance. The next example does the same routine for a gappy version of the data set (50% NaNs):

#gappy data (50%)
Cg <- cov(scale(zg), use="pair")
Eg <- eigen(Cg)
Sg <- svd(Cg)

#sum of lambda
sum(Eg$values)
sum(Sg$d)
sum(diag(Cg))

And here we see that the lambda values calculated by svd are greater than the sum of the diagonal of the covariance matrix. Those calculated by eigen are equal. However, because the covariance matrix is no longer positive definite, there are some negative trailing lambda values. In my previous question I showed that this tendency becomes greater with increasing gappiness.
So, I can live with this if need be, but now I'm concerned about how to correctly assign how much of the data set's variance is explained by each EOF. This should be lambda/sum(lambda). When I plot the cumulative explained variance of the EOFs, you will see the problem - because the eigen decomposition contains some negative eigenvalues, the slope of cumulative explained variance is steeper and bell-shaped:

#cumulative explained variance of the EOFs
E.cumexplvar <- cumsum(E$values/sum(E$values))
S.cumexplvar <- cumsum(S$d/sum(S$d))
Eg.cumexplvar <- cumsum(Eg$values/sum(Eg$values))
Sg.cumexplvar <- cumsum(Sg$d/sum(Sg$d))


###plot the cumulative explained variance
png("cumexplvar.png", width=8, height=4, units="in", res=200)
par(mfcol=c(1,2))
YLIM <- range(c(E.cumexplvar, S.cumexplvar, Eg.cumexplvar, Sg.cumexplvar))

plot(E.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="non-gappy")
points(S.cumexplvar, t="o", pch=2, col=2)
abline(h=1, col=8, lty=2)
legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)

plot(Eg.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="gappy")
points(Sg.cumexplvar, t="o", pch=2, col=2)
abline(h=1, col=8, lty=2)
legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)

dev.off()

The problem may be that I should be using the sum of the absolute eigenvalues to assign their explained variance, but this also leaves me to wonder how to interpret the explained variance of the negative eigenvalues. I would be very grateful for any insight, as this is not an issue that I have come across in any reference regarding EOF as applied to gappy data.