numpy 可以用实算术对角化斜对称矩阵吗？

Question

任何斜对称矩阵（A^T = -A）都可以变成厄米矩阵（iA) 并用复数对角化。但它是 also possible to bring it into block-diagonal form with a special orthogonal transformation 并且仅使用实算术找到它的特征值。这是在 numpy 中的任何地方实现的吗？

Answer 1

我们来看看numpy库的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()

dgeev()使用的方法在documentation of LAPACK. It requires the reduction of the matrix A to its real Schur formS中有描述。

任意实方阵A可以表示为：

A=QSQ^t

其中：

Q是实正交矩阵：QQ^t=I
S是实分块上三角矩阵。 S对角线上的块大小为1×1或2×2。

确实，如果 A 是斜对称的，则此分解似乎非常接近 A 的 block diagonal form obtained by a special orthogonal transformation。而且，真的是看到斜对称矩阵A的Schur形式S是...斜对称的！

确实，让我们计算 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

因此，若Q为特殊正交（det(Q)=1），则S为特殊正交变换得到的分块对角线形式。否则，一个特殊的正交矩阵 P 可以通过置换 Q 的前两列来计算，并且矩阵 A 的另一个 Schur 形式 Sd 通过改变符号来获得S_{12} 和 S_{21}。的确，A=PSdP^t。那么，Sd就是A经过特殊正交变换得到的分块对角形式

最后，即使numpy.linalg.eigvals()应用于实数矩阵returns复数，过程中涉及的复杂计算也很少！

如果您只想计算真正的 Schur 形式，请使用带参数 output='real' 的函数 scipy.linalg.schur()。

只是一段代码来检查：

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

Answer 2

是的，你可以通过在乘积中间加入单一变换来做到这一点，因此我们得到

A = V * U * V^-1 = V * T' * T * U * T' * T * V^{-1}.

一旦你有了想法，你就可以通过平铺来优化代码，但让我们通过明确地形成 T 来以天真的方式来做。

如果矩阵大小是偶数，则所有块都是复共轭。否则我们得到一个零作为特征值。特征值保证实部为零，所以首先要清除噪声，然后进行排序，使零点位于左上角（任意选择）。

n = 5
a = np.random.rand(n,n)
a=a-np.transpose(a)
[u,v] = np.linalg.eig(a)

perm = np.argsort(np.abs(np.imag(u)))
unew = 1j*np.imag(u[perm])

显然，我们也需要对特征向量矩阵重新排序以保持等价。

vnew = v[:,perm]

到目前为止，除了对特征值分解中的中间特征值矩阵进行重新排序外，我们什么也没做。现在我们从复数形式切换到实块对角线形式。

首先我们要知道有多少个零特征值

numblocks = np.flatnonzero(unew).size // 2
num_zeros = n - (2 * numblocks)

然后我们基本上形成另一个酉变换（这次是复杂的）并以同样的方式粘贴

T = sp.linalg.block_diag(*[1.]*num_zeros,np.kron(1/np.sqrt(2)*np.eye(numblocks),np.array([[1.,1j],[1,-1j]])))

Eigs = np.real(T.conj().T.dot(np.diag(unew).dot(T)))
Evecs = np.real(vnew.dot(T))

这为您提供了新的实值分解。所以代码都在一个地方

n = 5
a = np.random.rand(n,n)
a=a-np.transpose(a)
[u,v] = np.linalg.eig(a)

perm = np.argsort(np.abs(np.imag(u)))
unew = 1j*np.imag(u[perm])
vnew = v[perm,:]
numblocks = np.flatnonzero(unew).size // 2
num_zeros = n - (2 * numblocks)
T = sp.linalg.block_diag(*[1.]*num_zeros,np.kron(1/np.sqrt(2)*np.eye(numblocks),np.array([[1.,1j],[1,-1j]])))
Eigs = np.real(T.conj().T.dot(np.diag(unew).dot(T)))
Evecs = np.real(vnew.dot(T))
print(np.allclose(Evecs.dot(Eigs.dot(np.linalg.inv(Evecs))) - a,np.zeros((n,n))))

给出True。请注意，这是获得真实频谱分解的朴素方法。有很多地方需要跟踪数值误差的累积。

示例输出

Eigs
Out[379]: 
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        , -0.61882847,  0.        ,  0.        ],
       [ 0.        ,  0.61882847,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        , -1.05097581],
       [ 0.        ,  0.        ,  0.        ,  1.05097581,  0.        ]])

Evecs
Out[380]: 
array([[-0.15419078, -0.27710323, -0.39594838,  0.05427001, -0.51566173],
       [-0.22985364,  0.0834649 ,  0.23147553, -0.085043  , -0.74279915],
       [ 0.63465436,  0.49265672,  0.        ,  0.20226271, -0.38686576],
       [-0.02610706,  0.60684296, -0.17832525,  0.23822511,  0.18076858],
       [-0.14115513, -0.23511356,  0.08856671,  0.94454277,  0.        ]])

numpy 可以用实算术对角化斜对称矩阵吗？

Can numpy diagonalise a skew-symmetric matrix with real arithmetic?

python

numpy

matrix

linear-algebra

lapack