The PSD square root is continuous

Idea. A compactness-plus-uniqueness argument: bounded sequences have convergent subsequences, and the limit must satisfy the defining property of $\sqrt S$, which determines it uniquely.

Step 0: the definition. For $S$ PSD, the spectral theorem gives $S=Q\operatorname{diag}(\mu_i)Q^T$ with $\mu_i\ge0$, and $\sqrt S=Q\operatorname{diag}(\sqrt{\mu_i})Q^T$ is PSD with $(\sqrt S)^2=S$. *Uniqueness* of the PSD square root is a standard fact (used here without proof): if $T\ge0$ and $T^2=S$ then $T=\sqrt S$.

Step 1: boundedness. $T_k=\sqrt{S_k}$ is symmetric, so $\|T_k\|^2=\|T_k^2\|=\|S_k\|$ (for a symmetric matrix, $\|T\|=\rho(T)$ and $\rho(T^2)=\rho(T)^2$). Since $S_k\to S$ is convergent, $\|S_k\|$ is bounded, hence so is $(T_k)$.

Step 2: limit points. By Bolzano-Weierstrass, any subsequence of $(T_k)$ has a further subsequence $T_{k_j}\to T$. The cone of PSD matrices is closed, so $T$ is PSD. By continuity of multiplication (Problem 67), $T^2=\lim T_{k_j}^2=\lim S_{k_j}=S$. By uniqueness, $T=\sqrt S$.

Step 3: convergence. Every subsequence of $(T_k)$ has a sub-subsequence converging to the same limit $\sqrt S$; hence $T_k\to\sqrt S$ (if not, some subsequence would stay at distance $\ge\epsilon$ from $\sqrt S$, contradicting Step 2).

Remark.
- *Quantitative version.* The map is Hölder: $\|\sqrt A-\sqrt B\|\le\sqrt{\|A-B\|}$ for PSD $A,B$ (known; compare $n=1$, $\sqrt x$ at $0$). It is not Lipschitz at singular matrices: $\sqrt{\varepsilon}$ vs $\varepsilon$.
- On $\{A\ge cI\}$ with $c>0$ one has, in the Frobenius norm, $\|\sqrt A-\sqrt B\|_F\le\frac{\|A-B\|_F}{2\sqrt c}$. (With $X=\sqrt A-\sqrt B$ one has $\sqrt A\,X+X\sqrt B=A-B$, and in the eigenbases the coefficient of $x_{ij}$ is $\alpha_i+\beta_j\ge2\sqrt c$.) A numeric test confirmed the bound ($0.388\le1.551$).
- Used in Problem 90 to show the polar decomposition is continuous.