Задача о порядке перемножения матриц

Задача о порядке перемножения матриц — классическая задача динамического программирования, в которой дана последовательность матриц $A_1, A_2, ..., A_n$ и требуется минимизировать количество скалярных операций для вычисления их произведения. Матрицы предполагаются совместимыми по отношению к матричному умножению (то есть количество столбцов $A_{i - 1}$ совпадает с количеством строк $A_i$ матрицы).

Подробное описание задачи

Произведение матриц — ассоциативная операция, которая принимает на вход две матрицы размером k×m и m×n и возвращает матрицу размером k×n, потратив на это kmn операций умножения^[1].

Когда матрицы велики по одному измерению и малы по другому, количество скалярных операций может серьёзно зависеть от порядка перемножений матриц. Допустим, нам даны 3 матрицы $A_1, A_2, A_3$ размерами соответственно 10×100, 100×5 и 5×50. Существует 2 способа их перемножения (расстановки скобок): $((A_1A_2)A_3)$ и $(A_1(A_2A_3))$. В первом случае нам потребуется 10·100·5 + 10·5·50 = 7500 скалярных умножений, а во втором случае 100·5·50 + 10·100·50 = 75000 умножений — разница налицо. Поэтому может быть выгоднее потратить некоторое время на предобработку, решив, в каком порядке лучше всего умножать, чем умножать сразу в лоб.

Таким образом, даны n матриц: $A_1: \, p_0 \times p_1$, $A_2: \, p_1 \times p_2$, …, $A_n: \, p_{n-1} \times p_{n}$. Требуется определить, в каком порядке перемножать их, чтобы количество операций умножения было минимальным.

Решение задачи

Разберём 2 способа решения задачи, чтобы показать, насколько выгодно динамическое программирование в данном случае.

Перебор всех вариантов расстановки скобок

Оценим, сколько же нужно перебрать вариантов расстановки. Обозначим через P(n) количество способов расстановки скобок в последовательности, состоящей из n матриц. Когда матрица одна, то расставлять нечего, вариант один. Если n >= 2, то количество вариантов, которым можно расставить скобки является произведением количества вариантов, которым можно расставить скобки в составляющих результирующую матрицу произведениях.(e.g. Если $A_3 = A_1*A_2$, то количество вариантов, которым мы можем получить матрицу A3 равно произведению количества способов получить матрицу $A_1$ на количество способов получить $A_2$). Разбиение на матрицы, $A_1$ и $A_2$ может производиться на границе k-й и (k + 1)-й матриц для k = 1, 2, …, n — 1. Получаем рекуррентное соотношение: $P(n) = \left \{ \begin{array}{ll} 1,\ n=1 \\ \sum_{k=1}^{n-1} P(k)P(n-k),\ n>=2 \end{array} \right.$

Решением аналогичного рекуррентного соотношения является последовательность чисел Каталана, возрастающая как $\Omega$(${4^n}\over {n^{1.5}}$). Зависимость получается экспоненциальная, непригодная для практического применения в программах. Разберём более перспективный способ.

Динамическое программирование

1.Сведение задачи к подзадачам

Обозначим результат перемножения матриц $A_iA_{(i+1)}... A_j$ через $A_{i..j}$, где i<=j. Если i<j, то существует такое k, которое разбивает $A_{i..j}$ между матрицами $A_k$ и $A_{k+1}$, i<=k<j. То есть для вычисления $A_{i..j}$ надо сначала вычислить $A_{i..k}$, потом $A_{k+1..j}$ и затем их перемножить. Выбор k является аналогом расстановки скобок между матрицами. Выбрав некоторое k мы свели задачу к двум аналогичным подзадачам для матриц $A_{i..k}$ и $A_{k+1..j}$.

2.Рекурсивное решение

Обозначим через m[i, j] минимальное количество скалярных умножений для вычисления матрицы $A_{i..j}$. Получаем следующее рекуррентное соотношение: $m[i,j] = \left \{ \begin{array}{ll} 0, & i=j \\ \min\limits_{i \le k < j}\{m[i,k] + m[k+1,j] + p_{i-1}p_kp_j \}, & i < j \end{array} \right.$

Объясняется оно просто: для того, чтобы найти произведение матриц $A_{i..j}$ при i=j не нужно ничего делать — это и есть сама матрица $A_i$. При нетривиальном случае мы перебираем все точки разбиения матрицы $A_{i..j}$ на матрицы $A_{i..k}$ и $A_{k+1..j}$, ищем кол-во операций, необходимое чтобы их получить и затем перемножаем для получения матрицы $A_{i..j}$.(Оно будет равно кол-ву операций, потраченное на решение подзадач + стоимость умножения матриц $A_{i..k}A_{k+1..j}$). Считаем, что размеры матриц заданы в массиве $p$ и размер матрицы $A_i$ равен $p_{i-1} \times p_i$. Как обычно рекурсивный метод нельзя использовать напрямую — он будет экспоненциальным из-за большого кол-ва перекрывающихся подзадач.

3.Динамическое программирование

Будем запоминать в двумерном массиве m результаты вычислений для подзадач, чтобы избежать пересчета для уже вычислявшихся подзадач. После вычислений ответ будет в m[1,n](Сколько перемножений требуется для последовательности матриц от 1 до n — то есть ответ на поставленную задачу).Сложность алгоритма будет O$(n^3)$, так как у нас ${n \choose 2}$ вариантов выбора i, j : 1 <= i <= j <= n и O(N) точек разделения для каждого варианта. Детали станут понятны из реализации.

Реализация

Java

В методе main вбит пример из начала статьи, если запустить выведет 7500, как и ожидается.

public class MatrixChain {
       /**
         * Возвращает ответ на задачу об оптимальном перемножении матриц, используя
         * динамическое программирование.
         * Асимптотика решения - O(N^3) время и O(N^2) память.
         *  
         * @param p  массив размеров матриц, см.статью   
         * @return   минимальное количество скалярных умножений, необходимое для решения задачи
         */     
       public static int multiplyOrder(int[] p) {
                int n = p.length - 1;
                int[][] dp = new int[n + 1][n + 1];
 
                for (int i = 1; i <= n; i++) {
                        dp[i][i] = 0;
                }
 
                for (int l = 2; l <= n; l++) {
                        for (int i = 1; i <= n - l + 1; i++) {
                                int j = i + l - 1;
                                dp[i][j] = Integer.MAX_VALUE;
                                for (int k = i; k <= j - 1; k++) {
                                        dp[i][j] = Math.min(dp[i][j],
                                                   dp[i][k] + dp[k + 1][j] + p[i - 1] * p[k] * p[j]);
                                }
                        }
                }
                return dp[1][n]; 
        }
 
        public static void main(String[] args) {
                int[] test = { 10, 100, 5, 50 };
                System.out.println(MatrixChain.multiplyOrder(test));
        }
}

C#

Приведены только методы, которые непосредственно выполняют необходимые расчеты

dataStore — объект класса, который хранит все данные

Его атрибуты:

public List<List<int>> m;               //матрица m
public List<List<int>> s;                               //матрица s
public List<List<int>> result;                  //результат всех перемножений
public List<List<List<int>>> source;    //массив из 2-мерных матриц (A0,A1,...,An) которые нужно перемножить
public List<int> sizes = new List<int>();       //размеры матриц (записаны таким образом - 12,10,7,4 => значит 3 матрицы размерами 12x10,10x7,7x4)
public string order = new string('a', 0);       //правильное расположение скобок

Функциональные участки кода:

//© Paulskit 27.03.2010
//метод который находит матрицу m и s (там же под них и выделяется память)
private void matrixChainOrder(){
        int n = dataStore.sizes.Count - 1;
 
        //выделяем память под матрицы m и s
        dataStore.m = new List<List<int>>();
        dataStore.s = new List<List<int>>();
        for (int i = 0; i < n; i++){
            dataStore.m.Add(new List<int>());
            dataStore.s.Add(new List<int>());
            //заполняем нулевыми элементами
            for (int a = 0; a < n; a++) {
                dataStore.m[i].Add(0);
                dataStore.s[i].Add(0);
            }
        }
        //выполняем итерационный алгоритм
        int j;
        for (int l = 1; l < n; l++)
            for (int i = 0; i < n - l; i++) {
                j = i + l;
                dataStore.m[i][j] = int.MaxValue;
                for (int k = i; k < j; k++) {
 
                    int q = dataStore.m[i][k] + dataStore.m[k + 1][j] +
                        dataStore.sizes[i] * dataStore.sizes[k + 1] * dataStore.sizes[j + 1];
                    if (q < dataStore.m[i][j]) {
                        dataStore.m[i][j] = q;
                        dataStore.s[i][j] = k;
                    }
                }
            }
}
 
//метод - простое перемножение 2-х матриц
private List<List<int>> matrixMultiply(List<List<int>> A, List<List<int>> B) {
    int rowsA = A.Count;
    int columnsB = B[0].Count;
    //column count of A == rows count of B
    int columnsA = B.Count;
 
    //memory alloc for "c"
    List<List<int>> c = new List<List<int>>();
    for (int i = 0; i < rowsA; i++) {
        c.Add(new List<int>());
        for (int a = 0; a < columnsB; a++) {
            c[i].Add(0);
        }
    }
 
    //do multiplying
    for (int i = 0; i < rowsA; i++)
        for (int j = 0; j < columnsB; j++)
            for (int k = 0; k < columnsA; k++) 
                c[i][j] += A[i][k] * B[k][j];
 
    //return value
    return c;
}
 
//метод, который непосредственно выполняет перемножение в правильном порядке
//первоначально вызывается таким образом
//dataStore.result = matrixChainMultiply(0, dataStore.sizes.Count - 2); 
private List<List<int>> matrixChainMultiply(int i, int j) {
        if (j > i) {
        List<List<int>> x = matrixChainMultiply(i, dataStore.s[i][j]);
        List<List<int>> y = matrixChainMultiply(dataStore.s[i][j] + 1, j);
        return matrixMultiply(x, y);
    }
    else return dataStore.source[i];
}
 
//метод печатающий строку с правильной расстановкой скобок
 
private void printOrder(int i, int j){
    if(i==j) order += "A"+i.ToString();
    else {
        order+="(";
        printOrder(i,dataStore.s[i][j]);
        order+="*";
        printOrder(dataStore.s[i][j]+1,j);
        order+=")";
    }
}

Примечания

К данной задаче сводится задача оптимизации свободной энергии молекулы РНК в биоинформатике (здесь пара скобок в строке мономеров РНК определяет их спаривание). Подобное динамическое программирование реализовано в алгоритмах Nussinov или Zucker.

1. ↑ Существуют и более быстрые, чем kmn, алгоритмы умножения заполненных матриц, но они применяются крайне редко — прирост скорости наблюдается только на матрицах 100×100 и крупнее. Разреженные же матрицы умножают особыми алгоритмами, оптимизированными под ту или иную форму матрицы.

Литература

• Томас Х. Кормен и др. Алгоритмы: построение и анализ = INTRODUCTION TO ALGORITHMS. — 2-е изд. — М.: «Вильямс», 2006. — С. 1296. — ISBN 0-07-013151-1

• Sanjoy Dasgupta , Christos H. Papadimitriou, Umesh Vazirani Algorithms = Algorithms. — 1-е изд. — McGraw-Hill Science/Engineering/Math, 2006. — С. 336. — ISBN 0073523402

Для улучшения этой статьи желательно?: Исправить статью согласно стилистическим правилам Википедии. Викифицировать статью.