numpy可以用实数对角化斜对称矩阵吗？答案

【问题标题】：Can numpy diagonalise a skew-symmetric matrix with real arithmetic?numpy可以用实数对角化斜对称矩阵吗？
【发布时间】：2016-08-01 20:41:20
【问题描述】：

任何斜对称矩阵（A^T = -A）都可以转化为厄米特矩阵（iA) 并用复数对角化。但将其带入block-diagonal form with a special orthogonal transformation 并仅使用实数算术找到其特征值是also possible。这是在 numpy 的任何地方实现的吗？

【问题讨论】：

你为什么想要那个？无论哪种方式，特征值都是相同的。

标签： python numpy matrix linear-algebra lapack

【解决方案1】：

我们来看看LAPACK库的dgeev()函数。此例程计算任何实双精度方阵的特征值。而且，这个例程就在numpy库的python函数numpy.linalg.eigvals()后面。

dgeev() 使用的方法在documentation of LAPACK 中有描述。它需要将矩阵 A 减少到它的 real Schur form S。

任何实数方阵A可以表示为：

A=QSQ^t

地点：

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

确实，如果A 是斜对称的，那么这种分解似乎非常接近 block diagonal form obtained by a special orthogonal transformation 的A。而且，确实看到斜对称矩阵A的舒尔形式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是通过特殊正交变换得到的块对角形式。否则，可以通过置换Q 的前两列来计算一个特殊的正交矩阵P，并通过改变S_{12} 和S_{21} 的符号得到矩阵A 的另一个舒尔形式Sd。确实，A=PSdP^t。那么Sd就是A通过特殊正交变换得到的块对角形式。

最后，即使numpy.linalg.eigvals() 应用于实数矩阵返回复数，过程中也几乎没有复杂的计算！

如果您只想计算真正的 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

【讨论】：

【解决方案2】：

是的，您可以通过在产品中间进行统一转换来做到这一点，因此我们得到 p>

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.        ]])

【讨论】：