Define $\mathbf{f}(A) = A^2$, for $A \in \mathbb{R}^{n \times n}$.
(a) Applying the Inverse Function Theorem, show that every matrix $B$ in a neighbourhood of $I$ has (atleast) 2 square roots $A$, that is $A^2 = B$, each varying as a $C^1$ function of $B$.
(b) Can you decide if there are precisely 2 or more ? (Hint: in 2x2 case, what is $D\mathbf{f}\left ( \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \right )$ ?
Now I can apply the inverse function theorem to $\mathbf{f}(A)$ at $A =I$ to conclude the inverse exists. But I cannot see how to apply it and ``show that every matrix $B$ in a neighbourhood of $I$ has (atleast) 2 square roots $A$" is to be solved and the bit about number of square roots.
I would appreciate some hints. This is a problem from Ted Shifrin's book on Mutivariable Mathematics.
(a) Let $A=diag(I_p,-I_q)$ where $p+q=n$; then $f(A)=I_n$ and
$Df_A:H=\begin{pmatrix}P_p&Q\\R&S_q\end{pmatrix}\in M_n\rightarrow AH+HA=\begin{pmatrix}2P&0\\0&-2S\end{pmatrix}\in M_n$.
Then $rank(Df_A)=p^2+q^2$ and $Df_A$ is an isomorphism iff $p=0$ or $q=0$, that is $A=\pm I_n$. Finally, for every neighborhood $V$ of $I_n$, there is $U$ a neighborhood of $\pm I_n$ s.t. $f(U)\subset V$ and $f_{|U}:U\rightarrow f(U)$ is a $C^{\infty}$diffeomorphism.
Remark. Every solution of $X^2=I_n$ is similar to a matrix in the form $diag(I_p,-I_q)$; then a restriction $f_{|U}$ is a local diffeomorphism with values in a neighborhood of $I_n$ iff $U$ is a neighborhood of $\pm I_n$.
Conclusion. If $B$ is close to $I_n$ then there are at least two square roots of $B$, exactly ONE close to $I_n$ and exactly ONE another close to $-I_n$.
EDIT. (b) On the other hand, there are $B$ close to $I_n$ that have no other square roots than those above.
Consider $B=I_n+aJ_n$ where $a$ is a non-zero small real number and $J_n$ is the nilpotent Jordan block of dimension $n$ and solve $X^2=I+aJ_n$. Since $X$ and $J_n$ commute, $X$ is in the form $x_0I+x_1J_n+x_2J_n^2+\cdots+x_{n-1}J_n^{n-1}$. Since the system $I,J_n,J_n^2,\cdots,J_n^{n-1}$ is free, we obtain $x_0=\pm 1=\epsilon,x_1=\epsilon a/2,x_2=-\epsilon a^2/8\cdots$. Finally, there are exactly $2$ solutions $X=\pm (I+\dfrac{a}{2}J_n-\dfrac{a^2}{8}J_n^2+\cdots)$ (you can recognize the formal Taylor development of $(1+aJ_n)^{1/2}$).