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

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

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

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

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

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

Решение задачиПравить

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

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

Оценим, сколько же нужно перебрать вариантов расстановки. Обозначим через $P(n)$$P(n)$  количество способов расстановки скобок в последовательности, состоящей из n матриц. Когда матрица одна, то расставлять нечего, вариант один. Если $n\ge <!--\; \ge \; -->2$$n\geq 2$ , то количество вариантов, которым можно расставить скобки является произведением количества вариантов, которым можно расставить скобки в составляющих результирующую матрицу произведениях (т.е. если $A_3=A_1\ast <!--\; \ast \; -->A_2$$A_{3}=A_{1}*A_{2}$ , то количество вариантов, которым мы можем получить матрицу $A_3$$A_{3}$  равно произведению количества способов получить матрицу $A_1$$A_{1}$  на количество способов получить $A_2$$A_{2}$ ). Разбиение на матрицы, $A_1$$A_{1}$  и $A_2$$A_{2}$  может производиться на границе k-й и (k + 1)-й матриц для $k=1,2,\dots <!--\; \dots \; -->,n-<!--\; -\; -->1$$k=1,2,\dots ,n-1$ . Получаем рекуррентное соотношение: $P(n)=\{\begin{array}{l}1,n=1\\ \underset{k=1}{\overset{n-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}P(k)P(n-<!--\; -\; -->k),n>=2\end{array}$$P(n)=\left\{{\begin{array}{ll}1,\ n=1\\\sum _{k=1}^{n-1}P(k)P(n-k),\ n>=2\end{array}}\right.$ 

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

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

Сведение задачи к подзадачамПравить

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

Рекурсивное решениеПравить

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

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

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

Будем запоминать в двумерном массиве m результаты вычислений для подзадач, чтобы избежать пересчета для уже вычислявшихся подзадач. После вычислений ответ будет в m[1,n](Сколько перемножений требуется для последовательности матриц от 1 до n — то есть ответ на поставленную задачу).Сложность алгоритма будет O$\left(n^3\right)$$\left(n^{3}\right)$ , так как у нас $(\genfrac{}{}{0ex}{}{n}{2})$${n \choose 2}$  вариантов выбора i, j : $1\le <!--\; \le \; -->i\le <!--\; \le \; -->j\le <!--\; \le \; -->n$$1\leq i\leq j\leq n$  и $O(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][n];
		
		for (int i = 0; i < n; ++i) {
			dp[i][i] = 0;
		}
		
		for (int l = 1; l < n; ++l) {
			for (int i = 0; i < n - l; ++i) {
				int j = i + l;
				dp[i][j] = Integer.MAX_VALUE;
				for (int k = i; k < j; ++k) {
					dp[i][j] = Math.min(dp[i][j],
                                        dp[i][k] + dp[k + 1][j] + p[i] * p[k + 1] * p[j + 1]);
				}
			}
		}
		return dp[0][n - 1]; 
	}
	
	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.

ЛитератураПравить

• Томас Х. Кормен и др. Алгоритмы: построение и анализ = 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.