Most scalar functions can be generalized to square matrices, which requires to apply the function to the eigenvalues

Square root of a matrix

#!/usr/bin/env python
#coding:utf8

import numpy as np
import scipy.linalg as la

size = 3
A = np.random.random([size, size])

print "\n== Random square matrix can be subject to virtually any function =="
print A

print "\n== Eigenvalue decomposition of A =="
n, V = la.eig(A)       
print "eigenvalues: ", n
print "eigenvectors:", V

print "\n== Square root of eigenvalues =="
print "eigenvalue square roots:", n**.5

print "\n==  Square root of a matrix,  B = V . sqrt(D) . inv(V)  =="
B = np.real(np.dot(V, np.dot(np.diag(n)**.5, la.inv(V)))) 
print B

print "\n== Proof: B . B - A  is almost zero =="
print np.dot(B, B) - A


#print "\n-- SVD --"
#U, d, Vh = la.svd(A)     ## U, Vh are unitary
#print "eig U", la.eig(U)
#print "A=", A
#print "U*D*Vh =", np.dot(U, np.dot(np.diag(d), Vh))

#print "\n-- Schur decomposition --"
#T,U = la.schur(A)      ## U unitary
#print np.dot(U, np.dot(T, U.T.conj())) #  NOT WORKING ??

With 3x3 matrices I observed that sometimes the test works, and sometimes not.

Match was only obtained when two of the three eigenvalues were complex conjugate.