Any skew-symmetric matrix (A^T = -A) can be turned into a Hermitian matrix (iA) and diagonalised with complex numbers. But it is also possible to bring it into block-diagonal form with a special orthogonal transformation and find its eigevalues using only real arithmetic. Is this implemented anywhere in numpy?
Let's take a look at the dgeev() function of the LAPACK librarie. This routine computes the eigenvalues of any real double-precison square matrix. Moreover, this routine is right behind the python function numpy.linalg.eigvals() of the numpy library.
The method used by dgeev() is described in the documentation of LAPACK. It requires the reduction of the matrix A to its real Schur form S.
Any real square matrix A can be expressed as:
A=QSQ^t
where:
Q is a real orthogonal matrix: QQ^t=IS is a real block upper triangular matrix. The blocks on the diagonal of S are of size 1×1 or 2×2.Indeed, if A is skew-symmetric, this decomposition seems really close to a block diagonal form obtained by a special orthogonal transformation of A. Moreover, it is really to see that the Schur form S of the skew symmetric matrix A is ... skew-symmetric !
Indeed, let's compute the transpose of S:
S^t=(Q^tAQ)^t
S^t=Q^t(Q^tA)^t
S^t=Q^tA^tQ
S^t=Q^t(-A)Q
S^t=-Q^tAQ
S^t=-S
Hence, if Q is special orthogonal (det(Q)=1), S is a block diagonal form obtained by a special orthogonal transformation. Else, a special orthogonal matrix P can be computed by permuting the first two columns of Q and another Schur form Sd of the matrix A is obtained by changing the sign of S_{12} and S_{21}. Indeed, A=PSdP^t. Then, Sd is a block diagonal form of A obtained by a special orthogonal transformation.
In the end, even if numpy.linalg.eigvals() applied to a real matrix returns complex numbers, there is little complex computation involved in the process !
If you just want to compute the real Schur form, use the function scipy.linalg.schur() with argument output='real'.
Just a piece of code to check that:
import numpy as np
import scipy.linalg as la
a=np.random.rand(4,4)
a=a-np.transpose(a)
print "a= "
print a
#eigenvalue
w, v =np.linalg.eig(a)
print "eigenvalue "
print w
print "eigenvector "
print v
# Schur decomposition
#import scipy
#print scipy.version.version
t,z=la.schur(a, output='real', lwork=None, overwrite_a=True, sort=None, check_finite=True)
print "schur form "
print t
print "orthogonal matrix "
print z