Triangular and symmetric matrices

The 15 entries on and below the diagonal of a 5 × 5 matrix, packed into one strip. Watch each row land where the rows before it end.

⋯10static int store[CELLS];11 12int index_of(int i, int j)13{14#if COLUMN15    return j * N - j * (j - 1) / 2 + (i - j);16#else17    return i * (i + 1) / 2 + j;18#endif19}20 21int get(int i, int j)22{23    return j <= i ? store[index_of(i, j)] : 0;24}25 26int sym_get(int i, int j)27{28    if (j > i) {29        int t = i;30        i = j;31        j = t;32    }33    return store[index_of(i, j)];34}35 36int main(void)37{38    for (int i = 0; i < N; i++)39        for (int j = 0; j <= i; j++)40            store[index_of(i, j)] = 10 * (i + 1) + (j + 1);41    for (int i = 0; i < N; i++) {42        for (int j = 0; j < N; j++)43            printf("%2d ", get(i, j));44        printf("\n");45    }46    printf("store:");47    for (int k = 0; k < CELLS; k++)48        printf(" %d", store[k]);49    printf("\nL(4, 2) = %d at index %d; L(2, 4) = %d; S(2, 4) = %d\n",50           get(4, 2), index_of(4, 2), get(2, 4), sym_get(2, 4));51    return 0;52}

Half a grid in one strip

A lower triangular matrix is zero above its diagonal: L(i, j) = 0 whenever j > i. Row 0 has 1 entry that can be non-zero, row 1 has 2, and row n − 1 has n. So the triangle holds

1 + 2 + … + n = n(n + 1)/2

entries. For n = 5 that is 15 of 25. An upper triangular matrix is the mirror image, zero below its diagonal, with the same count.

The program keeps the 15 entries of a 5 × 5 lower triangle in one array, static int store[15], and needs a rule that says which cell holds L(i, j). Each entry’s value is its position, 10(i + 1) + (j + 1), so you can see in the output which entry landed in which cell.

The code and its cost

Row after row

Store the rows one after another. Rows 0 to i − 1 hold 1 + 2 + … + i = i(i + 1)/2 entries, so row i starts at that cell, and L(i, j) is j cells further on:

int index_of(int i, int j)
{
    return i * (i + 1) / 2 + j;
}

The runs start at 0, 1, 3, 6 and 10, the triangular numbers. L(4, 2) is at 4 × 5 / 2 + 2 = 12. The program prints the store in order:

store: 11 21 22 31 32 33 41 42 43 44 51 52 53 54 55
L(4, 2) = 53 at index 12; L(2, 4) = 0; S(2, 4) = 53

Column after column

Store the columns instead. Column 0 holds n entries, column 1 holds n − 1, and so on, so columns 0 to j − 1 hold n + (n − 1) + … + (n − j + 1) = jn − j(j − 1)/2 entries. L(i, j) is i − j cells into column j:

return j * N - j * (j - 1) / 2 + (i - j);

Now the runs are the columns, 5, 4, 3, 2 and 1 long. The program loops over the rows exactly as before; only index_of changes, so each row’s entries land in scattered cells. Set the control to “Column-major” to watch it. The same 15 values end up in another order:

store: 11 21 31 41 51 22 32 42 52 33 43 53 44 54 55
L(4, 2) = 53 at index 11; …

The two formulas give different cells for the same entry: 12 against 11 for L(4, 2). A program must read with the formula it wrote with.

Reading an entry

get supplies the zeros above the diagonal, as the diagonal matrix’s get did:

int get(int i, int j)
{
    return j <= i ? store[index_of(i, j)] : 0;
}

The program prints all 25 entries through get, rows 11 0 0 0 0 to 51 52 53 54 55.

Upper triangles

An upper triangle is the transpose of a lower one: U(i, j) = L(j, i). Store it by swapping the roles of i and j. Row-major, row i of an upper triangle holds n − i entries, so it starts at in − i(i − 1)/2 and U(i, j) is at in − i(i − 1)/2 + (j − i).

Symmetric matrices

A symmetric matrix has M(i, j) = M(j, i). Its upper half repeats its lower half, so store only the lower triangle and swap i and j when asked for an entry above the diagonal:

int sym_get(int i, int j)
{
    if (j > i) {
        int t = i;
        i = j;
        j = t;
    }
    return store[index_of(i, j)];
}

S(2, 4) becomes S(4, 2), and both read cell 12: the clip lights them together. The program prints S(2, 4) = 53, while the triangle itself has L(2, 4) = 0.

The cost

Storage drops from n² to n(n + 1)/2, about half. Reading an entry is one formula and one array read, as fast as a full grid.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic tri.c -o tri
./tri
gcc -std=c11 -Wall -Wextra -pedantic -DCOLUMN=1 tri.c -o tri-col
./tri-col

Worked example

Two formulas, one entry. L(3, 1): row-major, 3 × 4 / 2 + 1 = 7; column-major, 1 × 5 − 0 + 2 = 7, the same by chance. L(4, 2): row-major 12, column-major 2 × 5 − 1 + 2 = 11.

A large triangle. For n = 1,000 the triangle holds 1,000 × 1,001 / 2 = 500,500 cells, against 1,000,000 for the full grid.

Distances between cities. A table of distances between 5 cities is symmetric, since the road from A to B is as long as the road from B to A, and its diagonal is 0. The entries strictly below the diagonal, 4 + 3 + 2 + 1 = 10 numbers, hold the whole table.

Two checks. Where does row 3 start in row-major order? (Cell 6.) How many cells does a symmetric 10 × 10 matrix need? (55.)

Where it’s used

Distance tables. Road atlases and route planners keep one triangle of a symmetric table.

Covariance matrices. In statistics the covariance of a set of variables is symmetric, so libraries store one triangle.

Structures. The stiffness matrix of a bridge or frame in an engineering simulation is symmetric.

Solving equations. Linear solvers factor a matrix into a lower and an upper triangle (LU decomposition); symmetric matrices have special factorisations of their own.

Reference card

Next, banded matrices keep only a few diagonals.

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.