Toeplitz matrices

Seven numbers, one per diagonal, are enough for a 4 × 4 matrix. Watch each entry read from its diagonal’s number, then a filter built the same way.

⋯7static const int t[2 * N - 1] = { 7, 6, 5, 1, 2, 3, 4 };8static const int h[TAPS] = { 1, 2, 1 };9 10int get(int i, int j)11{12    return t[N - 1 + j - i];13}14 15int conv_entry(int i, int j)16{17    int k = i - j;18    return (k >= 0 && k < TAPS) ? h[k] : 0;19}20 21int main(void)22{23    for (int i = 0; i < N; i++) {24        for (int j = 0; j < N; j++)25            printf("%d ", get(i, j));26        printf("\n");27    }28    const int x[N] = { 1, 2, 3, 4 };29    printf("Mx =");30    for (int i = 0; i < N; i++) {31        int s = 0;32        for (int j = 0; j < N; j++)33            s += get(i, j) * x[j];34        printf(" %d", s);35    }36    const int u[LEN] = { 1, 2, 3, 4, 5 };37    printf("\nfiltered =");38    for (int i = 0; i < LEN + TAPS - 1; i++) {39        int s = 0;40        for (int j = 0; j < LEN; j++)41            s += conv_entry(i, j) * u[j];42        printf(" %d", s);43    }44    printf("\n");45    return 0;46}

One number per diagonal

A Toeplitz matrix is constant along every diagonal: M(i, j) depends only on j − i. In an n × n matrix j − i runs from −(n − 1) to n − 1, so the whole matrix is fixed by 2n − 1 numbers: its first column and its first row, which share the corner. The program’s 4 × 4 example is

1 2 3 4
5 1 2 3
6 5 1 2
7 6 5 1

and its seven numbers are stored in one array, the first column read upwards and then the first row:

static const int t[2 * N - 1] = { 7, 6, 5, 1, 2, 3, 4 };

A Toeplitz matrix need not be symmetric. This one is not: M(1, 0) = 5 but M(0, 1) = 2. It is symmetric only when the first row and the first column hold the same numbers.

The code and its cost

The index

Diagonal j − i is stored at t[n − 1 + j − i]. The offset n − 1 moves the lowest diagonal, j − i = −(n − 1), to cell 0:

int get(int i, int j)
{
    return t[N - 1 + j - i];
}

M(0, 0) is t[3] = 1, M(3, 0) is t[0] = 7, and M(0, 3) is t[6] = 4. The program prints all 16 entries through get; each row is the row above shifted one place right, with a new number coming in on the left.

The product

Mx needs every entry once, so it takes n² multiplications, as a full matrix does; the entries just come from 7 numbers instead of 16:

for (int i = 0; i < N; i++) {
    int s = 0;
    for (int j = 0; j < N; j++)
        s += get(i, j) * x[j];
    printf(" %d", s);
}

With x = (1, 2, 3, 4) the program prints Mx = 30 25 27 38. The saving is the memory: 2n − 1 numbers instead of n².

Convolution as a matrix

Convolution with a short filter is a Toeplitz matrix. Filter u = (1, 2, 3, 4, 5) with the taps h = (1, 2, 1): output i is the sum of h[k] × u[i − k]. Written as a matrix, entry (i, j) is h[i − j] when 0 ≤ i − j < 3, and 0 otherwise. It depends only on i − j, so the matrix is Toeplitz, and it is banded: only three diagonals are non-zero. A k-tap filter on L samples gives L + k − 1 outputs, so the matrix is (L + k − 1) × L, here 7 × 5:

int conv_entry(int i, int j)
{
    int k = i - j;
    return (k >= 0 && k < TAPS) ? h[k] : 0;
}

The program never stores the 35 entries; it works each one out from h when it needs it. Each column of the matrix holds the taps 1, 2, 1, one row lower than the column before. Multiplying by u gives

filtered = 1 4 8 12 16 14 5

which adds to 60, the input’s sum, 15, times the taps’ sum, 4.

The cost

A Toeplitz matrix takes 2n − 1 numbers instead of n². Reading an entry is one subtraction and one array read. The product with a vector still takes n² multiplications when done this way.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic toeplitz.c -o toeplitz
./toeplitz

Worked example

Row 2 of Mx. Row 2 of M is 6 5 1 2, so its entry of Mx is 6 × 1 + 5 × 2 + 1 × 3 + 2 × 4 = 6 + 10 + 3 + 8 = 27.

One filtered sample. Row 3 of the filter matrix has the taps in columns 1, 2 and 3: filtered[3] = 1 × 2 + 2 × 3 + 1 × 4 = 12.

A large matrix. For n = 1,000 the store holds 2 × 1,000 − 1 = 1,999 numbers against 1,000,000 entries.

Two checks. Which cell of t holds M(3, 1)? (3 + 1 − 3 = 1, the value 6.) How large is the matrix for a 4-tap filter on 10 samples? (13 × 10.)

Where it’s used

Digital filters. A filter applied to a block of samples is a banded Toeplitz matrix times the block; writing it that way lets a filter be analysed with the tools of linear algebra.

Speech coding. Linear prediction solves a system whose matrix, the autocorrelation matrix, is symmetric Toeplitz; fast solvers use the structure to work in about n² steps instead of n³.

Image blurring. Blurring an image with a small kernel is a block Toeplitz matrix: Toeplitz blocks arranged in a Toeplitz pattern.

Time series. Models of a series that depend only on the lag between samples give Toeplitz covariance matrices.

Reference card

Next, sparse matrices store only the entries that are not zero, wherever they are. Multiplying polynomials is the same convolution.

Platform: C11, gcc on x86-64 Linux: int 4 bytes, pointers 8, little-endian. Addresses in the pictures are short stand-ins, by one rule for the whole course; sizes are the real ones.