Trench Inverse Update

Initialization

<< FGN.m

For Debugging R Version

r = Table[FGNAcf[k, 0.7], {k, 0, 5}]

(GI5 = Inverse[toe[Take[r, 5]]]) // MatrixForm

(GI6 = Inverse[r // toe]) // MatrixForm

r5 = Take[r, 5]

rnew = r[[6]]

g = {r5, rnew} // Flatten // Reverse // Most

v = First[r5]

e = 1/(v - g . GI5 . g)

f = -GI5 . g*e

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

Inverse of Partitioned Matrix

So for updating the Toeplitz matrix we can write,

(((({{Γ, g}, {(g^′), v}})))^(-1))=(({{(Γ^(-1))(
                  I + (Γ^(-1))/a), -h}, {(h^′), b}}))

r = Table[FGNAcf[k, 0.7], {k, 0, 10}];

r5 = Take[r, 5];

R5 = toe[r5];

R5Inv = Inverse[R5];

GI = R5Inv; G = R5;

G2 = toe[Take[r, 6]];

MatrixForm[G2]

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

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

MatrixForm[Inverse[toe[Take[r, 6]]]]

e

fok = f = -GI . g*e

<< LinearAlgebra`MatrixManipulation`

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

Partition[{Partition[Flatten[{A, f}], Length[f]] // Transpose, f, e} // 
          Flatten, Length[f] + 1] // MatrixForm

Code Development

rA = Table[FGNAcf[k, 0.7], {k, 0, 6}];
r = Take[rA, 5];
GI = Inverse[G = toe[r]];
rnew = 
        rA[[6]];
MatrixForm[Inverse[toe[Take[rA, 6]]]]

g = {r, rnew} // Flatten // Reverse // Most;
v = G[[1, 1]];
e = 1/(v - g . GI . g);
a = g . g*e;
A = GI + (GI . Outer[Times, g, g] . GI)*e;
f = -GI . g*e;
H = BlockMatrix[{{A, Partition[f, 1]}, {{f}, {{e}}}}];

MatrixForm[H]

Using Compiler

Fast22 = 
        Compile[{{GI, _Real, 2}, {v, _Real}, {g, _Real, 
              1}}, 
1/(v - g . GI . g)];

e = Fast22[GI, v, g]

Fast12 = 
        Compile[{{GI, _Real, 2}, {e, _Real}, {g, _Real, 1}}, -GI . g*
            e];

f = Fast12[GI, e, g]

fok - f // Abs // Total

Fast11 = 
        Compile[{{GI, _Real, 2}, {e, _Real}, {g, _Real,  
              1}}, 
GI + (GI . 
                  Table[g[[i]]*
                      g[[j]], {i, 
                      Length[g]}, {j, Length[g]}] . GI)*e];

A = Fast11[GI, e, g]

A - Aok // Flatten // Abs // Total

H = 
      Partition[{Partition[Flatten[{A, f}], Length[f]] // Transpose, f, e} // 
          Flatten, Length[f] + 1]

H - Hok // Flatten // Abs // Total

Mathematica Function

\(\(TrenchInverseUpdate::usage = 
"TrenchInverseUpdate[GI,r,rnew] uses partitioned matrix result to obtain the order n+1 inverse given the order n inverse of G, GI; the first row of G, r; and the next element of the first row rnew.";\)\)

FastGIg = 
        Compile[{{GI, _Real, 2}, {g, _Real, 1}}, g . GI];

Fast22 = 
        Compile[{{GIg, _Real, 1}, {g, _Real, 
              1}, {v, _Real}}, 
1/(v - 
                Total[g*GIg])];

Fast12 = 
        Compile[{{GIg, _Real, 1}, {e, _Real}}, -GIg*e];

Fast11 = 
        Compile[{{GI, _Real, 2}, {e, _Real}, {GIg, _Real,  
              1}}, 
GI + 
            Table[GIg[[i]]*
                  GIg[[j]], {i, 
                  Length[GIg]}, {j, Length[GIg]}]*e];

TrenchInverseUpdate[GI_, r_, rnew_] := 
      Module[{v = First[r], 
          g = {r, rnew} // Flatten // Reverse // Most, e, f, 
          A}, 
GIg = FastGIg[GI, g]; 
e = 
          Fast22[GIg, g, v]; 
f = 
          Fast12[GIg, e]; 
A = 
          Fast11[GI, e, 
            GIg]; 
Partition[{Partition[Flatten[{A, f}], 
                  Length[f]] // Transpose, f, e} // Flatten, 
          Length[f] + 1]
]

Save Package

DeleteFile["d:/math/2006/Trench/code/update.m"]

\(Save["d:/math/2006/Trench/code/update.m", {TrenchInverseUpdate}]\)

Testing

rA = Table[FGNAcf[k, 0.7], {k, 0, 6}];
r = Take[rA, 5];
GI = Inverse[G = toe[r]];
rnew = 
        rA[[6]];
MatrixForm[Inverse[toe[Take[rA, 6]]]]

TrenchInverseUpdate[GI, r, rnew] // MatrixForm

Flatten[
        Inverse[toe[Take[rA, 6]]] - TrenchInverseUpdate[GI, r, rnew] // Abs] // 
      Total

Timing

rA = Table[FGNAcf[k, 0.7], {k, 0, 1000}]; n = 500;
GI2 = GI = TrenchInverse[toe[Take[rA, 500]]];
t1 = t2 = 0;
err = {};
Do[
G = toe[Take[rA, n + i]]; 
t1 += 
        First[Timing[GI1 = TrenchInverse[G]; ]]; 
GI = 
        GI2; 
r = 
        Take[rA, n + i - 1]; 
rnew = 
        rA[[n + 
            i]]; 
t2 += 
        First[Timing[
            GI2 = TrenchInverseUpdate[GI, r, 
                rnew]; ]]; 
err = {err, Flatten[
                GI1 - GI2] // Abs // 
            Total};
, 
{i, 
        100}]

{t1, t2}

err = Flatten[err]; {Mean[#], Max[#], Min[#]} &[err]