Toeplitz Inverse Updating

Introduction

The general result for partitioned matrices

Toeplitz case

So for updating the Toeplitz matrix from (Subscript[Γ, n]) to (Subscript[Γ, n + 1]) we can write,

(Subscript[Γ, n + 1]^(-1))= (((({{(Subscript[Γ, n]^(-1)), g}, {(g^′), γ(0) }})))^(-1))=(({{A, -f}, {(f^′), e}}))

where g,  v,  f,  e and A are defined in the Mathematica code below.

\(TraditionalForm\`Γ\_n = \((γ\_\(i - j\))\)\_\(n×n\)\)

\(TraditionalForm\`g = \((γ(n), \ γ(n - 1), \  ... \ , \ γ(1))\)\^′\)

e = 
      1/(γ(0) - 
            (g^′) (Subscript[Γ, n]^(-1)) g)

f = -(Subscript[Γ, n]^(-1)) g e

\(TraditionalForm\`A = Γ\_n\%\(-1\)\ \  + \ Γ\_n\%\(-1\)\ g\ g\^′\ Γ\_n\%\(-1\)\ e\)

Efficient Computation

For efficiency we just need to compute (α^′) =  
        (g^′) (Subscript[Γ, n]^(-1)) and then

e = 1/(γ(0) - (α^′) g)

f = -g  e

\(TraditionalForm\`A = Γ\_n\%\(-1\)\ \  + \ α\ α\^′\ e\)

Implementing this in my R drastically increased the speed.

It is interesting that updating requires O((n^2)) flops but nevertheless is usually much faster than the Trench algorithm.

Mathematica Code

This is just for prototyping.  I wrote an actual Mathematica function as well but this explains better the steps involved.  The code does not make use of the efficient method discussed above but is just a direct implementation.

<< LinearAlgebra`MatrixManipulation`

<< FGN.m

\(\(\(\(r = Table[FGNAcf[k, 0.9], {k, 0, 5}];\)

    \(r5 = Take[r, 5];\)

    \(R5 = toe[r5];\)

    \(R5Inv = Inverse[R5];\)

    \(\(GI = R5Inv\) \); \  (*inverse\ n = 5\ *) 
\(\(G = 
        R5\) \); \  (*\ 
      original\ matrix\ *) 
\(\(GI2 = 
        Inverse[toe[Take[r, 6]]]\) \);\)\(\ \ \ \)\)\)

r

The inverse for n = 6 obtained directly is

MatrixForm[GI2]

Now using the partitioned matrix formula

g = {r5, r[[6]]} // Flatten // Reverse // 
          Most;
v = G[[1, 1]];
e = 1/(v - g . GI . g);

MatrixForm[A = GI + (GI . Outer[Times, g, g] . GI)*e]

e

f = -GI . g*e

(H = BlockMatrix[{{A, Partition[f, 1]}, {{f}, {{e}}}}]) // 
      MatrixForm

Total[Flatten[Abs[H - GI2]]]

R Implementation

R Function

"ToeplitzInverseUpdate" <-
function( GI, r, rnew )
{   
    g <- rev(c(r,rnew)[-1])
    GIg <- c(crossprod(GI,g))
    e <- 1/(r[1]-sum(g*GIg))
    f <- -GIg*e
    A <- GI + outer(GIg,GIg)*e
    rbind(cbind(A,f),c(f,e))
}

Note:  GI is (Subscript[Γ, n]^(-1)), r is the first row of (Subscript[Γ, n]) and rnew is the last element of the first row of (Subscript[Γ, n + 1]), ie. rnew=γ(n)

Example Usage

Evaluating this script shows the correctness of the function

r<-FGNAcf(0:5,0.9)
r2<-r[1:5]
G<-toe(r2)
rnew<-r[6]
GI<-TrenchInverse(G)
GI2<-ToeplitzInverseUpdate(GI, r2, rnew)
GI2Direct<-TrenchInverse(toe(r))
err<-sum(abs(GI2-GI2Direct))

> err
[1] 1.135203e-14

Timing

The following script evaluates (Subscript[Γ, n + 1]^(-1)) for n = 501,   ...  ,  1000 and computes the total CPU time for the update and just using the direct method.  The updating method is faster by a factor of 3.

tim1<-tim2<-0
NREP<-500
err<-numeric(NREP)
rTOT<-FGNAcf(0:1000, 0.7)
n<-500
r<-rTOT[1:n]
G<-toe(r)
GI<-TrenchInverse(G)
for (i in 1:NREP){
  rnew<-rTOT[n+i]
  tim1<-tim1+system.time(GI2<-ToeplitzInverseUpdate(GI, r, rnew))[1]
  r<-c(r,rnew)
  G<-toe(r)
  tim2<-tim2+system.time(GI3<-TrenchInverse(G))[1]
  GI<-GI2
  err[i]<-sum(abs(GI2-GI3))
}
tottim<-(tim1+tim2)
rm(G,GI,GI2,GI3,r,rTOT)

> tim1
[1] 26.04
> tim2
[1] 90.06

The accuracy is good in this example.

> summary(err)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
5.798e-13 2.312e-12 3.225e-12 3.284e-12 4.252e-12 6.231e-12

As might be expected the accuarcy of the updating does slightly decrease as the number of updates increases

> apply(matrix(err, ncol=5), MARGIN=2, mean)
[1] 1.541542e-12 2.573623e-12 3.219167e-12 4.065311e-12 5.019128e-12