The notation $\Sigma$ is often used for covariance matrix Var($r_t$).
Using notation $\Psi_t = \text{Var}(\epsilon_t)$, a common modeling assumption is that the residuals are independent, $\epsilon_{i,t} \perp \epsilon_{j,t}, ~i\neq j$. Given this independence assumption, the covariance of residual returns matrix $\Psi$ is diagonal.
Some factor models (e.g. Barra) implement a concept of "linked specific risk" (LSR), perhaps to model two share classes of the same company. A residual return covariance matrix implementing a concept such as LSR would exhibit non-zero covariance for specific returns.
Idiosyncratic risk is asset specific. Does that imply that the covariance matrix ... is diagonal?
If your model assumes (your hypothesis asserts) idiosyncratic specific risk, the covariance of residual returns is diagonal. Should you collect sufficient data to reject this assumption, as the Barra model does, your residual returns would not be independent, and the residual return covariance matrix would not be diagonal.
would it make sense to impose this diagonality
Starting with a simple (maximum entropy) model "makes sense". Adding complexity under duress (rejecting the current hypothesis based upon data) as the Barra example illustrates makes sense. Modeling is an iterative process. Think deep fakes. If you can spot a feature that distinguishes the "real" from the "fake" image (or return series), impose a constraint that the synthesized return series' feature match the observed return series' feature.
the problem of inverting a large unrestricted estimated covariance matrix
A large estimation effort benefits from regularization. Two potential methods for your modeling efforts are:
- a Bayesian framework, where the prior is a simple model, and the observed data transforms this simple starting point; and,
- a maximum entropy framework, where features are added iteratively to a current model (as you "reject the null hypothesis" in response to statistically significant observed data features inconsistent with synthesized data features).
I wrestled these and related concepts in more detail in my thesis.
Edit: (response to first comment below)
Using Sherman-Morrison-Woodbury, and assuming you have a common factor
model that could be expressed as loadings on orthogonal common factor returns (rotate and scale original factor returns if necessary) plus residual returns, the effort to invert the covariance matrix is modest.
With a $k$-factor model, the only inversion required is of a $k \times k$ matrix $I_k + L^\textrm{T}\Psi^{-1}L$. The remainder of the process is matrix multiplication and matrix subtraction.
$$
\begin{align}
\Sigma &= L L^\textrm{T} + \Psi \\
\Sigma^{-1} &= \Psi^{-1} -
\Psi^{-1}L (I_k + L^\textrm{T}\Psi^{-1}L)^{-1}L^\textrm{T} \Psi^{-1} \\
\end{align}
$$
