Ant colony optimization method for solving parametric problems on vector SIMT accelerators

Full Text

Введение и обзор литературы

Задачи анализа и синтеза сложных систем основываются на решении проблем оптимизации их параметров по заданным критериям. В случае однокритериальной оптимизации широко применяются градиентные методы, обеспечивающие наискорейший поиск оптимальных значений. Однако градиентные методы демонстрируют высокую эффективность преимущественно на унимодальных функциях с достаточно простой формой. Для многоэкстремальных функций со сложной структурой эти алгоритмы, как правило, находят лишь локальный экстремум. Поиск глобального экстремума требует привлечения метаэвристических алгоритмов, которые стохастическим образом осуществляют параллельное исследование множества возможных решений [1]. Эволюционное происхождение данных эвристик предполагает их повышенную сходимость. Эвристическая природа таких алгоритмов обеспечивает простоту их модификации и адаптации к различным оптимизационным задачам, а возможность учета эвристической информации позволяет учитывать специфику решаемой проблемы, что порождает множество модификаций и практических применений [2].

В настоящей работе рассматриваются различные модификации метаэвристического бионического метода оптимизации муравьиных колоний [3]. Метод оптимизации муравьиных колоний (ACO, Ant Colony Optimization), первоначально предложенный M. Dorigo для решения задачи коммивояжера, в дальнейшем был распространен на широкий класс оптимизационных задач [3, 4]. Помимо классической задачи коммивояжера он применяется к задаче нескольких коммивояжеров (mTSP) [5], задачам размещения оборудования (QAP) [6, 7]. Основу метода составляет вероятностный выбор маршрутов муравьями-агентами с последующей оценкой найденных путей и коррекцией вероятностей за счет обновления феромонных следов на вершинах графа. Помимо вероятностного выбора (этап исследования) возможен переход в вершину с максимальной вероятностью (этап эксплуатации), что позволяет проводить поиск решений в окрестностях рациональных маршрутов. Хотя оригинальный алгоритм муравьиной колонии показал обнадеживающие результаты в решении задачи коммивояжера, несбалансированность этапов исследования и эксплуатации приводит к проблеме стагнации, когда все муравьи-агенты начинают следовать по одному и тому же маршруту. В процессе работы алгоритма осуществляется либо исследование новых областей пространства поиска, либо эксплуатация соседних высококачественных решений [8]. Если доминирует исследование, алгоритм может анализировать малоперспективные области, а при избыточной эксплуатации возможна преждевременная сходимость к субоптимальному результату. Различия между алгоритмами ACO заключаются в способах управления балансом между исследованием и эксплуатацией.

Модификации ACO, направленные на балансировку процессов исследования и эксплуатации, принято классифицировать на следующие группы: управление феромонами, параметризация и гибридизация [9]. К первой группе относятся различные механизмы, используемые в алгоритмах оптимизации муравьиных колоний, такие как элитизм в элитной системе муравьев (EAS), обновление следа в системе колонии муравьев (ACS), обучение на основе Q-обучения (AntQ), ранжирование в системе муравьев на основе рангов (RAS), ограничение следов в MMAS и вычитание феромонов в системе «лучший–худший» (BWAS) и другие [9]. Ко второй группе относятся методы, в которых исследование и эксплуатация поддерживаются за счет параметризации, при которой значения параметров алгоритма ACO изменяются в онлайн- или офлайн-режиме. Офлайн-методы реализуются либо по схеме проб и ошибок, либо с использованием методов машинного обучения [10]. Оба подхода применяются до запуска алгоритма и близки к методологии мультистарта. Онлайн-методы предполагают изменение значений параметров в процессе выполнения алгоритма, что существенно усложняет его структуру и снижает универсальность. К третьей группе относятся успешные комбинации базовых техник управления феромонами с алгоритмами локального поиска, которые значительно повышают качество решений, получаемых с помощью ACO. 
Применение локального поиска, такого как 2-opt, в сочетании с ACO продемонстрировало свою эффективность; также были предложены гибридные подходы, интегрирующие ACO с другими методами оптимизации. Наибольшее распространение получили комбинации генетического алгоритма и метода муравьиных колоний, а также метода муравьиных колоний и метода роя частиц [11, 12].

В методе муравьиных колоний муравьи-агенты перемещаются группами по $N$ агентов. Поскольку граф в процессе перемещения не изменяется, их движение может осуществляться параллельно. Параллельные алгоритмы подразделяются на мелкозернистые, когда каждый муравей-агент перемещается в отдельном потоке, и крупнозернистые, когда в одном потоке функционируют независимые колонии муравьев, взаимодействующие через специальные механизмы. Как правило, мелкозернистые алгоритмы уступают крупнозернистым из-за простоты алгоритма движения отдельного муравья-агента и высоких накладных расходов, связанных с управлением параллельными процессами. Для сокращения объема межпроцессорных коммуникаций был разработан метод частично асинхронной параллельной реализации (PAPI), в котором обмен информацией между колониями происходит после фиксированного числа итераций [13]. Подход «главный/подчиненный» к параллельной реализации метода муравьиных колоний для решения задач QAP и TSP использует концепцию главного потока, который агрегирует решения от параллельных колоний муравьев и транслирует информацию для обновления матриц феромонов в каждой колонии. Для реализации данного подхода необходимо, чтобы каждый поток, соответствующий отдельной колонии муравьев, обладал собственными данными и алгоритмами [14]. Дальнейшее развитие метода муравьиных колоний было связано с прогрессом вычислительной техники и появлением GPU-ядер (Graphics Processing Unit) NVIDIA, поддерживающих параллельные вычисления на базе технологии CUDA. Поскольку графические ускорители обеспечивают более высокую пиковую производительность по сравнению с многоядерными процессорами и отличаются доступной стоимостью, исследователи проявляют повышенный интерес к распараллеливанию алгоритма ACO именно на GPU, а не на центральных процессорах [15, 16]. Данный тип параллельных алгоритмов ACO основан на модели параллельного обмена данными (DPAM). Основная идея DPAM заключается в сопоставлении каждого муравья с группой потоков с общей памятью, например, блоком потоков CUDA. Наибольшее распространение параллельные реализации ACO получили для решения задачи коммивояжера. Назначая каждый поток одному или нескольким узлам графа, потоки в блоке могут совместно вычислять правило перехода между состояниями, тем самым повышая параллелизм на уровне данных [17]. При этом возможность эффективного распараллеливания на GPU для решения задач TSP и QAP ограничена необходимостью частого взаимодействия с общими данными, представленными в виде специализированных структур.

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

Подобные параметрические задачи также решались с применением метода муравьиных колоний. Отдельное внимание уделено комбинаторной оптимизации [11, 12] и ее применению для задач составления расписаний [18]. В таких задачах поиск решений модификациями ACO осуществляется в параметрическом графе (рис. 1), где каждому параметру соответствует множество допустимых значений (вершин графа), и метод муравьиных колоний выбирает ровно по одному значению для каждого параметра. Параметрический граф использовался, в частности, для формирования годовой заявки на производство запасных частей летательного аппарата в системе послепродажного обслуживания [19]. В области прогнозирования локализации месторождений полезных ископаемых предложен гибридный подход для оптимизации гиперпараметров метода опорных векторов (SVM), сочетающий алгоритм оптимизации муравьиных колоний и метод сеточного поиска (GS), который систематически оценивает все комбинации гиперпараметров для нахождения оптимальной конфигурации модели [20].

 

Рис. 1. Структура параметрического графа
[Figure 1. Parametric graph structure]

 

В работе исследуется разработанная авторами матричная модификация ACO, предназначенная для выполнения параллельных вычислений при решении параметрических задач. Предложенная реализация параллельной матричной модификации ACO ориентирована на использование SIMT-ускорителей (Single Instruction, Multiple Thread) с целью поиска оптимума целевой функции. Поскольку вычисление целевой функции в модели может быть ресурсоемким, возникает необходимость в организации промежуточного хранилища результатов и соответствующей адаптации алгоритма ACO для работы с этим хранилищем. В качестве базовых для матричного представления в работе выбраны модификации ACOCNI и ACOCCyI, которые обеспечивают улучшенную стратегию исследования за счет применения хэш-таблицы [21, 22]. В алгоритме ACOCCyI все найденные решения и соответствующие им значения целевых функций записываются в хэш-таблицу для организации поиска нового решения каждым муравьем-агентом. После построения маршрута муравей-агент выполняет поиск в хэш-таблице; если маршрут обнаруживается, это означает, что другие агенты его уже рассмотрели, и текущему агенту требуется найти альтернативный путь. Альтернативой может служить упрощенный алгоритм ACOCNI, в котором при обнаружении маршрута в хэш-таблице муравей-агент не осуществляет обновление феромонного следа. В ходе исследования рассматривались стандартные сложные функции, обладающие множеством экстремумов и, как правило, зависящие от небольшого количества вещественных параметров. Переход к дискретной задаче осуществляется путем дискретизации параметров с заданной точностью, что приводит к большому количеству близких дискретных значений для одного параметра. Сократить количество значений для отдельного параметра возможно за счет декомпозиции его области определения на подслои с последующей линейной сверткой значений из этих подслоев [23]. При такой постановке задачи фактически осуществляется переход от функции с двумя параметрами к функции со многими параметрами: количество параметров увеличивается, а количество значений для каждого параметра уменьшается. Данный подход позволяет эффективно организовать параллельные вычисления для каждого отдельного слоя (нового параметра).

Предложенная в работе матричная формализация ACO предназначена для применения на SIMT-ускорителях, в которых для каждого потока выполняется одна и та же инструкция над различными данными, представленными, как правило, в виде векторов и матриц. Указанная особенность архитектуры позволяет производить вычисления на более чем 1000 потоках без использования сложных механизмов синхронизации. Типичным примером таких вычислителей являются графические процессоры.

1. Постановка задачи и методология исследования

В методе муравьиных колоний муравьи-агенты, решающие задачу на графе, перемещаются между вершинами, руководствуясь его состоянием: весами на дугах или вершинах (в оригинальном алгоритме интерпретируемыми как феромонный след) и различными параметрами дуг и вершин (в оригинальном алгоритме используется длина дуги). Выбор дуги осуществляется в результате реализации дискретного распределения вероятностей переходов муравья-агента. Вероятности этого распределения определяются следующей системой уравнений:
\[\begin{equation} 
\tag{1}
\begin{array}{l}
P_{i,k}(t)=
\begin{cases}
\dfrac{\tau^{\alpha}_i(t) ( {1}/{\eta_i} )^{\beta}}{\sum\limits_{l\in J_{i,k}}{\tau^{\alpha}_l(t) ( {1}/{\eta_l} )^{\beta}}}, & j\in J_{i,k}, 
\\ 
\qquad \qquad 0,  & j\notin J_{i,k}; 
\end{cases}
\\ 
\tau_i(t+1)=\begin{cases} 
\lambda\tau_i(t)+\sum\limits^N_{k=1}{\dfrac{Q}{L_k(t)}}, & \text{если } L_k(t)\to\min, \\ 
\lambda\tau_i(t)+\sum\limits^N_{k=1}{Q L_k(t)}, & \text{если } L_k(t)\to\max, \end{cases}
\end{array}
\end{equation} \]
где $P_{i,k}(t)$ — вероятность перехода по дуге $i$ для муравья-агента $k$ на итерации $t$ (при условии нахождения агента в определенной вершине, где $i$ — индексы дуг, по которым возможно перемещение, образующие множество $J_{i,k}$); $\tau^{\alpha}_i(t)$ — значение, ассоциированное с $i$-той дугой в динамическом слое $T$ на итерации $t$; $ ({1}/{\eta_i})^{\beta}$ — значение, ассоциированное с $i$-той дугой в статическом слое ${H}$ (величина, обратная длине дуги); $L_k(t)$ — значение целевой функции для муравья-агента $k$ на итерации $t$; $\alpha$, $\beta$ — параметры мультипликативной свертки; $\lambda$ — коэффициент испарения феромона, позволяющий уменьшать влияние ранее найденных решений; $Q$ — параметр добавления феромона (вследствие нормализации в формуле для вероятности не оказывает влияния на работу алгоритма, рекомендуется устанавливать $Q=1$). Формула (1) представляет собой мультипликативную свертку длины дуги и феромонного веса, взятых с весовыми коэффициентами $\alpha$ и $\beta$. Введем результирующий вес $z_i(t)=\tau^{\alpha}_i(t)({1}/{\eta_i})^{\beta}$. Для вычисления вероятности с учетом $z_i(t)\in {\mathbb R}$ необходима нормализация: 
\[
P_{i,k}(t)=\dfrac{z_i(t)}{\sum\limits_{l\in J_{i,k}}{z_l(t)}}, \quad j\in J_{i,k},
\]
обеспечивающая $\forall k,t$ выполнение условия $\sum_i{P_{i,k}(t)}=1$.

В данной работе предлагается представление информации, требуемой для вычислений по вероятностной формуле (1), в виде набора слоев (рис. 2). Слои подразделяются на статические, не изменяющиеся между итерациями алгоритма, и динамические. Информация о длинах дуг ${\eta}_i$ образует статический слой значений ${H}\in {\mathbb R}$; величина ${\eta}_i$ не зависит от номера итерации $t$. На основе статического слоя задачу можно решать различными алгоритмами; для задачи коммивояжера, например, могут применяться методы ветвей и границ или динамическое программирование. Уникальным для метода муравьиных колоний является динамический слой ${T}\in {\mathbb R}$, в котором определяются и обновляются феромонные веса ${\tau}_i(t)$ на каждой дуге графа, где $t$ — номер итерации. На начальной итерации алгоритма значения динамического слоя, как правило, идентичны, поэтому данный слой не влияет на вероятность выбора агентом дуги. По мере работы алгоритма значения в динамическом слое изменяются, и его роль в выборе дуги агентом становится определяющей. Представление информации в методе муравьиных колоний в виде системы слоев позволяет вводить дополнительные слои, например, слой посещения вершин или слой дополнительных статистических параметров.

 

Рис. 2. Разделение информации о графе на слои в методе муравьиных колоний
[Figure 2. Graph information layer partitioning in the ant colony optimization method]

 

В предлагаемой матричной модификации ACO каждый слой представлен матрицей, а взаимодействие между слоями осуществляется посредством простейших матричных операций. Данный подход позволяет эффективно организовать матричные и параллельные вычисления на SIMT-ускорителях. Применение матричной модификации ACO ориентировано на класс дискретных параметрических задач, которые традиционно решаются последовательными версиями ACO. Дискретная параметрическая задача задается в виде множества параметров $P$. Если в результате дискретизации вещественного параметра с заданной точностью производится его разбиение на подслои, то каждый такой подслой трактуется как отдельный параметр, а общее количество параметров обозначается величиной $n$. Для каждого параметра $p_i$ из множества $P={\{p_1,p_2,\dots, p_i,\dots, p_n\}}$ определяется множество его допустимых значений $V_i=\{v_{1,i},v_{2,i},\dots, v_{j,i},\dots, v_{m_i,i}\}$, $i=\overline{1,n}$. Множество $ V_i$ является дискретным в силу специфики метода муравьиных колоний. Количество допустимых значений $m_i>0$ для $i$-того параметра зависит от его природы. 
Таким образом, для различных параметров может быть установлено различное количество возможных значений. В процессе оптимизации формируется вектор значений параметров, интерпретируемый как решение: 
\[
X_k=(x_{1,k}, x_{2,k}, \dots, x_{i,k}, \dots, x_{n,k}),\quad 
|X_k|=n,
\] 
где $\forall i\ \exists j : {(x}_{i,k}=v_{j,i})\wedge (v_{j,i}\in  V_i)$. Здесь $k$ — индекс муравья-агента, который формирует решение независимо от других агентов. На каждой итерации алгоритма муравьиных колоний путь вычисляется одновременно для $K$ агентов, что позволяет параллельно получать каждое решение $X_k$. Сформированный вектор значений параметров передается на вычислительный модуль, который возвращает значение целевой функции $f(X_k)$. Решения рассматриваются в дискретном пространстве значений параметров, при этом требования непрерывности и дифференцируемости функции $f(X_k)$ не накладываются. Целью решения параметрической задачи является нахождение оптимального дискретного вектора $X^{*}$, обеспечивающего минимум (или максимум) целевой функции: $f(X_k)\to\min$.

Для матричного представления параметрической задачи необходимо представить все множества допустимых значений $V_i$, $i=\overline{1,n}$, в единой матричной форме. С этой целью вычисляется максимальное количество значений среди всех параметров: $m=\max\limits_{i} m_i$. Затем для всех слоев $i$, у которых $m_i<m$, производится дополнение фиктивными значениями до получения матрицы $V$ размерности $(n\times m)$. Решения, найденные всеми муравьями-агентами на итерации, описываются матрицей $X$ размерности $(K\times n)$, а в результате вычисления целевой функции определяется вектор значений $Y=f(X)$. Среди компонент вектора $Y$ на каждой итерации определяется минимальное значение; по завершении всех $S$ итераций фиксируется найденный глобальный минимум целевой функции. Количество итераций $S$ задается как параметр алгоритма ACO.

2. Модификация метода муравьиных колоний для решения параметрической задачи на CPU

Ввиду отсутствия в параметрической задаче информации о статическом слое $H$ метод муравьиных колоний на ранних этапах демонстрирует склонность к стагнации вблизи первого найденного рационального решения. Для предотвращения преждевременной стагнации и улучшения стратегии эксплуатации предложено дополнение вероятностной формулы информацией, связанной с количеством посещений вершин графа агентами, что формирует дополнительный слой $\theta$:
\[\begin{equation}
\tag{2}
z_{i,j}(t)=\lambda_1\tau^{\alpha}_{{\rm norm},i,j}(t)+
\lambda_2(1/\eta_{i,j}(t))^{\beta}+
\lambda_3(\eta_{i,j}(t)/\eta_{\max,i})^{\gamma},
\end{equation}\]
где $i$ — номер параметра $p_i\in P $; $j$ — номер значения параметра $v_{j,i}\in V_i$; $t$ — номер итерации; $\lambda_1$, $\lambda_2$, $\lambda_3\in {\mathbb R}$ — коэффициенты аддитивной свертки, $\lambda_1$, $\lambda_2$, $\lambda_3\in[0,\infty)$; $\alpha$, $\beta$, $\gamma \in {\mathbb R}$ — показатели степени при слагаемых аддитивной свертки (учтены в связи с их наличием в (1)). Предполагается, что при применении модификации в реальных задачах будут использоваться только коэффициенты $\lambda_1$, $\lambda_2$, $\lambda_3$. Вычисление нормализованного значения феромонного веса для использования в аддитивной свертке выполняется по формуле
\[
\tau_{{\rm norm},i,j}(t)=\tau_{i,j}(t)/\sum^{m_i}_{l=1}{\tau_{i,l}(t)},
\]
где $\tau_{i,j}(t)\in T$ — нормированное значение динамического слоя $T$ для значения с номером $j$ параметра с номером $i$ на итерации $t$; $\eta_{i,j}(t)$ — количество использований в рассмотренных решениях значения с номером $j$ параметра с номером $i$ на итерации $t$ (слой $\theta$); $\eta_{\max,i}$ — максимально возможное количество использований значений для параметра с номером $i$.

3. Модификация метода муравьиных колоний для выполнения вычислений 
на SIMT-ускорителях

Для реализации вычислений на SIMT-ускорителях предлагается использовать матричную формализацию, предполагающую возможность параллельного выполнения операций. В рамках данного подхода вся информация в слоях хранится в матричной форме: матрица феромонных весов $T$  размерности $(n\times m)$ и матрица количества посещений вершин $\theta$  размерности $(n\times m)$. При приведении векторов значений параметров $ V_i$ к единой матрице $ V$ размерности $(n\times m)$ добавляются фиктивные значения параметров, которые могут быть произвольными, поскольку вероятность их выбора в методе муравьиных колоний будет нулевой из-за нулевого начального феромонного веса. Для обеспечения этого требования начальные значения матриц задаются следующим образом:
\[
\begin{array}{l}
\tau_{i,j}(0)=
\begin{cases}
1, & \text{если } j \leqslant m_i; \\
0, & \text{если } j > m_i; 
\end{cases}
\\
  \eta_{i,j}(0)=1,
\end{array}
\]
где $i$, $j$ — индексы строки и столбца матрицы $V$ соответственно. Отдельно следует отметить, что в данной постановке феромонные веса наносятся на вершины, а не на дуги параметрического графа.

На первом этапе работы алгоритма вычисляется матрица вероятностей выбора вершин. Поскольку сумма вероятностей должна удовлетворять условию нормировки $\sum_j{P_{i,j}(t)}=1$  $\forall i$, $t$, для определения матрицы вероятностей $P$ размерности $(n\times m)$ необходимо вычислить вектор-строку:
\[\begin{gather*}
Z'_i=\sum^m_{j=1} Z_{i,j},\quad |Z'|=n; 
\\
P= Z / Z',  \text{ где } Z=\lambda_1  T^{\alpha }+\lambda_2(1/\theta )^{\beta }+\lambda_3(\theta  / \theta ')^{\gamma },\quad Z_{ (n\times m)},
\end{gather*}\]
где ${\mathbf \theta }'$ — вектор максимально возможного количества посещений вершин, который может быть определен априори как $\theta '_i={\prod^n_{k=1}{m_k}}/{m_i}$, $|\theta '|=n$. Следует подчеркнуть необходимость нормализации матрицы $T$: $\tau_{{\rm norm},i,j} = \tau_{i,j} / \sum_{l=1}^{m} \tau_{i,l}$, чтобы обеспечить условие $\tau_{{\rm norm},i,j}\in [{0, 1}]$ для их использования в аддитивной свертке.

Выбор значений параметров, формирующих решение (вектор $X_k$), осуществляется на основе распределения вероятностей выбора вершин для каждого параметра (строки матрицы $P$) с применением метода обратной функции. Для работы этого метода требуется вычисление матрицы функций распределения $ F$ размерности $(n\times m)$, где $F_{i,j}=\sum^j_{l=1} P_{i,l}$. На первом этапе все матрицы имеют размерность $(n\times m)$, а векторы — размерность $n$.

На втором этапе определяется путь для каждого муравья-агента в группе и вычисляется значение целевой функции. Количество муравьев-агентов в группе задается параметром $K$. Для формирования пути каждого агента генерируется матрица случайных величин $R$ размерности $(K\times n)$, равномерно распределенных на интервале $(0, 1)$. Индекс $s$ для метода обратной функции определяется для каждого сгенерированного числа таким образом, чтобы выполнялось неравенство
\[\begin{equation} 
\tag{3}
F_{i,s-1} \leqslant r_{k,i} < F_{i,s} \quad  \forall i.
\end{equation}\]

Матрица решений (путей муравьев-агентов в параметрическом графе) $X$ размерности $(K\times n)$ определяется на основе матрицы $V$: $x_{k,i}=v_{i,s}$. Для каждого полученного решения вычисляется значение целевой функции $f(X_k)$, формирующее вектор $Y$, $|Y|=K$. В модификациях ACOCCyI и ACOCNI перед вычислением целевой функции выполняется проверка наличия решения в хэш-таблице. При отсутствии коллизий или их низкой вероятности поиск в хэш-таблице осуществляется за время $O(1)$. 
Если решение отсутствует, значение целевой функции вычисляется и в таблицу добавляется новая запись.

На третьем этапе выполняется обновление динамических слоев. После получения всех значений целевой функции производится обновление феромонных весов: 
\[
T(t+1)=\rho T(t)+\Delta T,
\]
где $\Delta T$ вычисляется на основе вектора $Y$. Разрешение несоответствия размерностей между  $T_{(n\times m)}$ и $Y_{(K)}$ осуществляется путем добавления элемента вектора $Y$ только к тем вершинам (значениям параметров), которые были выбраны муравьем-агентом, в соответствии с матрицей $X_{(K\times n)}$. Алгоритмически данная процедура может быть записана в виде 
\[
T_{x_{k,j}, j}(t+1) = T_{x_{k,j}, j}(t+1) + Q / Y_k
\]
для $j=\overline{1,n}$ и $k=\overline{1,K}$.

4. Реализация метода муравьиных колоний для работы на SIMT-ускорителях

Алгоритмическое обеспечение предложенной модификации метода муравьиных колоний базируется на архитектурных особенностях SIMT-ускорителей. В отличие от потоковых вычислителей, SIMT-ускоритель выполняет одну инструкцию над различными наборами данных. Для тестирования предложенного алгоритма была разработана программная реализация для графических ядер CUDA (NVIDIA), поэтому терминология параллельных решений опирается на концепции CUDA: загрузка данных и запуск функций осуществляются с центрального процессора (CPU), тогда как параллельные вычисления на GPU выполняются потоками (thread), объединенными в блоки (block).

После инициализации всех матриц и векторов в памяти GPU запускается процедура первого этапа. Поскольку на первом этапе все матрицы имеют размерность $(n\times m)$, а векторы — размерность $n$, все операции алгоритма выполняются параллельно на $n$ потоках. Каждый поток обрабатывает данные для одного параметра (или для одного подслоя, если параметр был разделен). 
Для каждого потока (с идентификатором $tx$) выполняется следующая последовательность действий.

1. Формируется нормализованный вектор феромонных весов $T_{{\rm norm},tx,j}$.  Общее количество арифметических операций составляет $2m$: $m$ операций для вычисления суммы $\sum_{l=1}^{m} \tau_{tx,l}$ и $m$ операций для выполнения поэлементного деления $\tau_{{\rm norm},tx,j} = \tau_{tx,j} / \sum_{l=1}^{m} \tau_{tx,l}$.
2. Вычисляется вектор $Z_{tx,j}$. Для этого выполняется $m$ операций по вычислению аддитивной свертки (2) для каждого $j$ с одновременным накоплением суммы ${Z'}_{tx}  = \sum_{j=1}^{m} Z_{tx,j}$.
3. Определяется вектор функции распределения $F_{tx,j}$. Для его расчета требуется $m$ операций для вычисления накопленной суммы $F_{tx,j} = \sum_{l=1}^{j} P_{tx,l}$, где $P_{tx,l} = Z_{tx,l} / {Z'}_{tx}$.

При разделении числовых параметров на отдельные подслои стандартный размер максимального слоя составляет $m=5$. Соответственно, на первом этапе в каждом потоке выполняется приблизительно $5 \cdot 4 + 4 = 24$ операций. По завершении первого этапа осуществляется синхронизация потоков и управление возвращается на CPU.

На втором этапе каждый блок $bx$ выполняет вычисления для одного муравья-агента $k$, а в каждом потоке блока $tx$ производится выбор значения $x_{k,tx}$ для соответствующего параметра. Таким образом, в потоке выполняются операции по вычислению значения
отдельного параметра.

1. Генерируется случайная величина $r_{bx,tx}$, равномерно распределенная на интервале $(0, 1)$.
2. В цикле определяется индекс $s_{bx,tx}$ для метода обратной функции. Максимальное количество итераций цикла равно $m$, однако возможен досрочный выход при выполнении условия (3). Для SIMT-вычислений досрочный выход из цикла не рекомендуется, поскольку он нарушает синхронность выполнения инструкций (если один поток завершит цикл на второй итерации, а другому потребуется четыре).
3. Выбранное значение параметра фиксируется: $x_{bx,tx}=v_{tx, s_{bx,tx}}$; при этом $x_{bx,tx}$ является элементом вектора решения $X_{bx}$.
4. Поскольку все предыдущие операции в блоке выполняются синхронно, после завершения пункта 3 становятся известны все выбранные значения параметров (от каждого потока $tx$) для каждого агента (блока $bx$). Это позволяет проверить наличие решения $X_{bx}$ в хэш-таблице.
5. Если решение $X_{bx}$ отсутствует в хэш-таблице, то вычисляется значение целевой функции $Y_{bx}=f(X_{bx})$. В противном случае, в зависимости от алгоритма, устанавливается $Y_{bx}=0$ (для ACOCNI) или выполняется возврат к пункту 1 для данного агента (для ACOCCyI).

Для выполнения второго этапа (при $m=5$) требуется примерно $1+5+1+1+1=9$ операций на поток. Следует отметить, что алгоритм ACOCCyI нарушает принципы SIMT-вычислений, поскольку для различных агентов может потребоваться разное количество операций. 
Указанную проблему можно устранить, используя алгоритм ACOCNI либо организуя дополнительный циклический поиск на последующих итерациях.

На третьем этапе каждый поток $tx$ выполняет обновление матриц $T$ и $\theta$ для закрепленного за ним параметра.

1. Выполняется «испарение» феромонных весов для всех $m$ вершин параметра: $T'_{tx,j} = \rho \cdot T_{tx,j}$.
2. Выполняется цикл по всем муравьям-агентам ($k = \overline{1,K}$). Каждый агент для параметра $tx$ выбрал одно значение. Определяется индекс выбранного значения $s = X_{k,tx}$, после чего обновляется соответствующий феромонный вес: 
   \[
   T''_{tx,s} = T'_{tx,s} + Q / Y_k.
   \]

Вычислительная сложность третьего этапа зависит от параметра алгоритма — количества агентов $K$ на итерации.

После завершения каждого этапа выполняется синхронизация всех потоков и управление возвращается с GPU на CPU. Данная операция вносит существенные накладные расходы и замедляет выполнение алгоритма. Для оптимизации быстродействия при реализации на CUDA целесообразно объединить этапы 3 и 1, что позволяет сократить количество переключений контекста до двух.

5. Обсуждение результатов работы и область их применения

Для тестирования предложенного алгоритма было реализовано программное обеспечение¹ на языке C/C++ на SIMD-вычислителе с использованием CUDA Toolkit.² Тестирование проводилось на задачах оптимизации, связанных с поиском минимума/максимума различных тестовых функций [24], например функции Шаффера:
\[
f(x_1,x_2)=\frac{1}{2}+\frac{\sin^2 (\sqrt{x_1^2+x_2^2} )-0.5}{1+0.001(x_1^2+x_2^2)}
\]
на интервале $[{-10; 10}]$ с точностью до $10^{-10}$. Каждый параметр был дискретизирован на 20 подслоев, что позволило сформировать матрицу $ V$ размерности $(40\times 5)$.

В табл. 1 представлены оценки математического ожидания времени выполнения различных этапов метода муравьиных колоний для модификации ACOCNI. Количество итераций алгоритма составило 500, для сбора статистики было выполнено 300 прогонов модели. Время выполнения этапа 2, а также время выделения памяти и копирования данных между CPU и GPU зависят от количества агентов на итерации, поскольку количество параллельных потоков определяется числом агентов. Результаты получены на графическом процессоре NVIDIA GeForce RTX 3060, содержащем 3840 унифицированных ядер.

 

Таблица 1. Оценка математического ожидания времени (в мс) выполнения различных этапов для модификации ACOCNI при изменении количества муравьев-агентов [Estimation of the mathematical expectation of execution time (in ms) for various stages of the ACOCNI modification when varying the number of ant agents]

Number of ants	Memory allocation	Synchronization and copying of data	Stage 1	Stage 2	Stage 3
1000	8.47	121.31	34.03	261.02	696.99
2000	8.44	195.33	35.39	464.04	1438.67
3000	7.63	257.49	36.16	676.24	2183.71
4000	8.10	311.39	37.80	872.19	2925.37
5000	8.26	356.13	37.54	1065.31	3661.39

 

Для алгоритма ACOCCyI результаты измерения времени приведены в табл. 2. Согласно полученным данным, при увеличении количества муравьев-агентов на итерации время выполнения этапа 2 резко возрастает, поскольку именно на этом этапе осуществляется циклический поиск нового решения.

 

Таблица 2. Оценка математического ожидания времени (в мс) выполнения различных этапов для модификации ACOCCyI при изменении количества муравьев-агентов [Estimation of the mathematical expectation of execution time (in ms) for various stages of the ACOCCyI modification when varying the number of ant agents]

Number of ants	Memory allocation	Synchronization and copying of data	Stage 1	Stage 2	Stage 3
1000	7.73	210.88	38.38	3530.61	721.12
2000	8.52	291.86	42.00	6667.19	1460.65
3000	7.99	319.20	38.92	9554.89	2200.23
4000	7.75	610.48	34.49	12324.47	2585.48
5000	8.20	342.28	39.34	19642.84	3698.99

 

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

 

Таблица 3. Оценка математического ожидания времени (в мс) выполнения различных этапов для модификации ACOCNI при изменении количества итераций [Estimation of the mathematical expectation of execution time (in ms) for various stages of the ACOCNI modification when varying the number of iterations]

Number of ants	Memory allocation	Synchronization and copying of data	Stage 1	Stage 2	Stage 3
1000	8.45	292.73	34.47	1040.51	0.77
1500	8.12	577.30	69.01	1948.56	1.55
2000	8.62	1089.04	107.92	2777.23	2.52
2500	8.52	1315.47	140.59	3567.15	3.14
3000	10.35	2010.78	180.08	4292.08	5.87

 

Полученные результаты демонстрируют высокую эффективность параллельного мелкозернистого варианта метода муравьиных колоний, формализованного в терминах матричных вычислений и адаптированного для выполнения на современных GPU. Использование SIMT-архитектуры позволяет существенно сократить время вычислений за счет ослабления требований к синхронизации потоков, что, в свою очередь, ускоряет работу алгоритма в целом. Активное развитие кластерных систем и суперЭВМ делает доступными не только вычисления на GPU, но и облачные вычисления с применением гетерогенных архитектур. В таких условиях параметрическая оптимизация методом муравьиных колоний позволяет автоматизировать процесс подбора гиперпараметров систем, функционирующих на вычислительных кластерах.

В табл. 4 приведены оценки времени работы многопоточной (GPU) и однопоточной (CPU) версий предложенного метода. Наблюдается, что использование графического процессора обеспечивает существенное ускорение вычислений.

Необходимо отметить, что матричная реализация метода муравьиных колоний демонстрирует высокую скорость работы и эффективность поиска решения даже на однопоточном процессоре, что достигается за счет оптимизации компилятором  C++ вычислительных процедур, в особенности на этапе 3. Кроме того, реализация алгоритма основана на работе с непрерывными областями памяти, где двумерные массивы представлены с помощью адресной арифметики. Благодаря эффективному кэшированию данных на CPU также достигается увеличение скорости выполнения этапа 3.

 

Таблица 4. Оценка времени выполнения (в мс) однопоточной реализации и реализации с применением CUDA [Estimated execution time (in ms) of single-threaded implementation and CUDA-based implementation]

	Time in ms ($10^3$ agents, $10^3$ iterations)
	ACOCNI	ACOCNI  + CUDA	ACOCCyI	ACOCCyI + CUDA
Memory allocation	0.07	8.47	0.06	7.73
Synchronization and copying of data	0.87	121.31	65.10	210.88
Stage 1	1.10	34.03	1.31	38.38
Stage 2	817.58	261.02	65039.84	3530.61
Stage 3	54.62	696.99	62.66	721.12
Result	873.30	992.04	65103.82	4290.11
CUDA acceleration	88.03%		1517.53%

 

Заключение

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

При этом алгоритм ACOCCyI в целом является предпочтительнее по сравнению с алгоритмом ACOCNI, поскольку быстрее находит оптимальное решение в параметрическом графе. Алгоритм ACOCNI, в свою очередь, демонстрирует лучшие временные показатели: 
на первом этапе требуется выполнение приблизительно 24 операций, на втором — 9 операций, а на третьем — $5+N$ операций. 
Алгоритм ACOCCyI может потребовать большего количества операций на втором этапе, что сильно зависит как от количества агентов на итерации, так и от сложности оптимизируемой функции и скорости сходимости алгоритма.

Дальнейшие исследования планируется направить на разделение участков алгоритма между вычислителями типа CPU или MIMD (Multiple Instruction, Multiple Data) и вычислителями типа GPU или SIMD. Такое разделение с последующим анализом времени выполнения отдельных блоков позволит спроектировать оптимальную структуру гетерогенного вычислителя, сочетающего SIMT- и MIMD-процессоры с высокоскоростными механизмами передачи данных. В настоящее время параметрическая оптимизация гиперпараметров методом муравьиных колоний применяется для задачи подбора оптимальных значений параметров моделей SARIMAX и LSTM-сетей, используемых для прогнозирования объемов пассажирских и грузовых авиаперевозок.

Конкурирующие интересы. Авторы заявляют об отсутствии конфликта интересов, который мог бы повлиять на объективность представленных результатов или процедуру рецензирования данной статьи.
Авторский вклад и ответственность. Все авторы внесли равный вклад в разработку концепции исследования, проведение вычислительных экспериментов, анализ полученных результатов и подготовку рукописи. Все соавторы несут полную ответственность за содержание публикации и одобрили финальную версию статьи перед подачей в печать.
Финансирование. Исследование выполнено при финансовой поддержке АНО «Аналитический центр при Правительстве Российской Федерации» в рамках договора № 70-2025-000814 от 04.06.2025.

---

¹https://github.com/kalengul/ACO_SIMD
²https://developer.nvidia.com/cuda-toolkit

About the authors

Vladimir A. Sudakov

Moscow Aviation Institute (National Research University); Keldysh Institute of Applied Mathematics of Russian Academy of Sciences

Author for correspondence.
Email: sudakov@ws-dss.com
ORCID iD: 0000-0002-1658-1941
SPIN-code: 1614-4760
Scopus Author ID: 7006296922
ResearcherId: I-2573-2018
https://www.mathnet.ru/rus/person112232

Dr. Techn. Sci., Associate Professor; Professor; Dept. of Computational Mathematics and Programming; Leading Researcher

Russian Federation, 125993, Moscow, Volokolamskoe Shosse, 4; 125047, Moscow, Miusskaya pl., 4

Yurii P. Titov

Moscow Aviation Institute (National Research University)

Email: kalengul@mail.ru
ORCID iD: 0000-0002-9093-6755
SPIN-code: 8866-5903
Scopus Author ID: 56549731400
ResearcherId: AAO-7748-2020
https://www.mathnet.ru/rus/person111602

Cand. Techn. Sci., Associate Professor; Associate Professor; Dept. of Computing Machines, Systems and Networks

Russian Federation, 125993, Moscow, Volokolamskoe Shosse, 4

References

Supplementary files