Алгоритм Тадао Такаока для нахождения максимальной подматрицы или Maximum Subarray Problem

Не так давно прошёл конкурс параллельного программирования Acceler8 2011. Суть задачи заключалась в поиске максимальной подматрицы в данной матрице (сумма элементов найденной подматрицы должна быть максимальной). После недолгого «гугления» было найдено, что некий алгоритм Тадао Такаока решает эту задачу быстрее других.

«Вызов принят!», и я начала искать этот алгоритм везде, где только можно, задавшись целью реализовать его. Не смотря на то, что распараллеливается он плохо и в своей сложности содержит немаленькую константу.

Однако всё, что удалось найти, — статьи на английском этого самого Тадао Такаоки (вот одна из этих статей). Пришлось переводить.

Сама идея алгоритма сначала казалась до безобразия простой:

Для начала нам необходимо посчитать все префиксные суммы s[i, j] для таких матриц, как a[1..i, 1..j], для всех i и j с ограничением s[i, 0] = s[0, j] = 0. Очевидно, это делается за время O(mn).
Основной алгоритм:
Если матрица оказывается одним элементом — возвращаем его значение.
Иначе:
  • если m > n, поворачиваем матрицу на 90 градусов.
Таким образом m ≤ n.
  • Пусть A_left — решение для левой половины.
  • A_right — решение для правой половины.
  • A_column — решение для column-centered задачи (нахождение такой максимальной подматрицы, которая бы пересекала центральную линию)
  • Общее решение — максимум из этих трёх.

Если A_left и A_right находятся через рекурсию, то находить A_column надо вручную.

Вот, что написано по поводу нахождения этого самого решения в оригинальной статье:

"Теперь column-centered задача может быть решена следующим образом.
A_column = max(k=1:i-1,l=0:n/2-1,i=1:m,j=n/2+1:n) {s[i, j] − s[i, l] − s[k, j] + s[k, l]}.
Мы фиксируем i и k, и максимизируем оставшееся выражение изменением l и j. Тогда это эквивалентно следующему выражению:
i = 1, ..., m и k = 1, ..., i − 1.
A_column [i, k] = max(l=0:n/2-1,j=n/2+1:n) {−s[i, l] + s[k, l] + s[i, j] − s[k, j]}
Пусть s∗[i, j] = −s[j, i]. Тогда данное выражение может быть сведено к:
A_column [i, k] = −min(l=0:n/2-1) {s[i, l] + s∗ [l, k]} + max(j=n/2+1:n) {s[i, j] + s ∗[j, k]}
"

Далее следует упоминуть о некой другой задаче Тадао Такаока: Distance Matrix Multiplication (не берусь переводить этот термин).

Суть такова:

Целью DMM является вычисление произведения расстояний(distance product) C = AB для двух n-мерных матриц A = а[i,j] and B = b[i,j], чьи элементы — вещественные числа.

c[i,j] = min(k=1:n) { a[i,k] + b[k,j] }

Суть c[i,j] — наименьшее расстояние из вершины i из первого слоя до вершины j в третьем слое в ациклическом ориентрированном графе, состоящем из трёх слоёв вершин. Эти вершины пронумерованы 1,...,n в каждом слое, и расстояние от i в первом слое до j во втором слое — a[i,j] и из i во втором слое до j в третьем слое — b[i,j]. Очевидно, что это выражение можно легко привести для подсчёта максимума вместо минимума — это будет задача о нахождении наибольших путей.

Так вот воспользовавшись этим определением можно получить следующее:
В формуле A_column [i, k] = −min(l=0:n/2-1) {s[i, l] + s∗ [l, k]} + max(j=n/2+1:n) {s[i, j] + s ∗[j, k]} видно, что первая часть формулы — DMM для минимума, вторая — DMM для максимума. Пусть S1 и S2 матрицы, чьи элементы (i, j) являются s[i, j − 1] s[i, j + n/2] соответственно. Для любой матрицы T, пусть T∗ — матрица, полученная из T транспонированием и отрицанием. Тогда верхнее выражение может быть посчитано «перемножением» S1 и S1∗ для минимума и получением нижне-треугольной матрицы, «перемножением» S2 и S2∗ для максимума и получением нижне-треугольной матрицы и, наконец, вычитанием из последней предыдущей матрицы. Вычислив в ней максимум, получим ответ:
A_column = max(DMM_max(S2*S2∗)-DMM_min(S1*S1∗)).

Далее преведен псевдокод для решения задачи DMM для матриц A и B размером nxn:

Copy Source | Copy HTML
  1. /Отсортировать n строк B в порядке убывания в списки list[i], т.о. в каждом списке list[i] будут находится индексы упорядоченных элементов строки i;
  2. /Заведём множество V = {1, ..., n};
  3. for i := 1 to n do begin
  4.   for k := 1 to n do begin
  5.     cand[k]:=first of list[k];
  6.     d[k] := a[i, k] + b[k, cand[k]];
  7.   end;
  8.   / Упорядочим V как очередь с приоритетом по ключам d[1], ..., d[n];
  9.   / Множество решений S - пустое;
  10.   /* Фаза 1 : До критической точки */
  11.   while |S| ≤ n − n log n do begin
  12.     /Найти v с минимальным ключем в очереди;
  13.   /Занести cand[v] в множество S;
  14.     c[i, cand[v]] := d[v];
  15.     /Пусть множество W = {w|cand[w] = cand[v]};
  16.   for w in W do
  17.       while cand[w] is in S do cand[w]:= next of list[w];
  18.   /Отсортировать очередь W по новым ключам d[w] = a[i, w]+b[w, cand[w]];
  19.   end;
  20.   U := S;
  21.   /* Фаза 2 : После критической точки */
  22.   while |S| < n do begin
  23.     /Найти v с минимальным ключём в очереди;
  24.     if cand[v] is not in S then begin
  25.       /Занести cand[v] в множество S;
  26.       c[v, cand[v]] := d[v];
  27.       /Пусть множество W = {w|cand[w] = cand[v]};
  28.     end else W = {v};
  29.     for w in W do
  30.       cand[w]:=next of list[w];
  31.     while cand[w] is in U do cand[w]:= next of list[w];
  32.     /Отсортировать очередь W по ключам d[w] = a[i, w]+b[w, cand[w]];
  33.   end;
  34. end.

В статье рекомендуется использовать бинарную кучу для очередей. В более поздних версиях было указано, что использование фиббоначиевой кучи могло бы ещё больше ускорить процесс, однако я, к своему стыду, не смогла её реализовать в полной мере. Поэтому использовала бинарную.

Однако даже реализовав этот псевдокод, я обнаружила, что решение не работает. Очевидная ошибка была найдена на этапе вычисления S1 и S2 — описанного правила для построения этих матриц было недостаточно для правильной работы алгоритма. И ни в одной статье я более подробного описания не нашла. Таким образом строить матрицы пришлось самой.

Итак, после долгих часов, потраченных на перебор и проверку, были получены следующие правила.

Пусть S — префиксная матрица для матрицы A.

Тогда матрица S1 размерами m на n/2 содержит в себе первый столбец из 0, а оставшиеся n/2-1 столбцы содержат первые столбцы префиксной матрицы S:


Матрица S1∗ размерами n/2 на m содержит первую строку и столбец из 0, оставшаяся подматрица представляет собой отрицание и транспонирование первых строк и столбцов префиксной матрицы S:

(k = n div 2 + n mod 2)

Матрица S2 размерами m на k содержит в себе столбцы префиксной матрицы S, начиная с k-го:


Матрица S2∗ размерами k на m: первый столбец из 0, далее — транспонирование и отрицание столбцов префиксной матрицы S, начиная с k-го столбца:


В целом алгоритм оказался очень сложным для понимания и реализации и не оправдал своих ожиданий на исключительную скорость(никаких дополнительных оптимизаций не осуществлялось).

В той же статье утверждается, что данный алгоритм работает за время O(n^3*log(log(n))/log(n)), а также упоминается, что алгоритм работает быстрее в сравнении с другими алгоритмами на матрицах большого размера (видимо, очень большого).

Если кому-то будет интересно или жизненно необходимо, здесь лежат исходники моей реализации (написано на плюсах, поэтому «понятности» не гарантирую). Буду очень рада, если эта статья кому-то поможет (не зря же я потратила столько времени на разбирательства с этим алгоритмом!).

Similar posts

AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 29

    +3
    > Не смотря на то, что распараллеливается он плохо

    Отлично распараллеливается, причём на нескольких уровнях: как task-based нерегулярный параллелизм на уровне left/right/column подзадач, так и сам distance matrix multiplication.
      +1
      Да? Правда интересно, как его можно было распараллелить, чтобы достичь лучших рез-тов, чем от эвристик, которыми все в итоге пользовались? Потому как у нас локально он вроде проигрывал альтернативному (с эвристиками).

      Плюс ко всему, забыла упомянуть, обнаружились огромные проблемы с памятью, которой катастрофически не хватало уже на средних матрицах (до больших не дошли).

      Хотя может все эти вопросы и решались бы, будь у нас больше времени (поздно вспомнили про конкурс). (=
        +1
        Тот факт, что на этих размерностях Кадан оказался быстрее, не говорит о том, что Такаока плохо распараллеливается.

        Мы тоже писали Такаоку bitbucket.org/gribozavr/max_subarray_2d и у нас вот такие результаты на тестах с русскоязычного форума:
        > 5000x5000 — 14.8 s на 40 ядрах, 517 s последовательно. маленький и средний — 27.8 s на 40 ядрах, 1092 s последовательно. Huge не пробовали.
          +1
          Извините, я имела ввиду в сравнении с Каданэ.
          Спасибо за ссылку, в свободное время обязательно посмотрю! Очень уж интересно! (=
          0
          >>Правда интересно, как его можно было распараллелить, чтобы достичь лучших рез-тов, чем от эвристик
          Тут проблема была не столько в распараллеливании с целью снижения итоговой сложности в расчете на ядро, а хотя бы в чисто аппаратных вопросах доступа к памяти. Не все знали, а уж тем более учитывали NUMA организацию, «бегая» всеми 40 потоками по матрице. А ведь при «честном» способе и полном хранении матрицы, на нее могло приходиться полгора гигабайта. Кеши, false sharing, тайминги памяти… Не сравнить с решением, которому нужны всего несколько мегабайт (к примеру).
        0
        «В целом алгоритм оказался очень сложным для понимания и реализации и не оправдал своих ожиданий на исключительную скорость(никаких дополнительных оптимизаций не осуществлялось).»

        Ты писала какой-нибудь другой алгоритм или в сравнении с чем он «не оправдал ожиданий»? Если был реализован ещё другой алгоритм — было бы интересно посмотреть разницу в папугаях.
          +1
          Был. (= Его писала не я, а мой сокомандник. И по сравнению с тем, который был гораздо проще, но использовал эвристики, связанные со спецификацией строения матриц, Такаока был не пригоден. Его не получилось как следует распараллелить и памяти он «жрал» очень много. )=
            +1
            Я слышал, что там задача упрощалась по сложности за счёт специфики формирования матриц. Мне кажется, сравнение чистого алгоритма с задаче-специфичным несправедливо.
              +3
              Не спорю, мне вообще жизнь казалась несправедливой, после потраченного времени (вместо учёбы) и осознания того, что едва ли мой алгоритм можно будет использовать в конкурсе. Но всё-равно разобраться с ним было забавно (=

              Зато как удачно ложилась реализация с Каданэ! Сколько там можно было эвристик и оптимизаций воткнуть! (= Жаль только, что времени не хватило.
                0
                Да, вместо алгоритма O(n*n*m) и подобных были решения вплоть до O(n*p), где p часто в пределах n*N, с малым натуральным N, а то и вообе меньшим чем m. Т.е. период не более нескольких строк :)
                  0
                  Видимо где-то в ваших словах закралась досадная ошибка, так как если p = n*N, при этом если рассматривать ситуацию обратную «а то и вообе меньшим чем m» (это слово «вообще»?), то получаем, что N может быть больше m. Тогда получаем p = n*(m+c) => O(n*n*(m+c)), что получается даже сложнее изначальной вашей оценки O(n*n*m).
                    0
                    p=m*N. Несколько строк, а не столбцов.

                    Да, «вообще». Извиняюсь за опечатку.
                      +2
                      Если кому-то будет интересно, тут решения всех команд с их же описаниями решений.

                      Вообще по Такаока сложилось аналогичное мнение, что алгоритм для матриц «слегка» бОльших, нежели предлагались в конкурсе (максимум 20k на 20k).
                        0
                        О! Хорошая компания собирается ;)
            +6
            К сожалению, не могу идентифицировать вас по нику на ISN и вспомнить название вашей команды… Рабочий комп не под рукой, да и прайваси…

            Но я верю, что жизнь — справедлива! Мне победитель конкурса говорил, что они потратили с сокомандником на полировку решения 3 (три!) недели в выделенном режиме.

            Мы ведь не пытаемся конкурировать с олимпиадным программированием, стараемся поставить во главу угла тонкую оптимизацию решения под конкретные платформы. Таким образом, хорошая стратегия для конкурсантов выглядит примерно так: 1) оценить быстродействие последовательного алгоритма 2) оценить возможности его распараллеливания 3) чем больше времени останется нв профилирование и низкоуровневую оптимизацию — тем лучше!

            Теперь посмотрим сюда: только 12 из 66 команд смогли корректно решить все матрицы, без учета скорости. Больше половины не уложились в лимит времени в одном или нескольких тестах. А ведь большинство использовало хорошо известный алгоритм.

            В отличие от большинства вы попробовали что-то новенькое, да и «разобраться с ним было забавно» :)

            А насчет победы — не переживайте, будут еще конкурсы… С наступившим старым новым годом, кстати.
              +1
              >>Теперь посмотрим сюда: только 12 из 66 команд смогли корректно решить все матрицы, без учета скорости.
              Вот лимит по времени работы там все-равно был :) Не успел — не правильно. По сути это и был учет скорости в неявном виде.
                0
                Согласен. Я пытался сформулировать мысль о важности как самого алгоритма, так и его оптимизации.
                +1
                Мы в рейтинге заняли кругленькое 10е место.
                Трёх недель выделить не получилось, поэтому и не сорвали джекпот ;)
                А вообще мне действительно было гораздо интереснее разобраться с этим «неподходящим» алгоритмом, чем с каданэ (который придумали за 10 минут) и всякими изощрёнными оптимизациями — не моё это. Поэтому я позорно дезертировала ковыряться с Такаокой, оставив сокомандника одного париться с эвристиками (=

                Главное не победа, а участие (=

                Спасибо, вас так же с наступившим (=
                  +2
                  Зато вы успеваете писать на Хабр, как я понимаю — в самую сессию :)
                    +1
                    Угу, завтра экзамен по матстату =/…
                      +2
                      Ни пуха! :)
                        +1
                        К чёрту (=
                0
                Кто-нибудь может мне объяснить, почему этот алгоритм напомнил мне расстояние Левенштейна?
                  0
                  Правильно ли я понял, что данный алгоритм предназначен для поиска максимальной подматрицы в рядом стоящих столбцах/строках?
                  Видели ли Вы где-то упоминание о его применении к задаче, когда допускается перестановка строк/столбцов?
                    +2
                    Да, вы правильно поняли задачу.
                    Встречный вопрос: что имеется в виду под вашим «допускается перестановка строк/столбцов»? Если это абсолютно любые перестановки в неограниченном количестве — то так можно переместить любой элемент матрицы на любую позицию, следовательно максимальная подматрица состоит из всех положительных элементов, собранных в прямоугольник. Или я что-то не так поняла?
                      +1
                      По-моему, не всегда это возможно.
                      Рассмотрим матрицу вида:
                      |1|1|0|
                      |1|0|1|
                      |0|1|1|

                      Посчитайте ее детерминант, он отличен от нуля. Если бы можно было переставить, получив блок из единиц размером 2x3, то детерминант был бы равен нулю (ибо ранг матрицы был бы 2).

                      Под перестановкой я имел в виду следующее. Нам дают матрицу. Мы можем сгенерировать любую перестановку(подстановку, permutation) для строк и столбцов, после чего задача в общем-то сводится к подряд стоящим столбцам/строкам. Вот только надо выбрать именно такую перестановку, чтобы последующее решение задачи было максимальным.
                    0
                    Извиняюсь, не нашел в правилах конкурса. Есть ли ограничения на размер используемой памяти?
                      0
                      В российском варианте практическое ограничение было в том, что программа будет запущена также на десктопной машине. На MTL же памяти 64 Gb, поэтому можно считать что не ограничено.
                        0
                        спасибо

                    Only users with full accounts can post comments. Log in, please.