Diagonal matrices

A 4 × 4 matrix with numbers only on its diagonal. Watch the diagonal drop into a strip of 4 cells while the zeros stay behind.

⋯10    int n;11    double *d;12};13 14bool diag_init(struct diag *m, int n)15{16    m->n = n;17    m->d = calloc((size_t)n, sizeof *m->d);18    return m->d != NULL;19}20 21double get(const struct diag *m, int i, int j)22{23    return i == j ? m->d[i] : 0.0;24}25 26bool set(struct diag *m, int i, int j, double x)27{28    if (i != j)29        return x == 0.0;30    m->d[i] = x;31    return true;32}33 34int main(void)35{36    struct diag m;37    if (!diag_init(&m, N))38        return 1;39    for (int i = 0; i < N; i++)40        set(&m, i, i, i + 2);41    for (int i = 0; i < N; i++) {42        for (int j = 0; j < N; j++)43            printf("%g ", get(&m, i, j));44        printf("\n");45    }46    printf("stored %d of %d cells\n", N, N * N);47    bool ok = set(&m, 0, 1, 9);48    printf("set(0, 1, 9): %s\n", ok ? "done" : "refused");49    double x[N], y[N];50    for (int i = 0; i < N; i++)51        x[i] = i + 1;52    for (int i = 0; i < N; i++)53        y[i] = m.d[i] * x[i];54    printf("y =");55    for (int i = 0; i < N; i++)56        printf(" %g", y[i]);57    printf("\n");58    free(m.d);59    return 0;60}

Only n numbers

A diagonal matrix is zero everywhere off its diagonal: M(i, j) = 0 whenever i ≠ j. An n × n one has n² entries, but only n of them can be anything other than zero. So the program stores those n and nothing else:

struct diag {
    int n;
    double *d;
};

d points at a block of n doubles, one per diagonal entry: d[i] is M(i, i). The program asks calloc for it, so every cell starts at 0:

m->d = calloc((size_t)n, sizeof *m->d);

For n = 4 that is 4 cells of 8 bytes, 32 bytes, where the full grid would need 16 cells and 128 bytes. The other 12 entries are not stored anywhere.

The code and its cost

Reading an entry

The zeros come from the code, not from memory. get checks where the entry is:

double get(const struct diag *m, int i, int j)
{
    return i == j ? m->d[i] : 0.0;
}

On the diagonal it reads d[i]; anywhere else it returns 0 without reading memory at all. The program prints the whole matrix through get, so all 16 entries appear on screen, zeros included:

2 0 0 0
0 3 0 0
0 0 4 0
0 0 0 5
stored 4 of 16 cells

Writing an entry

set writes the diagonal and nothing else:

bool set(struct diag *m, int i, int j, double x)
{
    if (i != j)
        return x == 0.0;
    m->d[i] = x;
    return true;
}

Off the diagonal there is no cell to write. A zero is accepted, because the matrix already holds 0 there and nothing changes. Any other value is refused, and set returns false. That way the store can never disagree with the matrix it claims to be. In the clip, set(0, 1, 9) is refused and the program prints set(0, 1, 9): refused.

Multiplying by a vector

Row i of M has a single non-zero entry, d[i], in column i. So entry i of y = Mx is just d[i] × x[i]:

for (int i = 0; i < N; i++)
    y[i] = m.d[i] * x[i];

That is n multiplications, where a full matrix would take n². With x = (1, 2, 3, 4) the program prints y = 2 6 12 20. A diagonal matrix scales each coordinate by its own factor, which is what diagonalisation is looking for.

The cost

Storing a diagonal matrix takes n cells instead of n². Reading an entry and writing one take one step each. Multiplying by a vector takes n steps instead of n².

Running it yourself

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

The size is a macro, so gcc -DN=8 … builds the 8 × 8 version, which prints stored 8 of 64 cells and y = 2 6 12 20 30 42 56 72.

Worked example

A large matrix. For n = 1,000, the strip holds 1,000 doubles, 8,000 bytes. The full grid would hold 1,000,000 doubles, 8,000,000 bytes, of which all but 1,000 are zeros.

Two diagonal matrices multiplied. If D has diagonal d and E has diagonal e, then DE is diagonal too, with entries d[i] × e[i]. That is n multiplications, and the product needs only n cells.

Two reads. For n = 4, get(2, 2) returns d[2] = 4, and get(2, 3) returns 0 without reading memory.

Two checks. Why does set(&m, 0, 1, 0) succeed? (Entry (0, 1) is already 0, so nothing changes.) How many cells does an 8 × 8 diagonal matrix store? (8, against 64.)

Where it’s used

Gains on channels. Scaling each channel of a signal by its own gain is multiplying by a diagonal matrix, one multiplication per channel.

Weighted fits. A weighted least-squares fit gives each measurement a weight; the weights sit on the diagonal of a matrix that is never stored in full.

Matrix libraries. Numerical libraries treat diagonal matrices as a special case, storing only the diagonal and skipping the zeros in every product.

Reference card

Next, triangular and symmetric matrices store half a grid with a formula for each entry’s place.

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