Why do large LMs use the transpose of the word embeddings matrix in the classification head? All literature, guides and tutorials describing the construction of language models have used two separate matrices for the input and output projections:
To project one-hot token IDs into hidden states, a matrix $U$ of size $\mathrm{VocabSize} \times \mathrm{HiddenDim}$ is used.
To project the hidden states back into vocab space prior to softmax, a matrix $V$ of size $\mathrm{HiddenDim} \times \mathrm{VocabSize}$ is used.
However, the source code for transformers in the Huggingface library shows a completely different process that uses only one matrix:
The classification head takes hidden states $h$ at the top of the Transformer stack and then performs $\operatorname{softmax}( \operatorname{MLP}(h) \cdot U^T+b )$ to recover the vocab probability distribution.

Why exactly does this work? Can the MLP and transpose of $U$ be considered some form of pseudo inverse process? And why have I not seen this in literature for modern Transformer architectures? Is this something that dates back to the first RNN-LMs and as such need not be described? And if so, why does every tutorial I've seen (both for encoder-decoder RNNs and Transformers) not describe this and instead opt for using two separate matrices? Could it be that this is harder to optimize due to weight sharing?
I would appreciate some insight on this and references to literature since I am unable to find anything myself.

Turning my comment into an answer.
This practice is called parameter tying. They’ve set $U$ and $V$ to be equal to each other as a constraint, which regularizes the model. The information about token relatedness is assumed to carry over, at least to some degree, between the source and target—especially for related languages.
You ask about whether this dates back to "the first RNN-LMs". Not quite. The paper that popularized the technique in the context of neural sequence models is Press and Wolf (2017).
You also ask about the transpose—that's so that the matrices are conformant. The original matrix $U$ maps from vocabulary to a latent space; the $V := U^{\intercal}$ matrix maps back from a latent space to the vocabulary.
In terms of implementation and how widespread it is, it’s offered as a flag in the popular seq2seq libraries.
