Удивительно быстрые алгоритмы

  • Tutorial

Изучая программирование я встречаю примеры невозможных алгоритмов. Интуиция говорит, что такого не может быть, но компьютер опровергает её простым запуском кода. Как такую задачу, требующую минимум кубических затрат по времени, можно решить всего за квадрат? А вон ту я точно решу за линию. Что? Есть гораздо более эффективный и элегантный алгоритм, работающий за логарифм? Удивительно!


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


Интересно? Добро пожаловать под кат!


Вычисление n-го элемента реккурентной последовательности за логарифм


Под "реккурентной" я понимаю последовательность, удовлетворяющую следующему уравнению:


$ a_{n+k+1} = \sum_{i = 1}^{k} c_{i} a_{n + i} $


Первые $k$ элементов последовательности считаются заданными. Число $k$ называется мощностью последовательности, а $c_i$ коэффициентами последовательности. Типичный пример: числа Фибонначи, где $a_1 = 0$, $a_2 = 1$, $c_1 = 1$, $c_2 = 1$. Получаются всем известные числа: 0, 1, 1, 2, 3, 5, 8, 13,… Вроде бы никакой сложности нет вычислять n-ный элемент за линию, но оказывается это можно делать за логарифм!


Идея: а что, если представить вычисление $a_n$ как возведение в $n$-ю степень какого-нибудь объекта? Возводить в степень можно за логарифм, только что это будет за объект? Вообще, что нужно знать для вычисления $a_{n+2}$? Только $a_n$ и $a_{n+1}$. Значит этот объект, чем бы он ни был, должен содержать эти два числа. Ещё должен быть какой-то "естесственный" способ перемножать эти объекты. Ничего не напоминает? Очень похоже на матрицы: они составлены из чисел и их можно перемножать! Но есть ли такая матрица, перемножая которую с самой собой, мы будет последовательно получать числа Фибонначи?


Оказывается, есть! Вот она:
wake up, Neo
Можете сами проверить и убедиться, что $a_{n+1} = A^n_{1,1}$. Удивительно, но это работает!


Теперь посмотрим, сможем ли мы обобщить этот случай на любую реккурентную последовательность? Почему это сработало в случае Фибонначи? Очевидно, выполняет следующее тождество:


$$display$$ \begin{equation*} \begin{pmatrix} f_{n+1} & f_n \\ f_n & f_{n-1} \end{pmatrix} \times \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} f_{n+2} & f_{n+1} \\ f_{n+1} & f_n \end{pmatrix} \end{equation*} $$display$$


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


$ t_1 = 1 \\ t_2 = 1 \\ t_3 = 1 \\ ... \\ t_{n+3} = t_{n+2} + t_{n+1} + t_n $


За генератор можно взять матрицу $A$:


$$display$$ \begin{equation*} \begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} \end{equation*} $$display$$


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


$$display$$ \begin{equation*} \begin{pmatrix} t_{n+2} & t_{n+1} & t_n \\ t_{n+1} & t_n & t_{n-1} \\ t_n & t_{n-1} & t_{n-2} \end{pmatrix} \times \begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} = \begin{pmatrix} t_{n+3} & t_{n+2} & t_{n+1} \\ t_{n+2} & t_{n+1} & t_n \\ t_{n+1} & t_n & t_{n-1} \end{pmatrix} \end{equation*} $$display$$


Выходит, в качестве первой матрицы $T$ надо взять:


$$display$$ \begin{equation*} \begin{pmatrix} t_5 & t_4 & t_3 \\ t_4 & t_3 & t_2 \\ t_3 & t_2 & t_1 \end{pmatrix} = \begin{pmatrix} 3 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 0 \end{pmatrix} \end{equation*} $$display$$


Тогда $t_{n+2} = (T \times A^{n - 1})_{1, 1}$. Благодаря ассоциативности матричного умножения матричного умножения, мы можем сначала вычислить возведение $A$ в $n - 1$ степень, а потом уже перемножить $T$ и $A$. В случае Фибонначи $T$ совпадало с $A$.


Получается, общий вид генератора для реккурентной последовательности мощности $k$ будет иметь вид:


$$display$$ \begin{equation*} \begin{pmatrix} c_1 & 1 & 0 & 0 & \cdots & 0 \\ c_2 & 0 & 1 & 0 & \cdots & 0 \\ c_3 & 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ c_{k-1} & 0 & 0 & 0 & \cdots & 1 \\ c_k & 0 & 0 & 0 & \cdots & 0 \end{pmatrix} \end{equation*} $$display$$


Начальная матрица, на которую умножается генератор, содержит первые $2k - 1$ элементов последовательности, расположенные в "пирамидном" порядке.


Понятно, что такой алгоритм будет работать за $O(k^3 \log n)$ по времени: каждое умножение матриц выполняется за $O(k^3)$, а всего $O(\log n)$ таких умножений. Кубический множитель играет важную роль. Например, в случае Фибонначи, преимущество нашего алгоритма над линейным достигается только при $n \ge 44$, а в случае Трибонначи вообще при $n \ge 208$. Однако значение множителя можно понизить если в рекурсивном алгоритме возведения в степень делить не на два случая чёт\нечёт как обычно, а на три. Тогда в основании логарифма будет тройка, что понизит преимущество до $n \ge 118$.


Следующий класс Matrix содержит методы умножения и возведения в степень квадратных матриц:


class Matrix:
    def __init__(self, n):
        self.n = n
        self.rows = [[0 for col in range(n)] for row in range(n)]

    def set(self, row, col, value):
        self.rows[row][col] = value

    def get(self, row, col):
        return self.rows[row][col]

    def __str__(self):
        result = ''
        for row in self.rows:
            result += ' '.join([str(col) for col in row])
            result += '\n'
        return result

    def __mul__(self, other):
        result = Matrix(self.n)
        for row in range(self.n):
            for col in range(self.n):
                s = sum([self.get(row, k) * other.get(k, col) for k in range(self.n)])
                result.set(row, col, s)
        return result

    def __len__(self):
        return self.n

    def __pow__(self, k):
        if k == 0:
            result = Matrix(len(self))
            for i in range(len(self)):
                result.set(i, i, 1)
        elif k == 1:
            result = self
        elif k == 2:
            result = self * self
        else:
            rem = k % 3
            prev = self.__pow__((k - rem) // 3)
            result = prev * prev * prev
            if rem:
                result *= self.__pow__(rem)
        return result

Внимание на метод __pow__: это определение операции M ** k, где M это объект класса Matrix. Он вычисляет её, рекурсивно разбивая каждый вызов на три кейса в зависимости от остатка от деления на 3. Это позволяет уменьшить константу времени работы алгоритма.


Определим наши $T$ и $A$ для Трибонначи, пользуясь Matrix:


A = Matrix(3)
A.set(0, 0, 1)
A.set(0, 1, 1)
A.set(1, 0, 1)
A.set(1, 2, 1)
A.set(2, 0, 1)
T = Matrix(3)
T.set(0, 0, 3)
T.set(0, 1, 1)
T.set(0, 2, 1)
T.set(1, 0, 1)
T.set(1, 1, 1)
T.set(1, 2, 1)
T.set(2, 0, 1)
T.set(2, 1, 1)
T.set(2, 2, 0)
n = int(sys.argv[1])
if n:
    print(T * A ** (n - 1))
else:
    print(T ** 0)

Этот код принимает аргумент $n$ — номер элемента последовательности Трибонначи, который мы хотим вычислить. Вот и всё, добавлю только, что такой алгоритм можно обобщить на нечисловые последовательности, нужно только чтобы элементы принадлежали какому-нибудь Кольцу. Упрощённо, это математическая структура, в которой определены сложение и умножение над произвольным множеством объектов.


Вычисление подмассива наибольшей длины за линию


Вторая задача звучит так: дан массив целых чисел A[1..n] (для простоты считаем что индексы идут с единицы). Надо найти подмассив A[i..j] с наибольшей суммой. В подмассиве присутствуют все элементы от i до j включительно. Эта задача интересна тем, что для неё существуют алгоритмы с временной сложностью $O(n^3)$, $O(n^2)$, $O(n \log n)$ и даже $O(n)$.


Вкратце идея каждого алгоритма:


  1. Куб. Пройтись по всем подмассивам, каждый раз высчитывая их сумму заново. Получается три вложенных цикла. Самый простой, но самый медленный вариант.
  2. Квадрат. Пройтись по всем подмассивам, сохраняя текущую сумму.
  3. $O(n \log n)$. Разбить массив пополам, рекурсивно рассмотреть три случая: искомый ответ полностью лежит в левой части, между ними, либо в правой.
  4. $O(n)$. Завести новый массив T[1..n], где i-тый элемент равен наибольшей сумме подмассива, оканчивающегося в i. Оказывается, можно считать $T$ за линию, а конечным ответом будет наибольшее из $T$. Это так называемый алгоритм Кадана

Нам интересен последний вариант. Давайте подумаем, как найти T[i + 1], зная T[i]? Если мы знаем сумму наибольшего подмассива, оканчивающегося в i, то следующий элемент может либо расширять его, либо начать новый подмассив. Действительно, T[i + 1] может быть либо T[i] + A[i + 1], либо A[i + 1], либо 0, если A[i + 1] < 0. Получаем реккурентную формулу:


T[0] = 0,
T[i + 1] = max{T[i] + A[i + 1], A[i + 1], 0} = max{T[i] + A[i + 1], 0}

Докажем последнее равенство. Понятно, что T[i] >= 0 для любого i. Пусть k = A[i + 1]. Рассмотрим три случая:


  1. k < 0. Тогда 0 превзойдёт k в первом max.
  2. k = 0. В первом max можно просто убрать второй аргумент.
  3. k > 0. Тогда max{T[i] + k, k, 0} = T[i] + k = max{T[i] + k, 0}.

Благодаря линейности и простоте уравнения, алгоритм получается совсем коротким:


def kadane(ints):
    prev_sum = 0
    answer = -1
    for n in ints:
        prev_sum = max(prev_sum + n, 0)
        if prev_sum >= answer:
            answer = prev_sum
    return answer

Заключение


В обоих задачах качественно повысить производительность нам помогла техника динамического программирования. Это не случайно, динамика часто даёт асимптотически оптимальные алгоритмы благодаря встроенной экономии: мы считаем всё только один раз.


А какие вы знаете удивительные алгоритмы?

Реклама
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее

Комментарии 28

    +1
    А зачем матрица T? Тут же вроде всегда было вектора достаточно…
      +3
      А какие вы знаете удивительные алгоритмы? — Знаем бессмысленный и беспощадный нормальный алгоритм Маркова для деления чисел habr.com/ru/post/110004
        +3
        Однажды у меня было удивление при чтении математического справочника 1950-1960-х годов. Там было много хаков, позволяющих просто посчитать вещи, которые я без компьютера не посчитаю.
          +1

          А справочник не подскажете?)

            0
            Давно у меня его нет. И уже тогда он был без обложки. Толстая книга размером А5. Характерной чертой эпохи была инструкция по работе с механическими счетными машинами.
          +1
          Алгоритм сортировки массива за линию Radix sort сильно ломает шаблоны
            +3
            Ну, он довольно условно за линию: он сортирует за O(N*число_разрядов). Если предположить, что все числа разные – то нам надо по крайней мере log(N) разрядов для их представления, и получаем всё тот же O(N*log(N)).
              0
              Мы же можем сортировать не поразрядно а побайтого, тоесть int64 будем сортировать по 8 байтам, в этом случае получим O(N). Другое дело, что на практике там слишком большой постоянный множитель + нужна доп память, сортировки сравнением оказываются быстрее по итогу.
                +2
                Получите N*log256(N), то же самое :-)

                На самом деле это один из примеров, показывающий, что асимптотики порой недостаточно для оценки. Допустим, у нас не больше миллиона чисел – тогда log256(N)<=2.5, log2(N)<=20. Т.е. логарифм можно смело заменять на (не такой уж большой) постоянный множитель.
                Ещё прикольнее алгоритмы, в сложности которых фигурирует обратная функция Аккермана, например, ru.wikipedia.org/wiki/Алгоритм_Краскала – она вроде и растёт от N, но настолько медленно, что проще забить.
            +2

            Так-то и матрицы можно быстрее чем за куб перемножать.

              0
              Как вариант, но константа всё равно будет значительной.
              +1

              С планшета часть формул не отрендерилась, не в первый раз уже

                0
                с компа тоже.
                ни одна
                +1
                Это хорошие алгоритмы, но это же не динамическое программирование.
                  0
                  Почему? Кадан чистая динамика
                  +2

                  С (целочисленными) последовательностями и Фибоначчи в частности есть засада, они быстро выходят за размер целого числа. Сотое число Фибоначчи уже не влезает в 64 бита. Так что если считать сложность алгоритма честно, то нужно учитывать что умножение чисел выполняется не за константу, а становится дороже при росте чисел :). Получится уже не логарифм.

                    +8

                    Мой любимый пример — алгоритм Шонхаге-Штрассена (?) для умножения больших чисел через преобразование Фурье за O(NlogN). Составляете последовательность из цифр числа A, применяете FFT, то же для числа B, потом их (FFT) перемножаете попарно, применяете обратное преобразование (IFFT), результат – цифры произведения А и B.


                    Пожалуй, самое сильное "ВАУ" за всю мою программистскую карьеру.


                    Деталей не помню, выглядит как-то так:


                    import numpy as np
                    
                    A = [0, 0, 0, 0, 0, 0, 0, 2]  #   2
                    B = [0, 0, 0, 0, 0, 1, 3, 2]  # 132
                    
                    FC = np.fft.fft(A) * np.fft.fft(B)
                    C = np.fft.ifft(FC)  # <-- Тут должна быть какая-то возня с С, типа переносов между разрядами
                    
                    print(" ".join([str(round(abs(c))) for c in C]))  # Prints "0 0 0 0 2 6 4 0"
                      0
                      Возможно стоило упомянуть о том, что некоторые рекуррентные последовательности можно вычислить за O(1). Например, произвольный член ряда вариации Фибоначчи, который начинается не с {1,1,2,...}, а с {a,b,a+b,...} можно вычислить через функцию

                      (Power((1+Sqrt(5))/2, n)*((-5+3*Sqrt(5))*a - (-5+Sqrt(5))*b) - Power(2/(1+Sqrt(5)), n)*((5+3*Sqrt(5))*a - (5+Sqrt(5))*b)*Cos(n*Pi))/10
                        0
                        Видимо вы имеете в виду Формулу Бине, которая также опирается на возведение в степень, и если эту операцию мы считаем O(1), то и представленный алгоритм тоже O(1), но вообще это некорректно.
                          0
                          Нет — в данном случае формула Бине является лишь частным случаем при a=1 и b=1, и, к слову, выводится она не через возведение в степень, а через производящие функции.
                            0
                            Выводится через производящие функции, но для подсчета конкретного числа Фибоначчи с её помощью нужно вычислять два возведения в степень
                              0
                              А возведение в степень — через умножение логарифмов. А логарифм можно посчитать даже на линейке.
                                +1
                                Сотрите, у меня было два основных тезиса
                                1) Вычисление степени матрицы можно свести к формулам наподобие того, что вы выписали
                                Диагонализация матрицы из статьи


                                умножая матрицу в степени на вектор (1, 1) получается последовательность Фибоначчи, умножая на (a, b) получается вариация, о которой Вы писали.
                                2) Некорректно считать, что такое возведение в степень — это O(1). Чтобы вычислить a^n по формуле exp(log(a) * n) нужно обеспечить достаточную точность при работе с вещественными числами, эта точность зависит от размера результата, точнее сказать не могу — не разбирался.
                                  0
                                  1) давайте пойдём дальше — как свести вычисления с матрицами к дробному отрицательному n или решить прямо противоположную задачу — зная значение члена ряда Фибоначчи найти его порядковый номер.

                                  2) O(1) значит независимость от n, а не сложность вычисления отдельно взятой операции для вычисления конечного результата. К тому же, в вещественных числах вычисление логарифма не зависит (сильно, какие-то колебания конечно же есть) от аргумента, это совершенно точно (а я разбирался и даже писал свою реализацию).
                                    0
                                    К тому же, в вещественных числах вычисление логарифма не зависит от аргумента

                                    Если Вы говорите об операциях скажем в 64-битном double или любом другом конечном типе хоть с плавающей, хоть с фиксированной точкой, то я согласен — вычисление логарифма не зависит от аргумента. Проблема в том, что в этих вычислениях будет погрешность, которая ограничена снизу в силу конечности представления. И этой точности не будет достаточно для вычисления произвольного числа Фибоначчи, так как последовательность Фибоначчи не ограничена, а следовательно не ограничено и количество значащих цифр в представлении чисел Фибоначчи.
                                      –1
                                      Даже избежав ограниченности числа с плавающей точкой, вы всё равно упрётесь в ограниченность памяти в компьютере, всех компьютеров мира, затем в ограниченность планеты Земля и всей галактики в целом, чтобы записать очередное вычисленное число.
                                        0
                                        Вроде бы Вы как раз говорите о том, что зависимость от n есть.
                        +5
                        «За линию» очень режет слух, почему нельзя было немного уточнить: «с линейной сложностью»?
                        Заголовок спойлера
                        image

                        Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                        Самое читаемое