скачать doc
Динамическое программирование
Задача о порядке перемножения матриц
Рассмотрим произведение матриц M1 M2 M3 ... Mn 1 Mn. Каждая матрица Mi имеет размер ri 1 ri. Вычисление произведения двух матриц – размер первой l m и размер второй m n – требует l m n умножений их элементов. Общее количество элементарных операций умножения, требуемое при вычислении произведения цепочки матриц, зависит от порядка, в котором производятся попарные умножения матриц. Требуется найти такой порядок перемножения матриц, который минимизирует общее количество элементарных операций умножения.
Рассмотрим пример M1 M2 M3 M4, где M1 имеет размер 1020, M2 2050, M3 501, а M4 1100. Если умножения происходят в порядке, соответствующем расстановке скобок M1 (M2 (M3 M4)), то потребуется 125000 умножений, а если это же произведение вычислять как (M1 (M2 M3)) M4, то потребуется всего 2200 умножений.
Пусть mij – оптимальное количество умножений, требуемое для вычисления произведения цепочки матриц M(i, j) = Mi Mi +1 ... Mj 1 Mj, где 1 i j n. Очевидно, что M(i, i) = Mi и mii = 0, а m1n – соответствует решению задачи для исходной цепочки M(1, n) . При 1 i j n справедливо следующее рекуррентное соотношение
mij = Min { mik + mk +1,j + ri 1 rk rj i k j}. (*)
Действительно, при оптимальном вычислении M(i, j) последним будет перемножение ранее вычисленных матриц M(i, k) и M(k +1, j) при некотором k, таком, что i k j. Это перемножение потребует ri 1 rk rj операций умножения. В свою очередь M(i, k) и M(k +1, j) также должны вычисляться оптимальным способом, т. е. за mik и mk +1,j умножений соответственно.
Заметим, что в правой части равенства (*) разности индексов i – k и j – k –1 у слагаемых mik и mk +1,j меньше, чем разность индексов i – j в mij. Таким образом, рекуррентное соотношение (*) следует решать, начиная с mii = 0 и последовательно увеличивая разность индексов i – j, до тех пор, пока не получим m1n.
Удобно представлять результаты вычислений в виде таблицы. В этой таблице строка с номером l состоит из ячеек T(i, j), индексы которых связаны соотношением i – j = l. Таким образом, в ячейках этой строки j = i + l и T(i, j) = T(i, i + l), при этом номер i ячейки в строке принимает значения от 1 до n – l и всего в строке имеется n – l ячеек: T(1, 1 + l), T(2, 2 +l), …, T(n l, n).
| l = 0 | Т(1, 1) | Т(2, 2) | Т(3, 3) | Т(4, 4) | … | Т(n1, n1) | Т(n, n) |
| l = 1 | Т(1, 2) | Т(2, 3) | Т(3, 4) | … | Т(n2, n1) | Т(n 1, n) | |
| l = 2 | Т(1, 3) | Т(2, 4) | … | Т(n3, n1) | Т(n2, n) | | |
| … | … | … | … | … | | | |
| l = n –3 | Т(1, n2) | Т(2, n1) | Т(3, n) | | | | |
| l = n –2 | Т(1, n1) | Т(2, n) | | | | | |
| l = n –1 | Т(1, n) | | | | | | |
В ячейках таблицы T(i, j) будем хранить вычисленное значение mij и то значение qij = k в диапазоне i k j, при котором был получен Min{mik + mk +1,j + ri 1 rk rj}. Следующий алгоритм вычисляет оптимальное значение m1n и заполняет таблицу T по строкам сверху вниз:
for i := 1 to n do m[i, i] := 0; {заполнение первой строки}
for l := 1 to n –1 do
for i := 1 to n – l do
begin
j := i + l;
{заполнение T(i, j):}
m[i, j] := +;
for k := i to j – 1 do
begin
s := m[i, k] + m [k +1,j] + ri 1 * rk * rj;
if s m[i, j] then
begin m[i, j] := s;
q[i, j] := k
end { if }
end { for k }
end { for i }
Этот алгоритм требует для хранения таблицы примерно n2/2 элементов памяти и около n3/3 выполнений тела внутреннего цикла.
Рассмотрим работу алгоритма для ранее приведенного примера вычисления M1 M2 M3 M4 (размеры матриц, как указано ранее). Для заполнения строки таблицы при l = 1 вычислим последовательно
m1,2 = m1,1 + m2,2 + r0 r1 r2= 102050=10000
m2,3 = m2,2 + m3,3 + r1 r2 r3= 20501=1000
m3,4 = m3,3 + m4,4 + r2 r3 r4= 501100=5000
Здесь фактически минимум находить не требуется, т. к. тело цикла по k выполняется лишь один раз (при k = i). Заполненная строка таблицы есть
| l = 1 | m1,2 = 10000 q1,2 = 1 | m2,3 = 1000 q2,3 = 2 | m3,4 = 5000 q3,4 = 3 |
Далее для заполнения строки таблицы при l = 2 вычислим последовательно
m1,3 = Min {m1k + mk +1,3 + r0 rk r3 k = 1, 2} =
Min {m1,1 + m2,3 + r0 r1 r3 , m1,2 + m3,3 + r0 r2 r3} =
Min {0 + 1000 + 200, 10000 + 0 + 500} =
Min{1200, 10500} = 1200 (при k = 1),
m2,4 = Min {m2k + mk +1,4 + r1 rk r4 k = 2, 3} =
Min {m2,2 + m3,4 + r1 r2 r4 , m2,3 + m4,4 + r0 r2 r3} =
Min {0 + 5000 + 100000, 1000 + 0 + 2000} =
Min{105000, 3000} = 3000 (при k = 3),
и дополним таблицу следующей строкой:
| l = 1 | m1,2 = 10000 q1,2 = 1 | m2,3 = 1000 q2,3 = 2 | m3,4 = 5000 q3,4 = 3 |
| l = 2 | m1,3 = 1200 q1,3 = 1 | m2,4 = 3000 q2,4 = 3 | |
Наконец, вычислим последнюю строку таблицы, состоящую из одной ячейки Т(1, 4):
m1,4 = Min { m1k + mk +1,4 + r0 rk r4 k = 1, 2, 3} =
Min { m1,1 + m2,4 + r0 r1 r4, m1,2 + m3,4 + r0 r2 r4,
m1,3 + m4,4 + r0 r3 r4, } =
Min {0 + 3000 + 20000, 10000 + 5000 + 50000, 1200 + 0 + 1000} =
Min {23000, 65000, 2200} = 2200 (при k = 3).
Итак, вся таблица вычислена и имеет вид
| l = 0 | m1,1 = 0 q1,1 = 1 | m2,2 = 0 q2,2 = 2 | m3,3 = 0 q3,3 = 3 | m4,4 = 0 q4,4 = 4 |
| l = 1 | m1,2 = 10000 q1,2 = 1 | m2,3 = 1000 q2,3 = 2 | m3,4 = 5000 q3,4 = 3 | |
| l = 2 | m1,3 = 1200 q1,3 = 1 | m2,4 = 3000 q2,4 = 3 | | |
| l = 3 | m1,4 = 2200 q1,4 = 3 | | ||
Пока мы вычислили лишь оптимальное количество требуемых умножений m1,4 = 2200 и записали в таблицу дополнительные данные qij, позволяющие получить собственно порядок перемножения матриц. В рассмотренном примере этот порядок можно получить, начиная от финальной ячейки и продвигаясь по закрашенным ячейкам. Действительно, поскольку q1,4 = 3, то минимум был получен при m1,3 + m4,4 + r0 r3 r4. Это значит, что последним по порядку должно производиться умножение матриц M(1, 3) M4. В свою очередь, поскольку q1,3 = 1, то на предыдущем шаге минимум был получен при m1,1 + m2,3 + r0 r1 r3, т. е. матрица M(1, 3) должна вычисляться как M(1, 3) = M1 M(2, 3), а M(2, 3) = M2 M3. Это в итоге дает оптимальный порядок перемножения матриц (M1 (M2 M3)) M4.
В общем случае порядок перемножений матриц легко определить рекурсивно. Пусть имеется функция перемножения двух матриц func Mult ( A, B: Matrix): Matrix. Тогда «набросок» функции перемножения цепочки матриц выглядит так
func MatrixSeqMult ( i, j: Index): Matrix;
{i j}
global q: Tab_q;
var k: Index; var A, B: Matrix;
begin
if i j then
begin
k := q[i, j];
A := MatrixSeqMult ( i, k);
B := MatrixSeqMult ( k +1, j);
Return Mult(A, B)
end
else {i = j} Return Mi
end {MatrixSeqMult}
Доступ к матрицам Mi заданной цепочки M1 M2 ... Mn здесь не детализирован.
Полезно сравнить решение, полученное методом динамического программирования, с решением методом ветвей и границ. В рассмотренном примере возможны следующие 5 вариантов перемножения матриц M1 M2 M3 M4, а именно:
M1 (M2 (M3 M4)),
M1 ((M2 M3) M4),
(M1 M2) (M3 M4),
(M1 (M2 M3)) M4,
((M1 M2) M3) M4.
Все они получаются в дереве перебора вариантов, изображенном далее на рисунке. Некоторые узлы дерева встречаются повторно (выделены цветом). Отметим, что в методе динамического программирования повторных вычислений не производится. Вычисления проводятся так, как будто дерево сканируется снизу вверх, а результаты вычислений сохраняются в таблице и при необходимости используются.
M
(1,4)
M1 M(2,4) M(1,2) M(3,4) M(1,3) M4
M2 M(3,4) M(2,3) M4 M1 M2 M3 M4 M1 M(2,3) M(1,2) M3
M3 M4 M2 M3 M2 M3 M1 M2
Оценку количества узлов дерева в общем случае можно получить, подсчитав все возможные варианты расстановок скобок в произведении матриц. Пусть pn – число вариантов расстановок скобок в произведении n сомножителей (считая самые внешние скобки). Считая, что «последнее» по порядку умножение может оказаться на любом из n –1 мест, запишем следующее рекуррентное соотношение
pn = p1 pn –1 + p2 pn –2 + … + pn –2 p2 + pn –1 p1.
Начальным условием является p1 = 1. Далее p2 = p1 p1 = 1, p3 = p1 p2 + p2 p1 =2, p4= p1 p3 + p2 p2 + p3 p1 = 4. Оказывается, что решением этого рекуррентного уравнения являются так называемые числа Каталана pn = Сn –1, где Сn =(2 k | k) / (k +1), а запись (n | m) обозначает биномиальный коэффициент (n | m) = n!/(m! (n – m)!). При больших значениях n справедливо Сn 4n / (n sqrt( n)), т. е. число узлов в дереве перебора есть экспоненциальная функция от n.
Приведем несколько первых чисел Каталана
| n | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| Cn | 1 | 1 | 2 | 5 | 14 | 42 | 132 | 429 | 1430 | 4862 | 16796 |