Calculate the variance explained in matrix Y by matrix X

I have two matrices corresponding to the same set of $n$ samples, with $j$ and $k$ variables, respectively ($j > 10000$, $k > 10000$). $X$ is an $n \times j$ matrix and $Y$ is an $n \times k$ matrix.

I'm interested to know what proportion of the overall variance of $Y$ can be explained by $X$. The ultimate goal is to compare $X$ to a similar matrix $X'$ and see which one explais $Y$ and what overlap in the explained variance these both matrices have.

My first approach was to try PLS, however I don't understand this method enough to know what I'm doing, and in any case, the implementation I tested (pls package in R) exited with an error indicating the failure to allocate half a terabyte of memory.

Furthermore, I was considering a following approach (possibly this is re-inventing the wheel; also, it might be completely wrong). PCA can partition the variance in both $X$ and $Y$, such that the resulting components have variances that sum up to 1, and are orthogonal. Therefore, I could try to calculate for the first component of $PCA(Y)$, what percentage of variance is explained using $X$; in R code, I'd try

set.seed(2708)
data(iris)
Y    <- iris[,1:4]
X    <- apply(Y, 2, function(x) x + rnorm( length(x), sd= 0.5))
pcaX <- prcomp(X, scale.= TRUE)
pcaY <- prcomp(Y, scale.= TRUE)

lm.m <- lm(pcaY$x[,1] ~ pcaX$x)
print( summary( lm.m )$r.squared) #  0.9319846 Then I could add these values for each component, multiplied by the fraction of variance in Y that this component holds: Yvars <- pcaY$sdev^2/sum(pcaY$sdev^2) foo <- apply(pcaY$x, 2, function(y) summary(lm(y ~ pcaX$x))$r.squared)
print(sum(foo * Yvars))
#  0.7993989

X2    <- apply(Y, 2, function(x) x + rnorm(length(x), sd= 5))
pcaX2 <- prcomp(X2, scale.= T)