Para utilizar la transformada de Fourier Discreta para invertir una matriz de covarianza

Question

Estoy trabajando en un problema que en su parte más difícil es para invertir una matriz de covarianza (en R). No pude uso habitual approches como la enfermedad vesicular porcina y Col. Entonces, decidí usar una transformada de Fourier Discreta (DFT). Pero yo no podía entender cómo aplicar el método, especialmente en R. Así, sus comentarios y posibles ejemplos en R, se agradece. muchas gracias de antemano.

Answer 1

Un circulantes es una matriz cuya primera columna es un vector $x$ y su posterior columnas se obtiene por la rotación de un elemento a la derecha. Aquí es el código R para producir cualquier circulantes de su primera columna 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

Por ejemplo,

> 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

Es invertida por el cambio a un eigenbasis. Los elementos de la diagonal son las entradas de la transformada de Fourier de x, así que invertir de forma individual y cambiar de nuevo a la base original:

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)

Por ejemplo, podemos demostrar que esto funciona por la multiplicación de su salida por el original circulantes y la comprobación de que la identidad de la matriz se obtiene (hasta insignificante de punto flotante de error):

> 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

---

Ser conscientes de que la enfermedad acondicionado todavía plaga este enfoque debido a de punto flotante de redondeo en fft. Es por eso que he implementado un reciprocal función: se niega a calcular $1/x$ al $x=0$. Como tal, inverse.circulant calcula un generalizada inversa, exactamente como en 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

Answer 2

Es el artículo de Wikipedia sobre matrices circulantes claro? Esto es algo que en un momento de la serie contexto se discute, por ejemplo, en las primeras páginas de Hannan de la Serie de Tiempo del libro. Los autovalores de una circulantes de la matriz está dada por la transformada de Fourier de lo que (de nuevo en una serie de tiempo de contexto), que sería el autocovariances. Para invertir la matriz, usted tiene que tomar los recíprocos de los valores y de pre - y post-multiplicar por la matriz cuyas columnas son los vectores propios.

Answer 3

¿Has probado una corrección por la adición de una pequeña $\epsilon$ de perturbación a la diagonal de la matriz que están tratando de invertir. Este es un estándar de procesamiento de rutina se utiliza para aplazar la singularidad problema, mientras que la inversión de una matriz de covarianza o de un estado de hesse.