Let $Y \in \mathbb R^{n \times c}$ and $r$ be an integer with $1 \le r \le \operatorname{rank}(Y)$. Consider the problem

$$\text{minimize} \|Y-X\|_{\text{Fro}}^2\text{ over }X \in \mathbb R^{n \times c} \text{ subject to rank}(X) \le r,$$ which seeks a $n$-by-$c$ matrix of rank $r$ which is closest to $Y$ w.r.t the Frobenius norm. The Eckart–Young–Mirsky theorem shows that the choice $X = \hat{Y}(r) := Y\Pi(r)$ solves this problem, where $\Pi(r) := V(r)V(r)^T$ the projection onto the subspace spanned by the first $r$ principal vectors $V(r) : = [V_1,\ldots,V_r] \in \mathbb R^{c \times r}$ of $Y$.

**Question:** What is $\hat{df}(r) := \operatorname{tr}\left(\frac{\partial vec(\hat{Y}(r))}{\partial vec(Y)}\right) = \sum_{i,j,k,l} \frac{\partial \hat{Y}_{ij}(r)}{\partial Y_{kl}}$ ?

**Motivation:** Such computations are necessary in applying SURE (Stein's Unbiased Risk Estimator) theory for tuning penalized least squares models, as a powerful alternative to cross-validation. Here "df" stands for "degrees of freedom" and is a measure of the complexity of the approximator model.

**Observation:** It's not hard to show that
$$\hat{df}(r) \ge \hat{df}_{\text{naive}}(r),$$

where $\hat{df}_{\text{naive}}(r):= nc - (n-r)(c-r)$ is the least number of free parameters needed to completely specify an $n$-by-$c$ matrix of rank $r$. Theorem 5.2 of this paper computes the correction $\hat{df}(r) - \hat{df}_{\text{naive}}(r)$ exactly to be

$$\hat{df}(r) - \hat{df}_{\text{naive}}(r) = 2\sum_{k=1}^r \sum_{l=r+1}^{\operatorname{rank}(Y)}\frac{\lambda_l}{\lambda_k-\lambda_l} \ge 0, $$ where $\lambda_1 \ge \lambda_2 > \cdots > \lambda_r > \lambda_{r+1} \ge \cdots \ge \lambda_{\operatorname{rank}(Y)}$ are the eigenvalues of $Y^TY$. However, the method used by the referred paper is direct "brute-force" and involves some very long formulae. I was wondering whether there could be a cleaner (i.e synthetic) way to obtain this.