What is your use for the root?  

For square matrices, many function can be defined using the Jordan canonical form of the matrix.  But not for your question.  

Checked the book by C. T. Chen.  

Article Linear system theory and design, by Chi-Tsong Chen, Oxford U...

You should not seek to apply your formula to matrix directly! Norm for matrices do exist though, as well as exponential and so on, but that would be a bit awkward.

Working on the set of complex numbers, you can readily demonstrate an isomorphism with the space of matrix of size 2x2. Your n roots do not define a basis, indeed.

Now, interpreting your roots as actual rotations of angle Arg(z)/n between successive roots, it is easy to define your set of rotation matrices (Think about the representation of the roots of unity on the circle). Their norm will be 1. Simply multiply each matrix by det(Z)^{1/n} to have matrices representing the n-th root of Z.

 Dear Alberto A. Díaz,

By using Spectral Decomposition. Let

A= U^{*}\diag(\lambda}{1}), \dots, \lambda{n})U,

where $U$ is unitary, and $\lambda_{i} \ge 0$ for all $i$.  We take

B = U^{*}\diag(\lambda_{1}^{1/2}},\dots, \lambda_{n}^{1/2})U,

Then the matrix $B$ is positive semidefinite, and B^{2} = A since

U^{*}U = I_{n}.

See for example:

W.-H. Steeb, Problems and solutions in introductory and advanced matrix calculus, World Scientific Publishing Co. Pte. Ltd., Toh Tuck Link, Singapore, 2006. pp.89-90

Dear Alberto,

Let A be a selfadjoint operator on a Hilbert space, with the spectrum contained in (0, infinity) and p>1 (not necessarily an integer). The problem is to solve the equation U^p -A = 0, in the unknown U. This can be done by using Newton's approximating method, for convex monotone operators (see my paper "Newton's method for convex operators and applications" on RG).

Best wishes, Octav