How is explained variance in sparse PCA calculated?

Question

Sparse PCA is a technique proposed by Zou et all in this paper. In usual PCA the obtained loadings are orthonormal, and the resulting scores are uncorrelated. However, in sparse PCA you give up these orthogonality and uncorrelatedness in exange of sparse loadings that ease interpretation in high dimension.

Since the loadings are no longer orthonormal in SPCA, it is not as straightforward as in PCA to compute the explained variance, and in order to solve this problem, in the paper mentioned above (section 3.4) they propose an alternative way of calculating the variance that takes into account the correlations.

Supose your data matrix is $X$ and you compute $P$ the matrix of sparse loadings. Then you can obtain $Z=XP$ the matrix of scores. Now, as stated in section (3.4), you can decompose $Z=QR$ where Q is orthonormal and $R$ is upper triangular, and the explained variance from the sparse component $j$ in $P$ can be calculated as $R_{jj}^2$. Therefore, the total variance explained using sparse PCA is equal to $\sum_{j=1}^k R_{jj}^2$.

Now, I used the elasticnet package in R that implements the sparse PCA proposed in the paper and I tried to replicate those calculations for the variance:

# Load data
data("BostonHousing")
x = BostonHousing[, -c(2, 4, 9, 10, 14)] # Select non-categorical variables
total_variance_in_x = sum(diag(var(x)))
n = dim(x)[1]

# Obtain sparse components
n_comp=9
spca_fit1 = spca(x, K=n_comp, type="predictor", sparse="penalty", para=rep(10, 9))
total_var_spca = spca_fit1$var.all/(n-1)

# Check that the variance in x is the same as in SPCA.
abs(total_variance_in_x - total_var_spca) < 1e-3 
[1] TRUE

# Obtain matrices P and Z
p_spca = spca_fit1$loadings
z_spca = as.matrix(x)%*%p_spca 

# QR decomposition of modified PCs
z_spca_qr = qr(z_spca)
q_spca = qr.Q(z_spca_qr)
r_spca = qr.R(z_spca_qr)

# compute adjusted variance
variance_z_spca_qr_aproach = diag(r_spca)^2/(n-1)
total_var_spca_qr_aproach = sum(variance_z_spca_qr_aproach)

# Variance calculated using QR aproach do not match with the variance calculated above
total_var_spca_qr_aproach
[1] 130942.1
total_var_spca
[1] 9308.785

# The ratio of explained variance using QR differes with the one offered by the package
cumsum(spca_fit1$pev)
[1] 0.9051107 0.9872324 0.9940043 0.9972498 0.9993186 0.9997213 0.9998830 0.9999122 0.9999122
spca_variance_ratio_qr_aproach = cumsum(variance_z_spca_qr_aproach / total_var_spca_qr_aproach)
spca_variance_ratio_qr_aproach
[1] 0.9874552 0.9989841 0.9994996 0.9997414 0.9998887 0.9999797 0.9999941 1.0000000 1.0000000

So my questions are:

1. Why is total_variance_in_x equal to total_var_spca ? Is there no loss on variability due to the correlatedness and non-orthogonality?
2. Why is total_var_spca_qr_aproach much larger than total_var_spca? I am following the steps described in the paper in which the package is based so these values should match.
3. Am I doing something wrong? Did I misunderstand anything?