Here is another proof adapted from my prof's lecture notes (maybe it's just a more noob-friendly explanation to Sheldon Axler's answer):
Step 1: Find a square root operator
Suppose that $A$ is the non-negative bounded linear operator on our Hilbert space, let
$$
p_n(x)\to\sqrt{x}
$$
be approximating polynomials due to Weierstrass Approximation Theorem, and let
$$
Y=\lim_{n\to\infty}p_n(A)=
\lim_{n\to\infty}\sum_{\lambda\in\sigma_p(A)}p_n(\lambda)P_\lambda ,
$$
where the last equal sign is due to Spectral Theorem:
- $\sigma_p(A)$ is the set of eigenvalues, the so-called point-spectrum of $A$;
- $P_\lambda$ is the orthogonal projection to the eigenspace of eigenvalue $\lambda$;
- We can write $p_n(A)$ like this, because $A^m=(\sum_{\lambda\in\sigma_p(A)}\lambda P_\lambda)^m=\sum_{\lambda\in\sigma_p(A)}\lambda^mP_\lambda$, since $P_\lambda P_\mu = \delta_{\lambda, \mu} P_\lambda$. Also see this Wikipedia page).
It can be proven that $Y$ fulfills $Y^2=A$.
Step 2: Show that it is unique
Suppose $X$ is another non-negative bounded linear operator such that $X^2=A$, then
$$
XA-AX=X^3-X^3=0 ,
$$
so the elements $A$ and $X$ commute with respect to function composition. This implies that
$$
Xp_n(A)-p_n(A)X=0 \qquad \forall n ,
$$
and
$$
XY-YX=\lim_{n\to\infty}\big(Xp_n(A)-p_n(A)X\big)=0 .
$$
Claim: $\mathbf{(X − Y)^3 = 0}$
Because
$$
(X − Y)X(X − Y) + (X − Y)Y(X − Y) \\
= (X^2-Y^2)(X-Y) = (A-A)(X-Y) = 0
$$
and both $(X − Y)X(X − Y)$ and $(X − Y)Y(X − Y)$ are non-negative operators, we have
$$(X − Y)X(X − Y) = (X − Y)Y(X − Y) = 0$$
and thus
$$
(X − Y)^3 = (X − Y)X(X − Y) − (X − Y)Y(X − Y) = 0 .
$$
Final Steps
We also have $(X − Y)^4=0$, and
$$
0 = \|(X − Y)^4\| = \|X − Y\|^4 .
$$
This is because $(X − Y)^*=X-Y$ self-adjoint, and $\|T^\ast T\|=\|T^\ast\|\|T\|$ for a bounded linear operator $T$ on a Hilbert space (the so-called $C^\ast$-identity).
So $X=Y$.