Why is computation of scores in Sparse PCA different from T=XP? I am recently learning Sparse PCA. From a lately published paper All sparse PCA models are wrong, but some are useful. Part I: Computation of scores, residuals and explained variance I learned that

since the loadings of sPCA are not orthogonal, we should use the correction formula below to compute the scores, where the superscript '+' indicates the Moore-Penrose inverse:

$$
\mathbf{T}=\mathbf{X}\mathbf{P}\left(\mathbf{P}^\text{T}\mathbf{P}\right)^+
$$
Actually I cannot figure out the reason behind it. As far as I know previously, $\mathbf{P}$ contains coefficient vectors which can construct Principal Component directly, so $\mathbf{T}=\mathbf{X}\mathbf{P}$ already stands for the scores.
I've tried to turn to Google, but didn't find any related topics. Could you please tell me why we should add $\left(\mathbf{P}^\text{T}\mathbf{P}\right)^+$ to the right side of the formula just because the loadings are not orthogonal?
Thanks in advance.
 A: $$P^TP$$
Is the inner product of the coefficient vectors (PCs, loadings ...) and returns their covariance matrix. This calculates the non-orthogonality of the coefficient vectors .
$$(P^TP)^+$$
Is the inverse of the covariance of the coefficient vectors, so if you post multiply by this then you correct for any covariance between the coefficient vectors.
This is critical because any sequential iterative algorithm for calculating PCA leads to an accumulation of errors between each PC. These errors compromise the orthogonality of the PCs. If you want to do a comparison between two algorithms that should give the same result compare SVD (linear algebra based so gives true PCA) and NIPALs (iterative estimation of sequential PCs) and you will see the errors building up and non-orthogonality. For sparse PCA, which actively shrinks variable contribution there is even more scope for non-orthogonality if shrinkage is tightly linked across PCs.
Here's some graphs from supporting information for a paper I have submitted that illustrate this. The data is a simulated dataset that has had no noise added. The first shows the correlation between equally ranked NIPALS and SVD against PC rank in the orange dots. You can see that initially the correlation between the two algorithms is high, but from 19 on it drops dramatically. However, vector identity is not an arbitrary rank label, but its shape and orientation, so we should match by closest correlation. We see that even with this correction significant divergence emerges.

So next I show how the rank of PC (identity based on correlation to SVD PCs, not rank) changes between the two algorithms. At first the ranks increase in tandem but become more disrupted. This is due to the accumulation of computation errors. Of particular note is the blip at rank 19 in the NIPALs model. This matches rank 1 in SVD, despite rank 1 in NIPALS already matching that!

Here we now compare PC 1 and 19 in the NIPALS algorithms - you can see they are indeed very similar apart from PC19 having significant contribution from computational noise.

finally now we compare the PC1 between the two algorithms - the difference is tiny (scale is $10^{-14}$ compared with $10^{-1}$ for the actual PCs), but as we accumulate successive small errors they grow unless renormalised.

