使用离散傅里叶变换反转协方差矩阵

机器算法验证 r 逆矩阵
2022-03-21 23:52:13

我正在研究一个问题,它的困难部分是反转协方差矩阵(在 R 中)。我不能使用像 SVD 和 Chol 这样的常用方法。然后,我决定使用离散傅里叶变换 (DFT) 方法。但我不明白如何应用该方法,尤其是在 R 中。因此,您的评论和 R 中可能的示例将不胜感激。提前谢谢了。

3个回答

循环是一个矩阵,其第一列是一个向量x其后续列是通过向右旋转一个元素获得的。这是从第一列生成任何循环的 R 代码x

rotate <- function(x,k) {c(tail(x,-k), head(x,k))}
circulant <- function(x) {
    n=length(x)
    apply(matrix(0:(n-1),1,n), 2, function(k) rotate(x,n-k))
} # Returns the circulant matrix of which x is the first column

例如,

> circulant(c(2,3,5,7))
     [,1] [,2] [,3] [,4]
[1,]    2    7    5    3
[2,]    3    2    7    5
[3,]    5    3    2    7
[4,]    7    5    3    2

它通过更改为特征基来反转。对角线元素是 的傅里叶变换的入口x,因此我们将它们单独反转并变回原来的基:

reciprocal <- function(x) {i <- which(x!=0); x[i] <- 1/x[i]; x}
inverse.circulant <- function(x) {
    n <- length(x)                  # x is the first column of the circulant
    i <- (0:(n-1)) %o% (1:n)        # Powers of exp(2 Pi I/n) in the eigenbasis q
    q <- matrix(exp(complex(real=-log(n)/2, imaginary=2*pi*i / n)), n, n)
    w <- reciprocal(fft(x))         # Reciprocals of nonzero eigenvalues
    Re(t(q) %*% diag(w) %*% Conj(q))# Convert back to the original basis
} # Returns a generalized inverse to circulant(x)

例如,我们通过将其输出乘以原始循环并检查是否获得了单位矩阵来证明这项工作(直至可忽略的浮点误差):

> a <- c(2,3,5,7)
> zapsmall(inverse.circulant(a) %*% circulant(a))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

请注意,由于fft. 这就是我实现一个reciprocal函数的原因:它拒绝计算1/x什么时候x=0. 因此,inverse.circulant计算广义逆,完全如MASS::ginv

# The following determines a nonsingular but ill-conditioned circulant:
> (a <- c(1, -200000/200001, -2500000/500001, 5000000/1000003))
[1]  1.000000 -0.999995 -4.999990  4.999985

> 1 / rcond(circulant(a)) # HUGE condition number!
[1]  5.404306e+16

> library(MASS)
> inverse.circulant(a) - ginv(circulant(a))
             [,1]          [,2]          [,3]          [,4]
[1,] 6.938894e-18 -2.081668e-17 -4.163336e-17  2.775558e-17
[2,] 1.387779e-17 -6.938894e-18 -3.469447e-18  1.387779e-17
[3,] 1.387779e-17  4.163336e-17 -6.938894e-18  1.040834e-17
[4,] 3.469447e-17 -2.775558e-17 -1.387779e-17 -2.081668e-17

您是否尝试过更正 - 通过添加一个小的ϵ对要反转的矩阵的对角线的扰动。这是一个标准处理例程,用于在反转协方差矩阵或粗麻布矩阵时推迟奇异性问题。

关于循环矩阵的维基百科文章是否清楚?这是在时间序列上下文中讨论的内容,例如,在Hannan 的时间序列书的第一页中。循环矩阵的特征值由自协方差(同样在时间序列上下文中)的傅里叶变换给出。因此,要反转矩阵,您必须取特征值的倒数,然后乘以列是特征向量的矩阵。