Why is this algorithm for PLS correct?

I am studying Partial Least Square Regression (PLS), and I am not able to understand how the algorithm for performing a PLS factorization works.

PLS is composed of two parts. First, we factor two matrices $$\mathbf{X}$$ and $$\mathbf{Y}$$, and then we perform the regression. I am stuck in the factorization part.

The hypotheses are that we have to matrices $$\mathbf{X} \in \mathbb{R}^{m\times p}$$ and $$\mathbf{Y} \in \mathbb{R}^{m \times q}$$, such that each column of the matrices has zero mean (i.e., they are mean-centered). Let us set $$k:= \max \left\lbrace p, q\right\rbrace$$.

If I understood correctly, we want to find four orthonormal matrices $$\mathbf{T} \in \mathbb{R}^{m \times k}$$, $$\mathbf{P} \in \mathbb{R}^{k\times p}$$, $$\mathbf{U} \in \mathbb{R}^{m \times k}$$ and $$\mathbf{Q} \in \mathbb{R}^{k \times q}$$ such that $$\begin{equation*} \mathbf{T} = \mathbf{X}\mathbf{P} \quad \text{and} \quad \mathbf{U} = \mathbf{Y}\mathbf{Q}. \end{equation*}$$

Furthermore, we want the sample variance-covariance matrices of $$\mathbf{Y}^{\top}\mathbf{X}\mathbf{P}$$ and of $$\mathbf{X}^{\top}\mathbf{Y}\mathbf{Q}$$ to be diagonal, and with eigenvalues listed in order of decreasing magnitude (the eigenvalues are all non-negative, as it follows from the definition of the variance-covariance matrix).

Since $$\mathbf{X}$$ and $$\mathbf{Y}$$ are mean-centered, the sample variance-covariance matrices are simply $$\begin{equation*} \text{var}(\mathbf{Y}^{\top}\mathbf{X}\mathbf{P}) = \frac{1}{k-1}\mathbf{P}^{\top}\mathbf{X}^{\top}\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X}\mathbf{P} \end{equation*}$$ and $$\begin{equation*} \text{var}(\mathbf{X}^{\top}\mathbf{Y}\mathbf{Q}) = \frac{1}{k-1}\mathbf{Q}^{\top}\mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}\mathbf{Y}\mathbf{Q}. \end{equation*}$$ So, by considering the spectral decompositions $$\begin{equation*} \mathbf{X}^{\top}\mathbf{Y}\mathbf{Y}^{\top}\mathbf{X} = \mathbf{U}\mathbf{D_1}\mathbf{U}^{\top} \end{equation*}$$ and $$\begin{equation*} \mathbf{Y}^{\top}\mathbf{X}\mathbf{X}^{\top}\mathbf{Y} = \mathbf{V}\mathbf{D_2}\mathbf{V}^{\top} \end{equation*}$$ (a spectral decomposition is a diagonalization where $$\mathbf{U}$$ and $$\mathbf{V}$$ are orthonormal matrices and where $$\mathbf{D_1}$$ and $$\mathbf{D_2}$$ are diagonal), where we can suppose that the eigenvalues are listed in order of decreasing magnitude, we have that the requirements of the PLS are satisfied by $$\begin{equation*} \mathbf{P} = \mathbf{U} \quad \text{and} \quad \mathbf{Q} = \mathbf{V}. \end{equation*}$$ Note that the non-zero eigenvalues of $$\mathbf{D_1}$$ and $$\mathbf{D_2}$$ form the same set, as observed here Relation between eigenvalues of $A^{\top}BB^{\top}A$ and $B^{\top}AA^{\top}B$.

So far so good, the problems arise when looking at the algorithm for computing the factorization. The algorithm, which is similar to the standard algorithm for computing eigenvalues and eigenvectors (the power method), proceeds in this way.

for i = 1:k
y = random vector
loop % iterate until t and u stop changing
p = X' u / norm(X' u) % X' means X transposed
t = X p
q = Y' t / norm(Y' t)
u = Y q
t(i) = t
u(i) = u
% deflation, as in the power method
X = X - t p'
Y = Y - u q'


My question: why is this algorithm correct?

Let us denote the $$i$$-th column of $$\mathbf{P}$$ with $$\mathbf{p_i}$$ (and let us use the same notation also for the other matrices). Does algorithm implicitly say that, for each $$i \in \min \left\lbrace 1, \ldots, k \right\rbrace$$ $$\begin{equation*} \mathbf{p_i} = \mathbf{X}^{\top}\mathbf{u}/\left\Vert \mathbf{X}^{\top} \mathbf{u}\right\Vert \quad \text{and} \quad \mathbf{q_i} = \mathbf{Y}^{\top}\mathbf{t}/\left\Vert \mathbf{Y}^{\top} \mathbf{t}\right\Vert? \end{equation*}$$

I have been trying proving these two identities for hours, but no way. So maybe the correctness of the algorithm does not imply those identities. Or maybe my formulation of PLS is not the one which the algorithm solves. In particular, maybe PLS does not require to have $$\text{var}(\mathbf{Y}^{\top}\mathbf{X}\mathbf{P}) = \mathbf{D_1}$$ and $$\text{var}(\mathbf{X}^{\top}\mathbf{Y}\mathbf{Q}) = \mathbf{D_2}$$, and that $$\mathbf{T}, \mathbf{P}, \mathbf{U}, \mathbf{Q}$$ are all orthonormal. In literature I have found the condition that for all $$i$$ we need to have \begin{equation*} \begin{aligned} &\mathbf{t_i}^{\top}\mathbf{t_i} = 1, \\ &\mathbf{u_i}^{\top}\mathbf{u_i} = 1, \\ &\mathbf{t_i}^{\top}\mathbf{u_i} \text{ maximum among the possible choice of \mathbf{t_i} and \mathbf{u_i}}. \end{aligned} \end{equation*} Is this requirement equivalent to the ones concerning the matrices $$\mathbf{D_1}$$ and $$\mathbf{D_2}$$ that I have listed above? And what can we say about $$\mathbf{t_i}^{\top}\mathbf{u_j}$$ for $$i \neq j$$? Are these value equal to zero?

$$\infty$$ thanks to those who help me!

Related: NIPALS for PCA and PLS, I have not found a good book either, I have been studying PLS putting together pieces from the web, and proceeding this way can be difficult. Also https://stats.stackexchange.com/questions/269306/understanding-nipals-algorithm-for-pls