To fast invert a real symmetric positive definite matrix that is almost similar to Toeplitz

Question

It is well known how to solve a Toeplitz system Ax = b, of a matrix A, n x n elements, through the use of circular embedding in 2n x 2n and use of fft and ifft algorithm (easy reference, fast toeplitz solution, good reference, fast algorithms). However I have a matrix that is slightly different from Toeplitz although it still has a diagonal-constant structure. Consider the matrix, (see the repetition in the 2nd column)

Special =
[ a, b, c, d, e, f, g, h, i]
[ b, a, b, e, d, e, h, g, h]
[ c, b, a, f, e, d, i, h, g]
[ d, e, f, a, b, c, d, e, f]
[ e, d, e, b, a, b, e, d, e]
[ f, e, d, c, b, a, f, e, d]
[ g, h, i, d, e, f, a, b, c]
[ h, g, h, e, d, e, b, a, b]
[ i, h, g, f, e, d, c, b, a]

toeplitz(Special(:,1)) =
[ a, b, c, d, e, f, g, h, i]
[ b, a, b, c, d, e, f, g, h]
[ c, b, a, b, c, d, e, f, g]
[ d, c, b, a, b, c, d, e, f]
[ e, d, c, b, a, b, c, d, e]
[ f, e, d, c, b, a, b, c, d]
[ g, f, e, d, c, b, a, b, c]
[ h, g, f, e, d, c, b, a, b]
[ i, h, g, f, e, d, c, b, a]

Now since the degree of freedom is still n , only n elements are needed to define the entire special matrix, we can therefore expect the solution would be easier and fast. The matrix is diagonalizable since it is symmetric positive definite, so it has real positive eigenvalues, so my question is how can I proceed to find the fast solution ?

I am not looking for a spoon-fed answer but a general direction/ hints/ comments as how one would go about solving such a problem.

I believe the derivation should be similiar to the toeplitz case, but in that case the solution was easy with circular embedding, in my case the matrix is not that amenable to be transformed into circulant matrix (atleast thats what it appears to me). I can also try MATLAB symbolic decomposition, but for large matrices with similar format (but not exactly) it would be difficult to get an answer.

Answer 1

Your matrix is block Toeplitz with Toeplitz blocks (BTTB), with $3\times 3$ blocks. I am not aware of direct algorithms for this kind of structure, but the natural thing to try is embedding it into a BCCB (block circulant with circulant blocks) matrix which can be inverted using a two-dimensional FFT.

Most of the literature I could find on BTTB focuses on the large case and iterative algorithms. If your matrix really is 9x9, I doubt that you can get any real speedup with respect to a good implementation of Gaussian elimination; this kind of algorithms usually starts paying off for larger dimensions.

Answer 2

This was answered by Federico Poloni. I am merely providing the details (some notations are changed). Consider a block toeplitz toeplitz blocks (BTTB) matrix (I am working on a signal processing problem where I encountered large matrices of this form to invert)of size 9 x 9 (i.e 3^2 x 3^2),

   BTTB =
    [ a_1, a_2, a_3, a_2, a_5, a_6, a_3, a_6, a_9]
    [ a_2, a_1, a_2, a_5, a_2, a_5, a_6, a_3, a_6]
    [ a_3, a_2, a_1, a_6, a_5, a_2, a_9, a_6, a_3]
    [ a_2, a_5, a_6, a_1, a_2, a_3, a_2, a_5, a_6]
    [ a_5, a_2, a_5, a_2, a_1, a_2, a_5, a_2, a_5]
    [ a_6, a_5, a_2, a_3, a_2, a_1, a_6, a_5, a_2]
    [ a_3, a_6, a_9, a_2, a_5, a_6, a_1, a_2, a_3]
    [ a_6, a_3, a_6, a_5, a_2, a_5, a_2, a_1, a_2]
    [ a_9, a_6, a_3, a_6, a_5, a_2, a_3, a_2, a_1]

Then the individual Toeplitz Blocks of size 3 x 3 are given by,

T1 =                     T2 =                       T3 = 
   [ a_1, a_2, a_3]          [ a_2, a_5, a_6]            [ a_3, a_6, a_9]
   [ a_2, a_1, a_2]          [ a_5, a_2, a_5]            [ a_6, a_3, a_6]
   [ a_3, a_2, a_1]          [ a_6, a_5, a_2]            [ a_9, a_6, a_3]

The Composite matrix (in terms of T1, T2, T3 for BTTB) is,

Composite =
            [ T1, T2, T3]
            [ T2, T1, T2]
            [ T3, T2, T1]     

By adding a special row and column to this composite matrix (which would amount to zero padding vectors x and b in Ax =b) we get a block circulant matrix with toeplitz blocks,

BlockCirculant=
               [ T1, T2, T3, T2]
               [ T2, T1, T2, T3]
               [ T3, T2, T1, T2]
               [ T2, T3, T2, T1]

Now the problem becomes to simply make the individual toeplitz blocks to circulant blocks, which is tackled easily by well known circular embedding. Consider T1 and its reflection coefficients matrix ref_T1,

T1 =                      ref_T1 = 
   [ a_1, a_2, a_3]               [   0, a_3, a_2]   
   [ a_2, a_1, a_2]               [ a_3,   0, a_3]  
   [ a_3, a_2, a_1]               [ a_2, a_3,   0]

Then the circulant block of doubled size in both dimensions becomes,

CirculantBlocks =                CirculantBlocks(expanded) =
        [    T1, ref_T1]                      [ a_1, a_2, a_3,   0, a_3, a_2]
        [ref_T1,     T1]                      [ a_2, a_1, a_2, a_3,   0, a_3]
                                              [ a_3, a_2, a_1, a_2, a_3,   0]
                                              [   0, a_3, a_2, a_1, a_2, a_3]
                                              [ a_3,   0, a_3, a_2, a_1, a_2] 
                                              [ a_2, a_3,   0, a_3, a_2, a_1]

Similarly we can do this for other toeplitz blocks, and finally we get a block circulant with circulant blocks matrix structure (BCCB) , which can be solved using fft2 routines as described in Computations with BCCB.

BCCB =
[     T1, ref_T1,     T2, ref_T2,     T3, ref_T3,     T2, ref_T2]
[ ref_T1,     T1, ref_T2,     T2, ref_T3,     T3, ref_T2,     T2]
[     T2, ref_T2,     T1, ref_T1,     T2, ref_T2,     T3, ref_T3]
[ ref_T2,     T2, ref_T1,     T1, ref_T2,     T2, ref_T3,     T3]
[     T3, ref_T3,     T2, ref_T2,     T1, ref_T1,     T2, ref_T2]
[ ref_T3,     T3, ref_T2,     T2, ref_T1,     T1, ref_T2,     T2]
[     T2, ref_T2,     T3, ref_T3,     T2, ref_T2,     T1, ref_T1]
[ ref_T2,     T2, ref_T3,     T3, ref_T2,     T2, ref_T1,     T1]                    

I hope, I am not abusing the answer by asking a little side question. But it develops as a concern of understanding the solution, since the size is increased in two stages, one stage where we add the special rows and columns to make entire structure go from BTTB to BCTB , then again when we add reflections and circulant embedding, to go from BCTB to BCCB. The matrices (A) I encounter grow in squared fashion with odd numbers (one dimension size) i.e. 3^2 , 5^2, 7^2,.., so consider another 25 x 25 (5^2 x 5^2) matrix with same structure,

BTTB =

[  a_1,  a_2,  a_3,  a_4,  a_5,  a_2,  a_7,  a_8,  a_9, a_10,  a_3,  a_8, a_13, a_14, a_15,  a_4,  a_9, a_14, a_19, a_20,  a_5, a_10, a_15, a_20, a_25]
[  a_2,  a_1,  a_2,  a_3,  a_4,  a_7,  a_2,  a_7,  a_8,  a_9,  a_8,  a_3,  a_8, a_13, a_14,  a_9,  a_4,  a_9, a_14, a_19, a_10,  a_5, a_10, a_15, a_20]
[  a_3,  a_2,  a_1,  a_2,  a_3,  a_8,  a_7,  a_2,  a_7,  a_8, a_13,  a_8,  a_3,  a_8, a_13, a_14,  a_9,  a_4,  a_9, a_14, a_15, a_10,  a_5, a_10, a_15]
[  a_4,  a_3,  a_2,  a_1,  a_2,  a_9,  a_8,  a_7,  a_2,  a_7, a_14, a_13,  a_8,  a_3,  a_8, a_19, a_14,  a_9,  a_4,  a_9, a_20, a_15, a_10,  a_5, a_10]
[  a_5,  a_4,  a_3,  a_2,  a_1, a_10,  a_9,  a_8,  a_7,  a_2, a_15, a_14, a_13,  a_8,  a_3, a_20, a_19, a_14,  a_9,  a_4, a_25, a_20, a_15, a_10,  a_5]
[  a_2,  a_7,  a_8,  a_9, a_10,  a_1,  a_2,  a_3,  a_4,  a_5,  a_2,  a_7,  a_8,  a_9, a_10,  a_3,  a_8, a_13, a_14, a_15,  a_4,  a_9, a_14, a_19, a_20]
[  a_7,  a_2,  a_7,  a_8,  a_9,  a_2,  a_1,  a_2,  a_3,  a_4,  a_7,  a_2,  a_7,  a_8,  a_9,  a_8,  a_3,  a_8, a_13, a_14,  a_9,  a_4,  a_9, a_14, a_19]
[  a_8,  a_7,  a_2,  a_7,  a_8,  a_3,  a_2,  a_1,  a_2,  a_3,  a_8,  a_7,  a_2,  a_7,  a_8, a_13,  a_8,  a_3,  a_8, a_13, a_14,  a_9,  a_4,  a_9, a_14]
[  a_9,  a_8,  a_7,  a_2,  a_7,  a_4,  a_3,  a_2,  a_1,  a_2,  a_9,  a_8,  a_7,  a_2,  a_7, a_14, a_13,  a_8,  a_3,  a_8, a_19, a_14,  a_9,  a_4,  a_9]
[ a_10,  a_9,  a_8,  a_7,  a_2,  a_5,  a_4,  a_3,  a_2,  a_1, a_10,  a_9,  a_8,  a_7,  a_2, a_15, a_14, a_13,  a_8,  a_3, a_20, a_19, a_14,  a_9,  a_4]
[  a_3,  a_8, a_13, a_14, a_15,  a_2,  a_7,  a_8,  a_9, a_10,  a_1,  a_2,  a_3,  a_4,  a_5,  a_2,  a_7,  a_8,  a_9, a_10,  a_3,  a_8, a_13, a_14, a_15]
[  a_8,  a_3,  a_8, a_13, a_14,  a_7,  a_2,  a_7,  a_8,  a_9,  a_2,  a_1,  a_2,  a_3,  a_4,  a_7,  a_2,  a_7,  a_8,  a_9,  a_8,  a_3,  a_8, a_13, a_14]
[ a_13,  a_8,  a_3,  a_8, a_13,  a_8,  a_7,  a_2,  a_7,  a_8,  a_3,  a_2,  a_1,  a_2,  a_3,  a_8,  a_7,  a_2,  a_7,  a_8, a_13,  a_8,  a_3,  a_8, a_13]
[ a_14, a_13,  a_8,  a_3,  a_8,  a_9,  a_8,  a_7,  a_2,  a_7,  a_4,  a_3,  a_2,  a_1,  a_2,  a_9,  a_8,  a_7,  a_2,  a_7, a_14, a_13,  a_8,  a_3,  a_8]
[ a_15, a_14, a_13,  a_8,  a_3, a_10,  a_9,  a_8,  a_7,  a_2,  a_5,  a_4,  a_3,  a_2,  a_1, a_10,  a_9,  a_8,  a_7,  a_2, a_15, a_14, a_13,  a_8,  a_3]
[  a_4,  a_9, a_14, a_19, a_20,  a_3,  a_8, a_13, a_14, a_15,  a_2,  a_7,  a_8,  a_9, a_10,  a_1,  a_2,  a_3,  a_4,  a_5,  a_2,  a_7,  a_8,  a_9, a_10]
[  a_9,  a_4,  a_9, a_14, a_19,  a_8,  a_3,  a_8, a_13, a_14,  a_7,  a_2,  a_7,  a_8,  a_9,  a_2,  a_1,  a_2,  a_3,  a_4,  a_7,  a_2,  a_7,  a_8,  a_9]
[ a_14,  a_9,  a_4,  a_9, a_14, a_13,  a_8,  a_3,  a_8, a_13,  a_8,  a_7,  a_2,  a_7,  a_8,  a_3,  a_2,  a_1,  a_2,  a_3,  a_8,  a_7,  a_2,  a_7,  a_8]
[ a_19, a_14,  a_9,  a_4,  a_9, a_14, a_13,  a_8,  a_3,  a_8,  a_9,  a_8,  a_7,  a_2,  a_7,  a_4,  a_3,  a_2,  a_1,  a_2,  a_9,  a_8,  a_7,  a_2,  a_7]
[ a_20, a_19, a_14,  a_9,  a_4, a_15, a_14, a_13,  a_8,  a_3, a_10,  a_9,  a_8,  a_7,  a_2,  a_5,  a_4,  a_3,  a_2,  a_1, a_10,  a_9,  a_8,  a_7,  a_2]
[  a_5, a_10, a_15, a_20, a_25,  a_4,  a_9, a_14, a_19, a_20,  a_3,  a_8, a_13, a_14, a_15,  a_2,  a_7,  a_8,  a_9, a_10,  a_1,  a_2,  a_3,  a_4,  a_5]
[ a_10,  a_5, a_10, a_15, a_20,  a_9,  a_4,  a_9, a_14, a_19,  a_8,  a_3,  a_8, a_13, a_14,  a_7,  a_2,  a_7,  a_8,  a_9,  a_2,  a_1,  a_2,  a_3,  a_4]
[ a_15, a_10,  a_5, a_10, a_15, a_14,  a_9,  a_4,  a_9, a_14, a_13,  a_8,  a_3,  a_8, a_13,  a_8,  a_7,  a_2,  a_7,  a_8,  a_3,  a_2,  a_1,  a_2,  a_3]
[ a_20, a_15, a_10,  a_5, a_10, a_19, a_14,  a_9,  a_4,  a_9, a_14, a_13,  a_8,  a_3,  a_8,  a_9,  a_8,  a_7,  a_2,  a_7,  a_4,  a_3,  a_2,  a_1,  a_2]
[ a_25, a_20, a_15, a_10,  a_5, a_20, a_19, a_14,  a_9,  a_4, a_15, a_14, a_13,  a_8,  a_3, a_10,  a_9,  a_8,  a_7,  a_2,  a_5,  a_4,  a_3,  a_2,  a_1]

The individual Toeplitz blocks matrices are,

T1 =
[ a_1, a_2, a_3, a_4, a_5]
[ a_2, a_1, a_2, a_3, a_4]
[ a_3, a_2, a_1, a_2, a_3]
[ a_4, a_3, a_2, a_1, a_2]
[ a_5, a_4, a_3, a_2, a_1]

T2 =
[  a_2, a_7, a_8, a_9, a_10]
[  a_7, a_2, a_7, a_8,  a_9]
[  a_8, a_7, a_2, a_7,  a_8]
[  a_9, a_8, a_7, a_2,  a_7]
[ a_10, a_9, a_8, a_7,  a_2]

T3 =
[  a_3,  a_8, a_13, a_14, a_15]
[  a_8,  a_3,  a_8, a_13, a_14]
[ a_13,  a_8,  a_3,  a_8, a_13]
[ a_14, a_13,  a_8,  a_3,  a_8]
[ a_15, a_14, a_13,  a_8,  a_3]

T4 =
[  a_4,  a_9, a_14, a_19, a_20]
[  a_9,  a_4,  a_9, a_14, a_19]
[ a_14,  a_9,  a_4,  a_9, a_14]
[ a_19, a_14,  a_9,  a_4,  a_9]
[ a_20, a_19, a_14,  a_9,  a_4]

T5 =
[  a_5, a_10, a_15, a_20, a_25]
[ a_10,  a_5, a_10, a_15, a_20]
[ a_15, a_10,  a_5, a_10, a_15]
[ a_20, a_15, a_10,  a_5, a_10]
[ a_25, a_20, a_15, a_10,  a_5]

The composite matrix which is 5 x 5 block and the circulant block which is 9 x 9 is:

Composite =
[ T1, T2, T3, T4, T5]
[ T2, T1, T2, T3, T4]
[ T3, T2, T1, T2, T3]
[ T4, T3, T2, T1, T2]
[ T5, T4, T3, T2, T1]

Block Circulant=
[ T1, T2, T3, T4, T5, T4, T3, T2]
[ T2, T1, T2, T3, T4, T5, T4, T3]
[ T3, T2, T1, T2, T3, T4, T5, T4]
[ T4, T3, T2, T1, T2, T3, T4, T5]
[ T5, T4, T3, T2, T1, T2, T3, T4]
[ T4, T5, T4, T3, T2, T1, T2, T3]
[ T3, T4, T5, T4, T3, T2, T1, T2]
[ T2, T3, T4, T5, T4, T3, T2, T1]

After circulant embedding it would become 18 x 18 (x 5 element in each) BCCB matrix, for small numbers this bloating is of no concern, however should this be of concern when I am dealing with matrix of size 1025^2 x 1025^2 ? Of-course I don't explicitly make this matrix A, whose entries I can create on the fly with a dedicated formula.