Interval DP: matrix chain multiplication

Why this matters

Multiplying two matrices is a fixed sequence of steps, and the cost depends on the sizes of the two matrices. When you need to multiply three or more matrices in a row, the order you choose changes the total cost, sometimes by a large amount. Matrix chain multiplication is the problem of finding the cheapest order.

It is a good first interval DP problem. Interval DP solves a problem about a range of items, such as matrices 2 to 5, by trying every way to split that range into two smaller ranges. DP stands for dynamic programming, the method of solving a big problem from answers to smaller parts, and each smaller part is solved once and saved.

Some terms used in this lesson:

• A matrix is a grid of numbers with rows and columns. A matrix with 2 rows and 3 columns is written 2 x 3.
• A scalar multiplication is one multiplication of two plain numbers. The cost of a matrix multiplication is counted as the number of these.
• A chain is a list of matrices that are multiplied one after another, such as A1 times A2 times A3.
• A split point is a place inside a chain where you choose to cut it into a left part and a right part. The last multiplication of the chain happens between the two parts.
• A range is a part of the chain from matrix i to matrix j, written as i..j.

The rest of the lesson uses one example, the chain A1 x A2 x A3 with sizes 10 x 30, 30 x 5 and 5 x 60. Two orders are possible, and they cost very different amounts.

Step 1: the cost of multiplying two matrices

To multiply a matrix of size p x q by a matrix of size q x r, each of the p rows of the first matrix is combined with each of the r columns of the second matrix. Each combination takes q scalar multiplications. So the total is p x q x r scalar multiplications. The result has size p x r.

The number of columns of the first matrix must match the number of rows of the second. That is why the chain is written with sizes that share a middle number.

Now take the chain A1 (10 x 30), A2 (30 x 5) and A3 (5 x 60). Two orders are possible.

Order 1: (A1 x A2) x A3.

1. A1 times A2 is 10 x 30 times 30 x 5. The cost is 10 x 30 x 5 = 1500. The result is 10 x 5.
2. The result times A3 is 10 x 5 times 5 x 60. The cost is 10 x 5 x 60 = 3000.

Total: 1500 + 3000 = 4500.

Order 2: A1 x (A2 x A3).

1. A2 times A3 is 30 x 5 times 5 x 60. The cost is 30 x 5 x 60 = 9000. The result is 30 x 60.
2. A1 times the result is 10 x 30 times 30 x 60. The cost is 10 x 30 x 60 = 18000.

Total: 9000 + 18000 = 27000.

The first order is six times cheaper. The matrices are the same and the answer is the same, only the cost changes. The job is to find the order with the smallest total for a chain of any length.

Step 2: name the cheapest cost for each range

Write the sizes as a list of numbers called dims. A chain of n matrices needs n + 1 numbers. Matrix i has size dims[i - 1] rows by dims[i] columns, using 1-based numbering for the matrices. For the example, dims is 10, 30, 5, 60, so A1 is 10 x 30, A2 is 30 x 5 and A3 is 5 x 60.

Define cost[i][j] as the smallest number of scalar multiplications needed to multiply matrices i to j. The answer for the whole chain is cost[1][n].

Look at a range of length 1, such as cost[2][2]. It is one matrix, so no multiplication happens and the cost is 0.

Now look at a longer range, such as cost[1][3]. The last multiplication of A1 to A3 splits the range at some point k, where k is between 1 and 2. The left part is matrices 1 to k. The right part is matrices k + 1 to 3. The cost of this split has three parts:

• the cost of the left part, cost[1][k]
• the cost of the right part, cost[k + 1][3]
• the cost of the final multiplication, which joins a left result of size dims[i - 1] x dims[k] with a right result of size dims[k] x dims[j]. That final step costs dims[i - 1] x dims[k] x dims[j].

Try each split point k and keep the smallest total. This gives the rule:

cost[i][j] = the minimum over k from i to j - 1 of
             cost[i][k] + cost[k + 1][j] + dims[i - 1] x dims[k] x dims[j]

For the example, the range 1..3 has two split points.

• k = 1: cost[1][1] + cost[2][3] + 10 x 30 x 60 = 0 + 9000 + 18000 = 27000.
• k = 2: cost[1][2] + cost[3][3] + 10 x 5 x 60 = 1500 + 0 + 3000 = 4500.

The smaller value is 4500, at k = 2. So the best order splits A1 to A3 between A2 and A3, which gives (A1 x A2) x A3.

Step 3: fill the table from short ranges to long ones

The cost of a range depends on the costs of shorter ranges inside it. For example, cost[1][3] uses cost[1][1], cost[1][2], cost[2][3] and cost[3][3]. So the shorter ranges must be filled first.

The table is filled by range length. First all ranges of length 1, which cost 0. Then all ranges of length 2, then length 3, and so on up to the full chain. Within one length, the order does not matter, because ranges of the same length do not use each other.

For the example, there are three matrices:

The range 1..2 costs 10 x 30 x 5 = 1500. The range 2..3 costs 30 x 5 x 60 = 9000. Both are used when the range 1..3 is filled.

Step 4: the code

The code below uses two tables. cost holds the smallest costs. split holds, for each range, the split point k that gave the smallest cost. The split table is used afterwards to print the best order.

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <vector>
 
struct ChainResult {
    std::vector<std::vector<std::int64_t>> cost;   // cost[i][j]: cheapest way to multiply i..j
    std::vector<std::vector<std::size_t>> split;   // split[i][j]: best k, so i..k and k+1..j
};
 
// dims has n + 1 numbers. Matrix i (1-based) has size dims[i-1] rows by dims[i] columns.
ChainResult matrix_chain(const std::vector<std::int64_t>& dims) {
    std::size_t n = dims.size() - 1;
    ChainResult r{
        std::vector<std::vector<std::int64_t>>(n + 1, std::vector<std::int64_t>(n + 1, 0)),
        std::vector<std::vector<std::size_t>>(n + 1, std::vector<std::size_t>(n + 1, 0)),
    };
    for (std::size_t len = 2; len <= n; ++len) {          // range length, shortest first
        for (std::size_t i = 1; i + len - 1 <= n; ++i) {  // range start
            std::size_t j = i + len - 1;                  // range end
            std::int64_t best = std::numeric_limits<std::int64_t>::max();
            for (std::size_t k = i; k < j; ++k) {         // try every split point
                std::int64_t cost = r.cost[i][k] + r.cost[k + 1][j]
                                  + dims[i - 1] * dims[k] * dims[j];
                if (cost < best) {
                    best = cost;
                    r.split[i][j] = k;
                }
            }
            r.cost[i][j] = best;
        }
    }
    return r;
}
 
std::string parenthesise(const ChainResult& r, std::size_t i, std::size_t j) {
    if (i == j) {
        return "A" + std::to_string(i);
    }
    std::size_t k = r.split[i][j];
    return "(" + parenthesise(r, i, k) + " " + parenthesise(r, k + 1, j) + ")";
}
 
// Cost of multiplying left to right, for comparison.
std::int64_t left_to_right(const std::vector<std::int64_t>& dims) {
    std::size_t n = dims.size() - 1;
    std::int64_t total = 0;
    std::int64_t rows = dims[0];
    for (std::size_t k = 2; k <= n; ++k) {
        total += rows * dims[k - 1] * dims[k];
    }
    return total;
}
 
void show(const std::vector<std::int64_t>& dims) {
    ChainResult r = matrix_chain(dims);
    std::size_t n = dims.size() - 1;
    std::cout << "dims:";
    for (auto d : dims) {
        std::cout << ' ' << d;
    }
    std::cout << "\nminimum cost: " << r.cost[1][n] << '\n';
    std::cout << "best order:   " << parenthesise(r, 1, n) << '\n';
    std::cout << "left to right: " << left_to_right(dims) << "\n\n";
}
 
int main() {
    show({10, 30, 5, 60});
    show({30, 35, 15, 5, 10, 20, 25});
 
    std::cout << "table for dims 10 30 5 60 (row i, column j):\n";
    ChainResult r = matrix_chain({10, 30, 5, 60});
    for (std::size_t i = 1; i <= 3; ++i) {
        for (std::size_t j = 1; j <= 3; ++j) {
            std::cout << std::setw(8) << (j >= i ? r.cost[i][j] : 0);
        }
        std::cout << '\n';
    }
}

Save the code as chain.cpp. Compile and run it with g++ -std=c++20 -Wall -o chain chain.cpp and then ./chain. The output is:

dims: 10 30 5 60
minimum cost: 4500
best order:   ((A1 A2) A3)
left to right: 4500
 
dims: 30 35 15 5 10 20 25
minimum cost: 15125
best order:   ((A1 (A2 A3)) ((A4 A5) A6))
left to right: 40500
 
table for dims 10 30 5 60 (row i, column j):
       0    1500    4500
       0       0    9000
       0       0       0

Some details of the code:

• matrix_chain uses 1-based matrix numbers, so the tables have one extra row and column that are not used. This keeps the formulas close to the math.
• The cost < best check keeps the first split that reaches the smallest value. Ties go to the smallest k.
• parenthesise rebuilds the order by reading split. It is a recursive function, meaning it calls itself on the smaller left and right parts until it reaches a single matrix.
• left_to_right computes the cost of the order that multiplies from the left, A1 times A2, then times A3, and so on. For the first example, that order happens to be the best one. For the second it costs 40500, against the best of 15125.

Example 2: six matrices by hand

The second example in the program uses six matrices, with dims 30, 35, 15, 5, 10, 20, 25. The matrices are A1 (30 x 35), A2 (35 x 15), A3 (15 x 5), A4 (5 x 10), A5 (10 x 20) and A6 (20 x 25). The program prints 15125 with the order ((A1 (A2 A3)) ((A4 A5) A6)). Here is the same total worked out by hand, using the order it found.

The top-level split is after A3. The left part is A1 to A3, and the right part is A4 to A6.

1. The left part is A1 x (A2 x A3).
   • A2 x A3 is 35 x 15 times 15 x 5. Cost: 35 x 15 x 5 = 2625. Result: 35 x 5.
   • A1 x (35 x 5) is 30 x 35 times 35 x 5. Cost: 30 x 35 x 5 = 5250. Result: 30 x 5.
   • Left part total: 2625 + 5250 = 7875.
2. The right part is (A4 x A5) x A6.
   • A4 x A5 is 5 x 10 times 10 x 20. Cost: 5 x 10 x 20 = 1000. Result: 5 x 20.
   • (5 x 20) times A6 is 5 x 20 times 20 x 25. Cost: 5 x 20 x 25 = 2500. Result: 5 x 25.
   • Right part total: 1000 + 2500 = 3500.
3. The final multiplication joins a 30 x 5 result with a 5 x 25 result. Cost: 30 x 5 x 25 = 3750.

The total is 7875 + 3500 + 3750 = 15125, which matches the program. In the formula, the final step is dims[0] x dims[3] x dims[6], which is 30 x 5 x 25.

Example 3: a chain of one matrix

If the chain has a single matrix, no multiplication happens, so the cost is 0. The table has only cost[1][1], and the loop over lengths starts at 2, so it never runs. The result is 0.

The input dims must still have at least two numbers. A chain with one matrix needs two numbers, for example {7, 3}, which is one 7 x 3 matrix. The input {7} has one number, which describes no matrix at all. The program does not check for this, and it fails. The mistakes section shows the message.

Common mistakes

Using the wrong size for the left end

Matrix i has size dims[i - 1] by dims[i]. The left end of a range has the size dims[i - 1], not dims[i]. Writing dims[i] gives a wrong total without any error message.

This is a version of matrix_chain with dims[i] in place of dims[i - 1], all other lines the same. The output for the first example is:

wrong index: 13500 (correct is 4500)

The wrong total is still a plausible number, so a test with a known answer is the best way to catch this. The example 10, 30, 5, 60 has a known answer of 4500.

Filling the table in the wrong order

If the outer loop goes over the starting position, and the inner loop goes over the ending position, then a longer range can be computed before one of its shorter parts. Its part is still 0 in the table, so the cost is too small. The loop over the range length in the code above avoids this.

Reading past the end of the table

The program takes dims with n + 1 numbers, and it reads cost[1][n]. If dims has only one number, n is 0 and the table has one row. Reading cost[1] is then past the end.

The g++ 15 build used for this lesson checks this index by default and stops with an assertion. This is the message from the run where show({7}) was called:

/usr/include/c++/15/bits/stl_vector.h:1263: constexpr std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = std::vector<long int>; _Alloc = std::allocator<std::vector<long int> >; reference = std::vector<long int>&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.

The program stops at this line and prints nothing else, because the text printed before the crash was still in the output buffer. The message shows the file, the line in the library, and the index that was out of range. A check that dims.size() >= 2 at the start of matrix_chain would give a clearer error.

Choosing the first multiplication by a rule of thumb

A rule such as "multiply the cheapest adjacent pair first" is a guess. It can give the right answer for one set of sizes and a wrong one for another. The table in this lesson tries every split point, so it does not depend on a guess.

Practice

Exercise 1

Find the cost of the best order for dims 1, 2, 3, 4, by hand. Then say which order gives it.

Show the answer

There are three matrices: A1 is 1 x 2, A2 is 2 x 3 and A3 is 3 x 4.

• Order (A1 x A2) x A3: 1 x 2 x 3 = 6, then 1 x 3 x 4 = 12. Total 18.
• Order A1 x (A2 x A3): 2 x 3 x 4 = 24, then 1 x 2 x 4 = 8. Total 32.

The best cost is 18, from the order ((A1 A2) A3).

Exercise 2

The program in this lesson gives 26000 for dims 40, 20, 30, 10, 30. Change show in a copy of the code to check this, and say which split point is chosen at the top level.

Show the answer

Running show({40, 20, 30, 10, 30}) prints:

dims: 40 20 30 10 30
minimum cost: 26000
best order:   ((A1 (A2 A3)) A4)

So the top-level split is after A1. Checking the order by hand: A2 x A3 costs 20 x 30 x 10 = 6000, giving a 20 x 10 result. A1 times that costs 40 x 20 x 10 = 8000, giving a 40 x 10 result. Then A4 costs 40 x 10 x 30 = 12000. The total is 6000 + 8000 + 12000 = 26000.

Exercise 3

Why must the code fill the table for ranges of length 2 before ranges of length 3?

Show the answer

The cost of a range of length 3 is found by adding the costs of two shorter ranges and a final multiplication. Those shorter ranges must already hold their final values. If the length-3 range were filled first, it would read zeros for its shorter parts and report a total that is too small.

Exercise 4

Write the rule for cost[i][j] in your own words, using the names i, j, k and dims.

Show the answer

The cheapest way to multiply matrices i to j is the smallest value of this sum, taken over every k from i to j - 1. The sum is the cost of the cheapest way to multiply matrices i to k, plus the cost of the cheapest way to multiply matrices k + 1 to j, plus the cost of joining the two results. Joining costs dims[i - 1] times dims[k] times dims[j].

Summary

• The order of a matrix chain changes the total cost of the scalar multiplications, though not the result.
• Multiplying a p x q matrix by a q x r matrix costs p x q x r scalar multiplications.
• For each range i..j, try every split point k and keep the smallest total of the two parts and the final join.
• Fill the table by range length, from length 2 up, so each shorter range is final before it is used.
• Store the best split point for each range, and use it to print the best order.