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.

Cross-posting (asking the same question on multiple sites) is frowned upon, because it can lead to duplication of effort and hence a waste of people's time. You should wait for an answer at M.SE before giving up and trying here. — user43208, Jul 05 '15 at 07:09
This is my very first post here. I will take care of this in future. But I didn't knew which forum was better for my post, both sounded alike , nor did I readily find the description of scope for both forums that I could compare promptly and decide. And finally in my defense, I did mention in the second post that I have already posted this else while expressing my hesitation to cross-post. However thanks for the comment. — Syed Alam Abbas, Jul 06 '15 at 01:00

score 1 · Accepted Answer · answered Jul 05 '15 at 10:50

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.

answered Jul 05 '15 at 10:50

Federico Poloni

3,465

Excellent observation. I would have really struggled to find this structure, I wasn't aware of it until now, that such a structure exists and has been formalized, studied and used (like image deblurring etc.). Thanks very much, this was a precise response. With your help I can now answer my question with more description. And your guess is right, I want this solution to be asymptotically fast, since I am dealing large matrices of similar structure. That's why I am looking for a fast solution. – Syed Alam Abbas Jul 06 '15 at 01:28

score 1 · Answer 2 · edited Apr 13 '17 at 12:20

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.

The extent of bloating to transform this symmetric positive definite block toeplitz matrix A of size n x n , where n is square of an odd number, is roughly 4n x 4n which includes operations to make block toeplitz to block circulant, and circulant embedding of individual blocks. Using BCCB , the solution of inversion can be obtained in order of O(n log(n)), rather than eigenvalues computed directly, O(n^3), when n is sufficiently large, a factor of 4 shouldn't matter. I think my assessment is correct, comments are welcome. — Syed Alam Abbas, Jul 06 '15 at 04:42

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

2 Answers2

Linked