Staircase spectrum: is there a known solution for this problem? •   Let $A$ be a (real) dense (non-sparse) matrix with size $rp\times M$  ($M < rp$) with orthonormal columns (i.e. $A^TA = I_M$).
•   Similarly, let $B$ be a dense matrix with size $rq\times N$  ($N < rq$) with orthonormal columns.
•   Let $V=\bigoplus_{i=1}^qA.$ Equivalently, $V=I_q\bigotimes A,$ where $\bigoplus$ denotes the direct sum, and $\bigotimes$ denotes the Kronecker
product.
•   Let $W=\bigoplus_{j=1}^pB.$ Equivalently, $W=I_p\bigotimes B.$
•   Let $R=V^TW,$ and therefore, 
$$
R=\big(\bigoplus_{i=1}^qA^T\big)\big(\bigoplus_{j=1}^pB\big)
$$
$$
=\big(I_q\bigotimes A^T\big)\big(I_p\bigotimes B\big) 
$$
•   Assume:  $rpq > qM + pN$ (so that the augmented matrix $\left[\begin{array}{c|c}V & W\end{array}\right]$ is a tall matrix). 
•   Note:  $p\neq q$, $p>1$ and $q>1$. 

Prove that the singular values of $R$ come in $n$-tuples, with
$$
n=\text{GCD}(p,q),
$$ where $\text{GCD}(p,q)$ is the greatest common divisor of $p$ and $q$. In other words, prove that a singular value of $R$ has
a multiplicity of $\text{GCD}(p,q)$.

Note: I am ONLY interested in proving that the largest singular value of $R$ is distinct when $p$ and $q$ are relatively
prime ($\text{GCD}(p,q) = 1$). However, I think giving the larger picture of the behavior of the singular values of $R$ is far more interesting, and probably even helpful in proving the specific case I'm interested in. Please refer to the Matlab code below to verify that $n=\text{GCD}(p,q)$.

Matlab code:
clear all; close all;

r = 3; % integer >= 1

p = 12; q = 16; %17 % for the singular values of R to be distinct, p and q must be relatively prime integers and larger than 1 (i.e. GCD(p,q) = 1 and p > 1 and q > 1) 

M = 8; % # of columns of A
N = 14; % # of columns of B. M and N MUST satisfy 1 < M < rp, 1 < N < rq AND they must satisfy assumption that r*p*q/(q*M+p*N) > 1.

[A,~] = qr(rand(r*p,M),0); % size = rp X M (M < rp, A'A = I) % Note: you can orthonormalize (e.g. via Gram Schmidt or SVD) any arbitrary 
                       % (full rank) data matrix other than noise, I'm only using
                       % noise as a data matrix because it's conevniently easy to generate.
[B,~] = qr(rand(r*q,N),0); % size = rq X N (N < rq, B'B = I)

V = kron(eye(q),A); % size = q(rp) X qM
W = kron(eye(p),B); % size = p(rq) X pN

R = V'*W; % i.e. R = direct sum (over q) of (A') X direct sum (over p) of (B)

% assumption
r*p*q/(q*M+p*N) % check is greater than 1 % This assumption is required for the augmented matrix [V W] to be tall. 
gcd(p,q) % check for relative primeness of p and q. % This assumption is required for the singular values of R to be distinct.

figure;stem(svd(R));   
title 'note how the singular values of R come in GCD(p,q)-tuples'

 A: Let me write 
\begin{align}
V &= I_q \otimes A = \text{blkdiag}(A, ..., A), \\
W &= I_p \otimes B = \text{blkdiag}(B, ..., B).
\end{align}
We can write
$$
V^TW = I_q \otimes (A^TB)=\text{blkdiag}(A^TB,...,A^TB)
$$
only if $p=q$. But if $p \neq q$, we can still do the following:
$$
R=V^TW = I_n \otimes (V_{\text{sub}}^TW_{\text{sub}})=\text{blkdiag}(V_{\text{sub}}^TW_{\text{sub}},...,V_{\text{sub}}^TW_{\text{sub}}),
$$
with
\begin{align}
n &= \text{GCD}(p, q), \\
V_{\text{sub}} &= I_{q/n} \otimes A = \text{blkdiag}(A,...,A), \\
W_{\text{sub}} &= I_{p/n} \otimes B = \text{blkdiag}(B,...,B).
\end{align}
So in the end, it turns out that $R$ is a block diagonal matrix with $n=\text{GCD}(p,q)$ identical matrices. To compute the singular values of a block diagonal matrix, it suffices to compute the singular values of the matrices. Hence, the singular values of $R$ are simply $n$ times the singular values of $V_{\text{sub}}^TW_{\text{sub}}$.
Note that the fact that $A$ and $B$ are orthogonal is not used for this derivation. It also holds in case $A$ and $B$ are non-orthogonal.
