我只是好奇

Operations on Hermitian Matrix

发表于 分类于
本文字数： 4.1k 阅读时长 ≈ 4 分钟

本篇文章拖了一段时间，说一下几个常见的复矩阵操作问题。在传统的信号处理方法中，频域是定义在复数域（$\mathbb{C}$）中的，所以我们经常会在算法中看到复数矩阵上的诸如求逆，求行列式，特征值分解等操作，比如adaptive beamforming（求逆，特征值），WPE（求逆）和CGMM（求逆，行列式）等等。最近几年，基于NN的前后端联合建模这块大家做的比较多，多通道上要结合自适应波束，WPE，或者使用一些最大似然估计的无监督损失函数，就需要将上述操作包络到网络中参与训练。因此如何在仅支持实数操作的DL框架（比如PyTorch）中进行这些复矩阵操作且保证梯度可传就成为了一个首先要被解决的问题。

目前最直观和方便的思路就是通过实数矩阵得到复矩阵操作的等价结果，像一些常见的加减乘除操作计算起来都十分直接。本篇文章主要说一下求逆，特征值，行列式这三个略微复杂的计算方法。一般情况下，我们进行这些操作的复矩阵对象为埃尔米特矩阵（满足共轭对称），因此可以充分利用此矩阵的相关性质简化计算。令埃尔米特矩阵$H \in \mathbb{C}^{N \times N}$，讨论如下：

求逆

假设$H$存在可逆矩阵$P$，根据定义有：

据此构建出的实矩阵$\bar{H}, \bar{P} \in \mathbb{R}^{2N \times 2N}$：

同样满足$\bar{H}*\bar{P} = I$：

即如果$H$可逆，那么其逆矩阵的实部和虚部可由其拓展实矩阵$\bar{H}$的逆$\bar{P} = \bar{H}^{-1}$求出。在PyTorch框架上，只需要调用一次对$\bar{H}$的求逆即可。

特征值 & 行列式

根据埃尔米特矩阵的性质，其特征值都为实数，加之前端算法中处理的相关矩阵一般都要求其可逆，因此$H$也满足正定性，即特征值为正。前端的一些算法中其实对特征值这个概念没有直接的需求，波束里面多用的是特征向量，但是有些情况下我们需要求解$H$矩阵的行列式的值，什么情况下呢，比如我们要优化cacgmm的目标函数：

$\mathbf{R}_f^{-1}$项我们已经解决，剩下$\det \mathbf{R}_f$待解。根据行列式等于特征值的乘积的性质可将行列式的计算转化为特征值的求解问题。假设$\lambda$是$H$的一个特征值，$p$是其对应的特征向量，根据定义有：

对于$(2)$中的$\bar{H}$，有：

以及

即$\lambda$是$\bar{H}$的二重特征值，其对应的两个特征向量为$[p_r, p_i]$和$[-p_i, p_r]$。如果只是想求行列式，那么可以暂且不关心特征值和特征向量的对应关系，得到$N$个非重特征值$\lambda_{0 \cdots N - 1}$之后，由：

得到行列式的值即可。

特征向量

如果想进一步得到特征向量的表示，需要看一下其和特征值之间的对应关系。将$(6)(7)$式合并写成矩阵形式：

其中$\Lambda = \text{diag}(\lambda_0, \cdots, \lambda_{N - 1})$。在numpy和pytorch的eig函数得到的特征值默认将重值放在一起，因此需要对应的把$(9)$式改写成

其中

不过由于存在重根，因此$\bar{P}’$矩阵中的第$i \times N$列和$i \times N + 1$列即使发生交换，$(10)$式也同样成立，这样就带来了实部，虚部和正负号上的混淆。我暂时还没有找到一种确定性的方法能够完全匹配复矩阵下的特征向量计算结果，所以也就没有将特征向量的计算纳入到神经网络的优化中去。后续如果发现好的方案继续在此更新。

测试

最后可以在numpy上简单写个测试程序验证一下上述结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#!/usr/bin/env python

# wujian@2020

import numpy as np
import torch as th

def random_cmat(N):
    R = np.random.rand(N, N)
    I = np.random.rand(N, N)
    M = R + I * 1j
    M = 0.5 * (M + M.T.conj())
    eig_val, eig_vec = np.linalg.eigh(M)
    eig_val[eig_val < 0.0001] = 0.0001
    M = eig_vec @ np.diag(eig_val) @ eig_vec.T.conj()
    return M
  
def pack_cmat(cmat):
    n = cmat.shape[0]
    rmat = np.zeros([n * 2, n * 2])
    rmat[:n, :n] = cmat.real
    rmat[:n, n:] = -cmat.imag
    rmat[n:, :n] = cmat.imag
    rmat[n:, n:] = cmat.real
    return rmat

def test_cplx_inv():
    M = random_cmat(4)
    M_inv = np.linalg.inv(M)
    N = pack_cmat(M)
    N_inv_fake = np.linalg.inv(N)
    N_inv = np.zeros_like(M)
    N_inv.real = N_inv_fake[:4, :4]
    N_inv.imag = N_inv_fake[4:, :4]
    assert np.allclose(M_inv, N_inv)

def test_cplx_eig():
    M = random_cmat(4)
    cplx_eig_val, _ = np.linalg.eigh(M, UPLO="U")
    cplx_det = np.linalg.det(M).real
    N = pack_cmat(M)
    real_eig_val, _ = np.linalg.eigh(N, UPLO="U")
    real_eig_val = real_eig_val[::2]
    real_det = np.cumprod(real_eig_val)[-1]
    assert np.allclose(cplx_eig_val, real_eig_val)
	assert np.allclose(cplx_det, real_det)
    
if __name__ == "__main__":
    test_cplx_inv()
    test_cplx_eig()