[proofplan]
We conjugate the generally non-symmetric matrix $W^{-1}B$ to the symmetric positive semidefinite matrix $A := W^{-1/2}BW^{-1/2}$. This reduces the determinant and trace computations to diagonal entries of $A$. The same conjugation identifies $(B+W)^{-1}B$ with $(I_p+A)^{-1}A$, whose eigenvalues are $\theta_k/(1+\theta_k)$. Finally, the Roy statistic is already defined as the largest nonzero root.
[/proofplan]
[step:Conjugate the MANOVA root matrix to a symmetric positive semidefinite matrix]
Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $W \in \mathbb{R}^{p \times p}$ is symmetric positive definite, the matrix square-root theorem gives a unique symmetric positive definite matrix $W^{1/2} \in \mathbb{R}^{p \times p}$ such that $(W^{1/2})^2=W$; let $W^{-1/2}$ denote its inverse. Define
\begin{align*}
A := W^{-1/2}BW^{-1/2} \in \mathbb{R}^{p \times p}.
\end{align*}
The matrix $A$ is symmetric because
\begin{align*}
A^\top = (W^{-1/2}BW^{-1/2})^\top = W^{-1/2}B^\top W^{-1/2} = W^{-1/2}BW^{-1/2}=A,
\end{align*}
and it is positive semidefinite because, for every $x \in \mathbb{R}^p$,
\begin{align*}
x^\top A x
= x^\top W^{-1/2}BW^{-1/2}x
= (W^{-1/2}x)^\top B(W^{-1/2}x)
\geq 0.
\end{align*}
Moreover,
\begin{align*}
W^{1/2}(W^{-1}B)W^{-1/2}
= W^{-1/2}BW^{-1/2}
= A.
\end{align*}
Thus $W^{-1}B$ is similar to $A$, so the nonzero eigenvalues of $A$ are exactly $\theta_1,\dots,\theta_s$, counted with multiplicity.
Since $A$ is real symmetric positive semidefinite, the finite-dimensional spectral theorem for real symmetric matrices gives an orthogonal matrix $Q \in \mathbb{R}^{p \times p}$ and numbers $\lambda_1,\dots,\lambda_p \in [0,\infty)$ such that
\begin{align*}
A = QDQ^\top,
\qquad
D := \operatorname{diag}(\lambda_1,\dots,\lambda_p).
\end{align*}
After reordering the diagonal entries, we may take
\begin{align*}
\lambda_k = \theta_k \quad \text{for } 1 \leq k \leq s,
\qquad
\lambda_k = 0 \quad \text{for } s < k \leq p.
\end{align*}
(citing a result not yet in the wiki: finite-dimensional spectral theorem for real symmetric matrices)
[guided]
The matrix $W^{-1}B$ need not be symmetric, so we first replace it by a similar symmetric matrix. Let $I_p \in \mathbb{R}^{p \times p}$ denote the identity matrix. Since $W \in \mathbb{R}^{p \times p}$ is symmetric positive definite, the matrix square-root theorem applies and gives a unique symmetric positive definite matrix $W^{1/2} \in \mathbb{R}^{p \times p}$ satisfying $(W^{1/2})^2=W$; let $W^{-1/2}$ denote its inverse. Define
\begin{align*}
A := W^{-1/2}BW^{-1/2} \in \mathbb{R}^{p \times p}.
\end{align*}
This is the correct conjugate because it preserves eigenvalues while making symmetry visible. Indeed,
\begin{align*}
A^\top
= (W^{-1/2}BW^{-1/2})^\top
= W^{-1/2}B^\top W^{-1/2}
= W^{-1/2}BW^{-1/2}
= A,
\end{align*}
where we used the symmetry of both $B$ and $W^{-1/2}$. It is also positive semidefinite: for each $x \in \mathbb{R}^p$,
\begin{align*}
x^\top A x
= (W^{-1/2}x)^\top B(W^{-1/2}x)
\geq 0,
\end{align*}
because $B$ is positive semidefinite.
Now compare $A$ with $W^{-1}B$. Multiplying $W^{-1}B$ on the left by $W^{1/2}$ and on the right by $W^{-1/2}$ gives
\begin{align*}
W^{1/2}(W^{-1}B)W^{-1/2}
= W^{-1/2}BW^{-1/2}
= A.
\end{align*}
Therefore $A$ and $W^{-1}B$ are similar matrices, and similar matrices have the same eigenvalues with algebraic multiplicity. Hence the nonzero eigenvalues of $A$ are precisely the MANOVA roots $\theta_1,\dots,\theta_s$.
Because $A$ is real symmetric positive semidefinite, the finite-dimensional spectral theorem for real symmetric matrices applies: there is an orthogonal matrix $Q \in \mathbb{R}^{p \times p}$ and a diagonal matrix
\begin{align*}
D := \operatorname{diag}(\lambda_1,\dots,\lambda_p)
\end{align*}
with $\lambda_k \geq 0$ for every $k$ such that
\begin{align*}
A = QDQ^\top.
\end{align*}
Since the nonzero eigenvalues are exactly $\theta_1,\dots,\theta_s$, we reorder the diagonal entries so that
\begin{align*}
\lambda_k = \theta_k \quad \text{for } 1 \leq k \leq s,
\qquad
\lambda_k = 0 \quad \text{for } s < k \leq p.
\end{align*}
(citing a result not yet in the wiki: finite-dimensional spectral theorem for real symmetric matrices)
[/guided]
[/step]
[step:Compute Wilks' Lambda from the determinant factorization]
Using $B+W = W^{1/2}(I_p+A)W^{1/2}$, we obtain
\begin{align*}
\det(B+W)
&= \det(W^{1/2})\det(I_p+A)\det(W^{1/2})\\
&= \det(W)\det(I_p+A).
\end{align*}
Therefore
\begin{align*}
\Lambda
= \frac{\det W}{\det(B+W)}
= \frac{1}{\det(I_p+A)}.
\end{align*}
Since $A=QDQ^\top$ and $Q$ is orthogonal,
\begin{align*}
\det(I_p+A)
&= \det\bigl(Q(I_p+D)Q^\top\bigr)\\
&= \det(I_p+D)\\
&= \prod_{k=1}^p (1+\lambda_k)\\
&= \prod_{k=1}^s (1+\theta_k).
\end{align*}
Thus
\begin{align*}
\Lambda
= \prod_{k=1}^s \frac{1}{1+\theta_k}.
\end{align*}
For each $1 \leq k \leq s$, the definition $\hat{\rho}_k^2=\theta_k/(1+\theta_k)$ gives
\begin{align*}
1-\hat{\rho}_k^2
= 1-\frac{\theta_k}{1+\theta_k}
= \frac{1}{1+\theta_k}.
\end{align*}
Hence
\begin{align*}
\Lambda
= \prod_{k=1}^s(1-\hat{\rho}_k^2).
\end{align*}
[/step]
[step:Compute the Pillai trace from the transformed eigenvalues]
First rewrite $(B+W)^{-1}B$ using $A$. Since
\begin{align*}
B+W = W^{1/2}(I_p+A)W^{1/2},
\end{align*}
we have
\begin{align*}
(B+W)^{-1}B
&= W^{-1/2}(I_p+A)^{-1}W^{-1/2}B.
\end{align*}
Multiplying on the left by $W^{1/2}$ and on the right by $W^{-1/2}$ gives
\begin{align*}
W^{1/2}\bigl((B+W)^{-1}B\bigr)W^{-1/2}
&= (I_p+A)^{-1}W^{-1/2}BW^{-1/2}\\
&= (I_p+A)^{-1}A.
\end{align*}
Thus $(B+W)^{-1}B$ is similar to $(I_p+A)^{-1}A$. Since $A=QDQ^\top$,
\begin{align*}
(I_p+A)^{-1}A
&= Q(I_p+D)^{-1}DQ^\top.
\end{align*}
The diagonal entries of $(I_p+D)^{-1}D$ are $\lambda_k/(1+\lambda_k)$, so the trace is
\begin{align*}
V
&= \operatorname{tr}\bigl((B+W)^{-1}B\bigr)\\
&= \operatorname{tr}\bigl((I_p+A)^{-1}A\bigr)\\
&= \sum_{k=1}^p \frac{\lambda_k}{1+\lambda_k}\\
&= \sum_{k=1}^s \frac{\theta_k}{1+\theta_k}.
\end{align*}
By the definition of $\hat{\rho}_k^2$,
\begin{align*}
V = \sum_{k=1}^s \hat{\rho}_k^2.
\end{align*}
[guided]
The Pillai statistic involves $(B+W)^{-1}B$, so we need to express this matrix in the same diagonal coordinates as $A$. The identity
\begin{align*}
B+W = W^{1/2}(I_p+A)W^{1/2}
\end{align*}
implies
\begin{align*}
(B+W)^{-1}
= W^{-1/2}(I_p+A)^{-1}W^{-1/2}.
\end{align*}
Therefore
\begin{align*}
(B+W)^{-1}B
= W^{-1/2}(I_p+A)^{-1}W^{-1/2}B.
\end{align*}
To identify the eigenvalues, conjugate by $W^{1/2}$:
\begin{align*}
W^{1/2}\bigl((B+W)^{-1}B\bigr)W^{-1/2}
&= (I_p+A)^{-1}W^{-1/2}BW^{-1/2}\\
&= (I_p+A)^{-1}A.
\end{align*}
Thus $(B+W)^{-1}B$ and $(I_p+A)^{-1}A$ are similar and have the same trace.
Now use the diagonalization $A=QDQ^\top$. Since $I_p+A=Q(I_p+D)Q^\top$, its inverse is $Q(I_p+D)^{-1}Q^\top$, and hence
\begin{align*}
(I_p+A)^{-1}A
= Q(I_p+D)^{-1}DQ^\top.
\end{align*}
This matrix has eigenvalues
\begin{align*}
\frac{\lambda_1}{1+\lambda_1},\dots,\frac{\lambda_p}{1+\lambda_p}.
\end{align*}
Taking the trace, which is invariant under similarity and equals the sum of diagonal entries for a diagonal matrix, gives
\begin{align*}
V
&= \operatorname{tr}\bigl((B+W)^{-1}B\bigr)\\
&= \operatorname{tr}\bigl((I_p+A)^{-1}A\bigr)\\
&= \sum_{k=1}^p \frac{\lambda_k}{1+\lambda_k}\\
&= \sum_{k=1}^s \frac{\theta_k}{1+\theta_k}.
\end{align*}
The final equality uses $\lambda_k=\theta_k$ for $1 \leq k \leq s$ and $\lambda_k=0$ for $s<k\leq p$. Since $\hat{\rho}_k^2=\theta_k/(1+\theta_k)$, this is exactly
\begin{align*}
V = \sum_{k=1}^s \hat{\rho}_k^2.
\end{align*}
[/guided]
[/step]
[step:Compute the Hotelling Lawley trace from the MANOVA roots]
Since $W^{-1}B$ is similar to $A$, the trace is invariant under this similarity:
\begin{align*}
U
= \operatorname{tr}(W^{-1}B)
= \operatorname{tr}(A).
\end{align*}
Using $A=QDQ^\top$,
\begin{align*}
\operatorname{tr}(A)
= \operatorname{tr}(D)
= \sum_{k=1}^p \lambda_k
= \sum_{k=1}^s \theta_k.
\end{align*}
Therefore
\begin{align*}
U = \sum_{k=1}^s \theta_k.
\end{align*}
[/step]
[step:Identify Roy's largest root statistic]
The roots were ordered as
\begin{align*}
\theta_1 \geq \theta_2 \geq \cdots \geq \theta_s > 0.
\end{align*}
By definition, Roy's largest-root statistic $\Theta$ is the largest nonzero eigenvalue of $W^{-1}B$. Hence
\begin{align*}
\Theta = \theta_1.
\end{align*}
Combining this identity with the determinant and trace computations above gives all four asserted root expressions.
[/step]
