Let $\rho$ be $n\times n$ symmetric matrix, thus the spectral decomposition of $\rho$ is $$ \rho=\sum _{i=1}^n e_i\left|\psi _i\right\rangle \left\langle \psi _i \right| $$ where $e_i$ and $\psi _i$ are the eigenvalues and the corresponding eigenvectors of $\rho$. I saw in some sources
$$ \sqrt{\rho }=\sum _{i=1}^n \sqrt{e_i}\left|\psi _i\right\rangle \left\langle \psi _i \right| $$ Is that correct? If so, how can I prove it?
Note that in a spectral decomposition, the $\psi_i$ form an orthonormal basis. From the fact that you said a "symmetric" matrix, I assume that $\rho$ and its eigenvectors have real entries.
With that said, it is correct that the operator $\sigma = \sum _{i=1}^n \sqrt{e_i}\left|\psi _i\right\rangle \left\langle \psi _i \right|$ satisfies $\sigma^2 = \rho$. We compute \begin{align} \left(\sum _{i=1}^n \sqrt{e_i}\left|\psi _i\right\rangle \left\langle \psi _i \right|\right)^2 &= \left(\sum _{i=1}^n \sqrt{e_i}\left|\psi _i\right\rangle \left\langle \psi _i \right|\right) \left(\sum _{j=1}^n \sqrt{e_j}\left|\psi _j\right\rangle \left\langle \psi _j \right|\right) \\ & = \sum_{i,j = 1}^n (\sqrt{e_i}\left|\psi _i\right\rangle \left\langle \psi _i \right|)\cdot (\sqrt{e_i}\left|\psi _j\right\rangle \left\langle \psi _j \right|) \\ & = \sum_{i,j = 1}^n \sqrt{e_i}\sqrt{e_j}\left|\psi _i\right\rangle \left\langle \psi _i | \psi _j\right\rangle \left\langle \psi _j \right| \\ & = \sum_{i = 1}^n \sqrt{e_i}^2\left|\psi _i\right\rangle \left\langle \psi _i \right| = \sum_{i = 1}^n e_i\left|\psi _i\right\rangle \left\langle \psi _i \right|. \end{align}
Note that if any of the $e_i$ are negative, then $\sigma$ will be a non-Hermitian complex symmetric.
If you meant to say "Hermitian" instead of symmetric, then the square root computed is neither symmetric nor Hermitian in the case of negative eigenvalues. It will, however, be a normal operator.
The expression $\sqrt{\rho}$ is only typically applicable to positive semidefinite operators $\rho$, which is to say that none of the eigenvalues of $\rho$ are negative. In this case, $\sqrt{\rho}$ is indeed defined as being the $\sigma$ described above, with all square roots of eigenvalues referring to the principal (non-negative) square root.