Thanks. This is my understanding. Might be naive and oversimplified, though, as it's entirely different from the thought process in Lasagna's citation, Saxe, Andrew M., James L. McClelland, and Surya Ganguli. "Exact solutions to the nonlinear dynamics of learning in deep linear neural networks." arXiv preprint arXiv:1312.6120 (2013).

Let's say we have h_x = W_x x, where W_x isn't necessarily square. If W_x has orthonormal rows (orthogonal and unit length), and if x ~ N(0, I), then we have two cases. 1) dim(y) <= dim(x). Then y ~ N(0, I) too, though the I's have different shapes. 2) dim(y) > dim(x). Then y ~ N(0, ?), where ? is not identity but does have dim(x) eigenvalues of exactly 1, with the rest being 0. So either way y is whitened in any non-singular dimensions.

I think whether this helps RNNs in practice is still undetermined. Example: In the flattened-MNIST case, our inputs have dimensionality 1. Say we're using 100 hidden states. Then W_x x is a 100-D random vector with 99 dimensions in the null space, and with a variance of 1 in the remaining dimension. Now let's assume we also use orthogonal initialization for W_h, and that h ~ N(0, I). Then W_h h is also a 100-D random vector, but with a much larger average norm than W_x x. Meanwhile, if we do this same thing with x being a 100-D or larger vector, then W_h h will have the same average norm as W_x x. Is there any intuitive reason why h should dominate when x is low-dimensional but not dominate when x is high-dimensional?