Csdrhrt

I am trying to invert a matrix using Woodbury identity. The inversion using Cholesky decomposition has the following pseudo-code:

For $t=1,2,...$

$(1);; text{Read};x_tinmathbb{R}^n$

$(2);;D_{t-1}=diag(|theta_t^1|,...,|theta_t^n|)$

$(3);;A_t=A_t+x_tx_t'$

$(4);;A_{t}^{-1}=sqrt{D_{{t-1}}}left(amathbf{I}+sqrt{D_{{t-1}}}A_tsqrt{D_{{t-1}}}right)^{-1}sqrt{D_{{t-1}}}$

$(5);;text{Read};y_tinmathbb{R}$

$(6);;b=b+y_tx_t$

$(7);;theta_t=A_{t}^{-1}b$

End For

The well known application of Sherman-Morrison on Recursive Least Squares is as follows:

$$A_t^{-1}=(aI+x_tx_t')^{-1}=A_{t-1}^{-1} - frac{(A_{t-1}^{-1}x_t)(A_{t-1}^{-1}x_t)'}{1+x_t'A_{t-1}^{-1}x_t}$$

where $A_0^{-1}=frac{1}{a}I$ and we can set $A_t = A_{t-1}+sum_{t=1}^Tx_tx_t'$, which will lead to time complexity of $O(n^2)$.

The above technique is mentioned here. The two implementation in $texttt{R}$ are as follows:

X <-matrix(runif(1000),20,10)

Y<-rnorm(20)

a<- 0.1



Cholsky<-function(X,Y,a){

  X <- as.matrix(X)

  Y <- as.matrix(Y)

  T <- nrow(X)

  N <- ncol(X)

  aI<- diag(a,N)

  bt<- matrix(0,ncol=1,nrow=N)

  for (t in 1:T){

    xt<-X[t,]

    At <- aI + (xt %*% t(xt))

    InvA<-chol2inv(chol(At))

    bt <- bt + (Y[t] * xt)

    theta<- InvA %*% bt

  }

  return(theta)

}

Cholsky(X,Y,a)





Morrison<-function(X,Y,a){

  X <- as.matrix(X)

  Y <- as.matrix(Y)

  T <- nrow(X)

  N <- ncol(X)

  At<-diag(1/a,N)

  bt<- matrix(0,ncol=1,nrow=N)

  for (t in 1:T){

    xt<-X[t,]

    At <- At + (xt %*% t(xt))

    InvA <- At - ((t(xt%*%At)%*%(as.matrix(xt%*%At)))

                    /as.numeric(xt%*%At%*%xt+1))

    bt <- bt + (Y[t] * xt)

    theta<- InvA %*% bt

  }

  return(theta)

}

Morrison(X,Y,a)

They don't give the same result. So, perhaps I should not expect the implementations to be equivalent.

I was wondering if I could invert the following (for the above case) more efficiently:

$$A_t^{-1}=left(D_{t-1}^{frac{1}{2}}A_tD_{t-1}^{frac{1}{2}}+aIright)^{-1}$$

where $A_t=A_{t-1}+x_tx_t'$ and $A_0=mathbf{0}$.

Essentially, I want to invert:

$$M_t=left(amathbf{I}+sqrt{D_{{t-1}}}x_tx_t'sqrt{D_{{t-1}}}right)$$

I say:

$$M_{t}^{-1}=M_{t-1}^{-1}-frac{(M_{t-1}D_{t-1}^{frac{1}{2}}x_t)(M_{t-1}D_{t-1}^{frac{1}{2}}x_t)'}{1+x_t'D_{t-1}^{frac{1}{2}}M_{t-1}D_{t-1}^{frac{1}{2}}x_t}$$

where $D_0=diag(mathbf{1}$). So, $M^{-1}_0=frac{1}{a}I$

RLS identity given above $(aI+u_tv_t)^{-1}$ uses $u=A_{t-1}^{-1}x_t,v_t=u_t'$,
I am using $u=M_{t-1}^{-1}D_{t-1}^{frac{1}{2}}x_t,v_t=u_t'$

One may write the implementations in $texttt{R}$ as follows:

X <-matrix(runif(1000),20,10)

Y<-rnorm(20)

a<- 0.1



#Cholesky implementation

X <- as.matrix(X)

Y <- as.matrix(Y)

T <- nrow(X)

N <- ncol(X)

bt<- matrix(0,ncol=1,nrow=N)

At<- diag(0,N)

I<- diag(a,N);Mt<-diag(1/a,N)

theta0<- rep(1,N)

for (t in 1:2){

  xt<-X[t,]

  Dt <- diag(sqrt(abs(as.numeric(theta0))))

  At <- At + (xt %*% t(xt))

  Mt <-  I + (Dt%*%At%*%Dt)

  InvA <- chol2inv(chol( Mt )) 

  AAt<- Dt %*%InvA%*% Dt

  bt <- bt + (Y[t] * xt)

  theta0 <- AAt %*% bt

  print(theta0)

}

Above is the correct implementation of the pseudo code. If I swap the following lines. I don't get the same answer.

Mt <-  Mt + (Dt%*%At%*%Dt)

InvA<- Mt - ((t(xt%*%Dt%*%Mt)%*%(as.matrix(xt%*%Dt%*%Mt)))

                   /as.numeric(xt%*%Dt%*%Mt%*%t(as.matrix(xt%*%Dt))+1))

Why is that?