Журнал вычислительной математики и математической физики, 2021, T. 61, № 5, стр. 845-864

Моделирование структуры данных с помощью блочного тензорного разложения: разложение объединенных тензоров и вариационное блочное тензорное разложение как параметризованная модель смесей

И. В. Оселедец 12*, П. В. Харюк 123**

1 Сколковский институт науки и технологий
121205 Москва, Большой бульвар, 30, стр. 1, Россия

2 Институт вычислительной математики им. Г.И. Марчука РАН
119333 Москва, ул. Губкина, 8, Россия

3 МГУ им. М.В. Ломоносова, Факультет вычислительной математики и кибернетики
119991 Москва, Ленинские горы, 1, стр. 52, Россия

* E-mail: ivan.oseledets@gmail.com
** E-mail: kharyuk.pavel@gmail.com

Поступила в редакцию 24.12.2020
После доработки 24.12.2020
Принята к публикации 14.01.2021

Полный текст (PDF)

Аннотация

Развивается идея использования тензорных разложений в качестве параметрической модели группового анализа данных. Представлены две модели на основе блочного тензорного разложения: детерминистическая и вероятностная, с использованием различных форматов слагаемых. Установлена связь между блочным тензорным разложением и смесями непрерывных латентных вероятностных моделей: на основе блочного тензорного разложения построена модель смеси распределений со структурированным представлением. Модели были протестированы в задаче кластеризации набора цветных изображений и данных электрической активности мозга. Результаты показывают, что предложенные подходы способны к выделению релевантной индивидуальной составляющей данных. Библ. 54. Фиг. 4. Табл. 5.

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

1. ВВЕДЕНИЕ

В процессе поиска зависимостей в данных важно учитывать их природу. Ряд собранных наборов данных может быть представлен в тензорном виде (как многомерные массивы), и учет их размерности представляется важным не только для эффективной обработки, но и для потенциальной интерпретации найденных параметров в конкретном прикладном исследовании. Более того, подобные данные часто являются избыточными и могут быть приближены меньшим числом параметров. Тензорные разложения предоставляют возможность построения подобных приближений. Различные тензорные форматы, такие как канонический (полилинейный, CPD) или формат Таккера (см. [1], [2]), нашли успешное применение в задачах анализа данных (см. [3], [4]). В данной работе в качестве основы для моделей мы использовали так называемое блочное тензорное разложение (BTD), которое является одним из способов обобщения классических тензорных форматов и было предложено в [5]–[7].

В условиях ряда исследований требуется извлекать общую и индивидуальную информацию из данных. Например, подобные исследования возникают в когнитивных, медицинских и иных задачах изучения мозга по группе из нескольких человек (см., например, [8]–[10]); другим примером может служить химическая экспертиза набора образцов (см. [11]–[13]). Естественно ожидать от подобных наборов данных наличие общих (групповых) и индивидуальных частей, и реконструкцию соответствующих частей из данных мы будем называть групповым анализом данных. Предложенные в данной работе модели структуры данных основаны на BTD со слагаемыми различных типов, что позволяет моделировать групповые и индивидуальные части различными форматами.

Предложенные тензорные модели структуры данных основаны на концепте “связанного многомерного анализа компонент” (см. [14], [15]), реализованного в виде BTD с дополнительными условиями; кроме того, реализована модель смеси вероятностных распределений, параметризованная через BTD и названная вариационным блочным тензорным разложением (VBTD). В предыдущих работах по групповому анализу данных (см. [16], [17]) специальный случай первой (детерминистической) модели (а именно, условное BTD со слагаемыми в (Lr, 1) формате) был рассмотрен как составляющий элемент обработки данных для задачи классификации. В настоящей работе модели были протестированы в условиях задачи кластеризации. Сведения об использованных в работе обозначениях приведены в п. 10.1.

2. СТРУКТУРИРОВАНИЕ НАБОРОВ ДАННЫХ ДЛЯ ГРУППОВОГО АНАЛИЗА

Одна из главных целей группового анализа данных может быть сформулирована как выделение общей информации из комбинированного набора данных. Вопрос в том, что может быть названо “общей информацией”. В данной работе мы придерживаемся распространенной в области обработки сигналов точке зрения на природу данных: предположения о том, что наблюдения $\{ {{X}_{i}}\} _{{i = 1}}^{N}$ генерируются в результате преобразования активности нескольких источников, представленных в векторном виде. Адаптация такой точки зрения к групповому анализу данных ведет к разделению источников на общие, ${{S}_{{{\text{com}}}}}$, и индивидуальные, $\{ {{S}_{i}}\} _{{i = 1}}^{N}$.

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

Предположение группового анализа данных. Наблюдения $\{ {{X}_{i}}\} _{{i = 1}}^{N}$ могут быть приближены аддитивной функциональной зависимостью от общих, ${{S}_{{{\text{com}}}}}$, и индивидуальных, $\{ {{S}_{i}}\} _{{i = 1}}^{N}$, источников

(1)
${{X}_{i}} \approx {{F}_{i}}({{S}_{{{\text{com}}}}},{{S}_{i}}) = {{f}_{i}}({{S}_{{{\text{com}}}}}) + {{g}_{i}}({{S}_{i}}),\quad i = \overline {1,N} .$

В случае, если наблюдения ${{X}_{i}}$ двумерны (являются матрицами), распространенным способом моделирования указанной зависимости является линейное смешивание:

(2)
${{f}_{i}}({{S}_{{{\text{com}}}}}) = {{S}_{{{\text{com}}}}}B_{{{{f}_{i}}}}^{{\text{T}}},\quad {{g}_{i}}({{S}_{i}}) = {{S}_{i}}B_{{{{g}_{i}}}}^{{\text{T}}},\quad i = \overline {1,N} ,$
где ${{B}_{{{{f}_{i}}}}}$, ${{B}_{{{{g}_{i}}}}}$ – матрицы коэффициентов общих и индивидуальных источников соответственно для образца ${{X}_{i}}$. Модель факторизации данных ставит вопрос о том, как найти требуемые параметры разложения. В [18] предложен ответ для общей части в виде группового метода независимых компонент (group ICA). Zhou и др. (см. [19]) построили алгоритм вычисления обеих частей разложения через приближение матрицы общих источников, ортогональной всем индивидуальным источникам. В [20] общая часть оценивается через малоранговое разложение склеенных матричных данных $X = [{{X}_{1}} \ldots {{X}_{N}}] \approx {{S}_{{{\text{com}}}}}\mathop {[{{B}_{{{{f}_{1}}}}} \ldots {{B}_{{{{f}_{N}}}}}]}\nolimits^{\text{T}} $ и каждая индивидуальная часть вычисляется как малоранговое разложение данных после удаления групповой части: ${{X}_{i}} - {{S}_{{{\text{com}}}}}B_{{{{f}_{i}}}}^{{\text{T}}} \approx {{S}_{i}}B_{{{{g}_{i}}}}^{{\text{T}}}$.

В случае многомерных данных ${{X}_{i}}$ размерами ${{n}_{1}} \times \ldots \times {{n}_{d}}$, как правило, используются предварительно матризованные наблюдения. Однако в этом случае рассматриваемая модель игнорирует размерность входных данных. Детерминистическая BTD модель, предложенная в данной работе, обобщает подход для больших размерностей следующим образом. Рассмотрим матрицу развертки тензорных данных по выбранной моде источников $k \in \{ 1,2, \ldots ,d\} $ (например, отвечающей пикселям или времени) и применим к ней сформулированные ранее предположения (1), (2):

(3)
${\text{unfol}}{{{\text{d}}}_{k}}\left( {{{X}_{i}}} \right) = \mathop {\left( {{{X}_{i}}} \right)}\nolimits_{(k)} \approx {{S}_{{{\text{com}}}}}B_{{{{f}_{i}}}}^{{\text{T}}} + {{S}_{i}}B_{{{{g}_{i}}}}^{{\text{T}}},\quad i = \overline {1,N} .$
Далее учтем тензорную структуру слагаемых в правой части (3). Также вместо вычисления набора связанных разложений объединим данные $\{ {{X}_{i}}\} _{{i = 1}}^{N}$ вдоль новой групповой оси, получив общий тензор данных $X$ размерами ${{n}_{1}} \times \ldots \times {{n}_{d}} \times N$, и наложим условия на специальную структуру параметров, отвечающих групповой моде. Объединение данных требует выполнения условий однородности для $\{ {{X}_{i}}\} _{{i = 1}}^{N}$: все образцы должны иметь равную размерность с одинаковыми размерами мод.

3. ДЕТЕРМИНИСТИЧЕСКАЯ ГРУППОВАЯ МОДЕЛЬ НА ОСНОВЕ BTD РАЗЛОЖЕНИЯ ОБЪЕДИНЕННЫХ ДАННЫХ

После учета предположений из разд. 2 модельная структура входных данных примет следующий вид:

(4)
$X \approx \sum\limits_{i = 1}^N {\left( {{{X}_{i}}({{S}_{{{\text{com}}}}},{{\theta }_{{{\text{com}},i}}}) + {{X}_{i}}({{S}_{i}},{{\theta }_{{{\text{ind}},i}}})} \right)} \circ {{{\mathbf{e}}}_{i}},$
где ${{{\mathbf{e}}}_{i}}$ есть $i$-й столбец единичной матрицы ${{I}_{N}}$, $ \circ $ обозначает внешнее произведение, ${{\theta }_{{{\text{com}},i}}}$, ${{\theta }_{{{\text{ind}},i}}}$ – общие и индивидуальные параметры, соответствующие тому или иному тензорному формату и отвечающие общим, ${{S}_{{{\text{com}}}}}$, и индивидуальным, ${{S}_{i}}$, источникам. Сосредоточимся на двух тензорных форматах: формате Таккера и (Lr, 1) (специальном каноническом с параметром $P$, отвечающим числу мод с полноразмерными фактор-матрицами), который был предложен в [21]. В случае (Lr, 1) формата для ${{T}_{i}} = {{X}_{i}}({{S}_{i}},{{\theta }_{{{\text{ind}},i}}})$ соответствующая сумма может быть переписана в следующем виде:
(5)
$\sum\limits_{i = 1}^N \,{{T}_{i}} \circ {{{\mathbf{e}}}_{i}} = \sum\limits_{i = 1}^N \,\left| {\left[ {C_{1}^{{[i]}}, \ldots ,C_{P}^{{[i]}},{\mathbf{c}}_{{P + 1}}^{{[i]}}{{E}^{{[i]}}}, \ldots ,{\mathbf{c}}_{d}^{{[i]}}{{E}^{{[i]}}}} \right]} \right| \circ {{{\mathbf{e}}}_{i}} = \left| {\left[ {{{C}_{1}}, \ldots ,{{C}_{d}}E,{{I}_{N}}E} \right]} \right|,$
где умножение фактор-матриц $\{ {{C}_{l}}\} _{{l = P + 1}}^{d}$ на $E = [{{E}^{{[1]}}}, \ldots ,{{E}^{{[N]}}}]$ дублирует необходимое количество раз их столбцы. Выражение справедливо для обеих частей разложения: индивидуальных и общей. В последнем случае дополнительно предполагаем, что ${{X}_{i}}({{S}_{{{\text{com}}}}},{{\theta }_{{{\text{com}},i}}}) = {{p}_{i}}{{T}_{G}}$, с дополнительными параметрами $\sum\nolimits_{i = 1}^N \,{{p}_{i}} = 1$, ${{p}_{i}} \geqslant 0$, которые позволяют регулировать присутствие общей части ${{T}_{G}}$. Из этого предположения следует, что групповые части модели могут быть представлены как ${{T}_{G}} \circ {{\boxed1}_{{N \times 1}}}$, где ${{\boxed1}_{{N \times 1}}}$ – вектор из единиц размерами $N \times 1$. Сводя все предположения воедино, а также задав число компонент (ранги) для индивидуальных и общих слагаемых, $L = [{{L}_{1}} \ldots {{L}_{N}}{{L}_{{N + 1}}}]$, приходим к следующей модели:
(6)
$\begin{gathered} \hfill \\ X \approx \left| {[{{C}_{1}}, \ldots ,{{C}_{d}},{{C}_{{d + 1}}}]} \right|,\quad {{C}_{k}} = \left\{ \begin{gathered} [\overbrace {C_{k}^{{[1]}}}^{{{L}_{1}}}, \ldots ,\overbrace {C_{k}^{{[N + 1]}}}^{{{L}_{{N + 1}}}}],\quad k = \overline {1,P} , \hfill \\ [{\mathbf{c}}_{j}^{{[1]}}, \ldots ,{\mathbf{c}}_{j}^{{[N + 1]}}]E,\quad k = \overline {P + 1,d} ,\quad j = k - P, \hfill \\ [{{I}_{N}}{\mathbf{p}}]E,\quad k = d + 1, \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} $
где $E = \left[ {{{{\boxed1}}_{{1 \times {{L}_{1}}}}} \otimes {{{\mathbf{e}}}_{1}}, \ldots ,{{{\boxed1}}_{{1 \times {{L}_{{N + 1}}}}}} \otimes {{{\mathbf{e}}}_{{N + 1}}}} \right]$, $ \otimes $ – произведение Кронекера, p – вектор параметров, определенный ранее. Однако вариативность групповой части в данной модели может быть слишком ограниченной. Для моделирования более сложной зависимости от общей части можно использовать формат Таккера со схожими условиями на фактор-матрицу групповой оси. Таким образом, полностью (Lr, 1) групповую модель (далее GLRO) дополним похожей моделью, в которой общая часть представлена в формате Таккера (далее GTLD):
(7)
$\begin{gathered} X \approx \left| {\left[ {G;{{A}_{1}}, \ldots ,{{A}_{d}},{{A}_{{d + 1}}}} \right]} \right| + \left| {\left[ {{{C}_{1}}, \ldots ,{{C}_{d}},{{C}_{{d + 1}}}} \right]} \right|, \\ \\ {{A}_{{d + 1}}} = \operatorname{diag} ({\mathbf{p}}),\quad {{C}_{k}} = \left\{ \begin{gathered} [\overbrace {C_{k}^{{[1]}}}^{{{L}_{1}}}, \ldots ,\overbrace {C_{k}^{{[N]}}}^{{{L}_{N}}}],\quad k = \overline {1,P} , \hfill \\ [{\mathbf{c}}_{j}^{{[1]}}, \ldots ,{\mathbf{c}}_{j}^{{[N]}}]E,\quad k = \overline {P + 1,d} ,\quad j = k - P, \hfill \\ {{I}_{N}},\quad k = d + 1, \hfill \\ \end{gathered} \right. \\ \end{gathered} $
где последняя мода отвечает групповой оси, как и в (6).

В силу того что как канонический и (Lr, 1) форматы, так и формат Таккера одинаковым образом разделяют оси тензора, имеет смысл указать общий случай BTD со слагаемыми во всех трех форматах:

(8)
$X \approx \underbrace {\sum\limits_{m = 1}^M {\left| {\left[ {{{G}^{{[m]}}};A_{1}^{{[m]}}, \ldots ,A_{{d + 1}}^{{[m]}}} \right]} \right|} }_{{\text{в}}\;{\text{формате}}\;{\text{Таккера}}} + \underbrace {\left| {\left[ {{{B}_{1}}, \ldots ,{{B}_{{d + 1}}}} \right]} \right|}_{{\text{канонический}}} + \underbrace {\left| {\left[ {{{C}_{1}}, \ldots ,{{C}_{{d + 1}}}} \right]} \right|}_{({\text{Lr}},1)}.$
Отметим, что программная реализация предложенных моделей позволяет вычислять параметры такого разложения более общего вида. В групповом анализе данных также можно моделировать индивидуальные части в формате Таккера, но это потребует увеличения числа гиперпараметров (рангов) для каждого отдельного слагаемого. В случае же (Lr, 1) формата потребуется указать лишь один вектор гиперпараметров $L$, который задает ранги всех полноразмерных фактор-матриц, отвечающих индивидуальным слагаемым.

Помимо гиперпараметров, упомянутых ранее ($P$, $L$, и/или ранги Таккера), предложенный подход требует выбора некоторого подмножества среди первых $P$ осей, $\Omega \subseteq \{ 1,2, \ldots ,P\} $, для которых вводится условие отделимости соответствующих общих и индивидуальных фактор-матриц (аналогично [19]) с тем, чтобы убрать пересечения между общей и индивидуальными активностями:

(9)
$\begin{gathered} \mathop {(C_{\gamma }^{{[N + 1]}})}\nolimits^{\text{T}} C_{\gamma }^{{[i]}} = 0,\quad {\text{для}}\;{\text{модели}}\;(6), \hfill \\ \mathop {({{A}_{\gamma }})}\nolimits^{\text{T}} C_{\gamma }^{{[i]}} = 0,\quad {\text{для}}\;{\text{модели}}\;(7),\quad \forall i = \overline {1,N} ,\quad \forall \gamma \in \Omega . \hfill \\ \end{gathered} $
В данной работе мы остановились на выборе множества $\Omega $, состоящего из одной моды.

4. ВАРИАЦИОННОЕ БЛОЧНОЕ ТЕНЗОРНОЕ РАЗЛОЖЕНИЕ КАК МОДЕЛЬ СМЕСИ РАСПРЕДЕЛЕНИЙ

Рассмотрим альтернативный подход к параметризации данных с помощью BTD. Первым отличием является вероятностный характер модели: векторизация наблюдаемых данных предполагается реализацией случайного вектора ${\mathbf{x}}$, зависящего от латентного (скрытого) представления; процесс моделирования состоит в детализации функции плотности распределения $p({\mathbf{x}})$, а результат обучения потенциально применим также и для генерации новых данных (генеративная модель). Другое важное отличие состоит в том, что индивидуальные параметры модели специфичны для каждой отдельной группы наблюдений (кластера) в предположении, что известно точное число кластеров $K$, а не для каждого отдельного представителя.

4.1. Параметризованные смеси распределений

В общих чертах, предложенная вероятностная модель на основе BTD возникает в результате параметризации смеси непрерывных латентных моделей в предположении о структурированности матриц смешивания. Частным случаем подобных смесей является смесь вероятностных главных компонент (mixture of PPCAs, см. [22]). Базовая модель вероятностных главных компонент (PPCA) (см. [23]) представляет собой параметризацию наблюдений ${\mathbf{x}}$ с помощью скрытого (латентного) представления, имеющего стандартное нормальное распределение, ${\mathbf{z}} \sim N({\mathbf{0}},I)$: $p({\mathbf{x}}) = \int {p({\mathbf{x}}\,{\text{|}}\,{\mathbf{z}})p({\mathbf{z}})d{\mathbf{z}}} $, $p({\mathbf{x}}\,{\text{|}}\,{\mathbf{z}}) = N(\mu + W{\mathbf{z}},{{L}^{{ - 1}}})$, где $\mu $ – постоянная составляющая среднего (сдвиг), $W$ – матрица смешивания, ${{L}^{{ - 1}}} = {{\sigma }^{2}}I$ – ковариационная матрица (в более общем случае факторного анализа ${{L}^{{ - 1}}}$ не обязательно является скалярной). В модели смеси вероятностных непрерывных скрытых представлений предполагается, что $p({\mathbf{x}})$ является мультимодальным распределением (см. [24]), порожденным вероятностной смесью нескольких распределений, каждая мода которого параметризована своим непрерывным скрытым представлением:

(10)
$p({\mathbf{x}}) = \sum\limits_{k = 1}^K {{{\pi }_{k}}{{p}_{k}}({\mathbf{x}})} ,\quad {{p}_{k}}({\mathbf{x}}) = \int {{{p}_{k}}({\mathbf{x}}\,{\text{|}}\,{{{\mathbf{z}}}_{k}})p({{{\mathbf{z}}}_{k}})d{{{\mathbf{z}}}_{k}}} ,\quad 0 \leqslant {{\pi }_{k}} \leqslant 1,\quad \sum\limits_{k = 1}^K {{{\pi }_{k}} = 1} ,$
где ${{p}_{k}}({\mathbf{x}}\,{\text{|}}\,{{{\mathbf{z}}}_{k}}) = N({{\mu }_{k}} + {{W}_{k}}{{{\mathbf{z}}}_{k}},L_{k}^{{ - 1}})$, $p({{{\mathbf{z}}}_{k}}) = N({\mathbf{0}},I)$. Если проводить аналогии с разд. 3, то в среднем наблюдения ${\mathbf{x}}$ предполагаются сгенерированными в результате взвешенной суммы нескольких процессов линейного смешивания. Коэффициенты $\{ {{\pi }_{k}}\} _{{k = 1}}^{K}$ выступают в роли априорных вероятностей принадлежности к тому или иному кластеру. Соответствующие апостериорные вероятности вычисляются по правилу Байеса и имеют вид
(11)
${{\gamma }_{k}}({\mathbf{x}}) = p(k\,{\text{|}}\,{\mathbf{x}}) = \frac{{{{\pi }_{k}}{{p}_{k}}({\mathbf{x}})}}{{\sum\nolimits_l \,{{\pi }_{l}}{{p}_{l}}({\mathbf{x}})}}.$
Обратим также внимание, что априорные вероятности принадлежности к тому или иному кластеру в построенной модели были параметризованы через исходные наблюдения и параметры кластера ${{\pi }_{k}} \equiv {{\pi }_{k}}({\mathbf{x}},{{W}_{k}})$ (подробнее см. следующий п. 4.2.).

Дальнейшая связь с BTD устанавливается следующим образом: для каждого $k = \overline {1,K} $ положим ${{\mu }_{k}} = 0$, обозначим реализацию вектора ${{{\mathbf{z}}}_{{\mathbf{k}}}}$ как ${{{\mathbf{Z}}}_{{\mathbf{k}}}}$ и рассмотрим ${{W}_{k}}{{{\mathbf{Z}}}_{k}}$ как векторизацию некоторого тензора, структура которого зависит от выбора того или иного тензорного формата, например:

– канонический (полилинейный, CP):

(12)
$\begin{array}{*{20}{c}} {{{{\mathbf{Z}}}_{k}} = {\text{vec}}\left( {\operatorname{diag} ({{{\widehat \lambda }}_{k}})} \right) \in {{\mathbb{R}}^{{{{R}_{k}} \times 1}}},\quad {{W}_{k}} = C_{k}^{{(d)}} \odot \ldots \odot C_{k}^{{(1)}}} \end{array};$

– (Lr, 1) с $1 < P < d$, ${{{\mathbf{L}}}_{k}} = ({{L}_{{k;1}}}, \ldots ,{{L}_{{k;{{R}_{k}}}}})$, ${{M}_{k}} = \sum\nolimits_i \,{{L}_{{k;i}}}$ (CP специального вида):

(13)
$\begin{gathered} \begin{array}{*{20}{c}} {{{{\mathbf{Z}}}_{{\mathbf{k}}}} \in {{\mathbb{R}}^{{{{M}_{k}} \times 1}}},\quad {{W}_{k}} = \widetilde C_{k}^{{(d)}} \odot \ldots \odot \widetilde C_{k}^{{(1)}},\quad \widetilde C_{k}^{{(p)}} = \left\{ \begin{gathered} C_{k}^{{(p)}} \in {{\mathbb{R}}^{{{{n}_{p}} \times {{M}_{k}}}}},\quad p \leqslant P, \hfill \\ \underbrace {C_{k}^{{(p)}}{{E}_{k}}}_{},\quad p > P; \hfill \\ \end{gathered} \right.} \end{array} \hfill \\ _{{{\text{повтор}}\;{\text{столбцов}}\;C_{k}^{{(p)}} \in {{\mathbb{R}}^{{{{n}_{p}} \times {{R}_{k}}}}}}} \hfill \\ \end{gathered} $

– Таккера с ядром в роли латентного представления (“Tucker-core”):

(14)
$\begin{array}{*{20}{c}} {{{{\mathbf{Z}}}_{k}} = \operatorname{vec} ({{{\widehat G}}_{k}}) \in {{\mathbb{R}}^{{{{r}_{1}} \cdot \ldots \cdot {{r}_{d}} \times 1}}},\quad {{W}_{k}} = A_{k}^{{(d)}} \otimes \ldots \otimes A_{k}^{{(1)}}} \end{array};$

– Таккера с фактор-матрицей в роли латентного представления (“Tucker-factor”), $G \in {{\mathbb{R}}^{{{{r}_{0}} \times {{r}_{1}} \times \ldots \times {{r}_{d}}}}}$:

(15)
$\begin{array}{*{20}{c}} {{{{\mathbf{Z}}}_{k}} = \operatorname{vec} ({\mathbf{\hat {a}}}_{k}^{{(0)}}) \in {{\mathbb{R}}^{{{{r}_{0}} \times 1}}},\quad {{W}_{k}} = \left( {(A_{k}^{{(d)}} \otimes \ldots \otimes A_{k}^{{(1)}})\mathop {({{{({{G}_{k}})}}_{{(0)}}})}\nolimits^{\text{T}} } \right) \otimes I;} \end{array}$

– тензорного поезда (TT) [25] с ядрами $G_{k}^{{(i)}} \in {{\mathbb{R}}^{{{{r}_{{i - 1}}} \cdot {{n}_{i}} \times {{r}_{i}}}}}$, ${{I}_{{{{N}_{{\{ \ldots i\} }}}}}} = {{I}_{{{{n}_{{i + 1}}} \ldots {{n}_{d}}}}}$:

(16)
$\begin{array}{*{20}{c}} {{{{\mathbf{Z}}}_{k}} \in {{\mathbb{R}}^{{{{r}_{d}} \times 1}}},\quad {{W}_{k}} = ({{I}_{{{{N}_{{\{ \ldots 1\} }}}}}} \otimes G_{k}^{{(1)}}) \ldots ({{I}_{{{{N}_{{\{ \ldots d - 1\} }}}}}} \otimes G_{k}^{{(d - 1)}})G_{k}^{{(d)}}.} \end{array}$

Таким образом, правдоподобие модели вероятностного BTD-смешивания принимает следующий вид:

(17)
$p({\mathbf{x}}\,{\text{|}}\,{{{\mathbf{z}}}_{1}}, \ldots ,{{{\mathbf{z}}}_{K}}) = \sum\limits_{k = 1}^K {{{\pi }_{k}}N({\mathbf{x}}\,{\text{|}}\,{{W}_{k}}{{{\mathbf{z}}}_{k}},L_{k}^{{ - 1}})} ,\quad {{W}_{k}}{{{\mathbf{z}}}_{k}} = \operatorname{vec} ({{T}_{k}}({{{\mathbf{z}}}_{k}})),$
ковариационные матрицы $L_{k}^{{ - 1}}$ полагаются либо имеющими скалярный вид (изотропный шум), либо диагональными. Обратная связь модели с классическим BTD разложением может быть прослежена через условное математическое ожидание ${\mathbf{x}}$ при известных реализациях скрытых представлений $\{ {{{\mathbf{z}}}_{k}}\} _{{k = 1}}^{K}$, результат соответствует развертке тензора в BTD формате с ${{\pi }_{k}}$ в роли весов для слагаемых:
(18)
$E[{\mathbf{x}}\,{\text{|}}\,{{{\mathbf{z}}}_{1}}, \ldots ,{{{\mathbf{z}}}_{K}}] = \sum\limits_{k = 1}^K {{{\pi }_{k}}{{W}_{k}}{{{\mathbf{z}}}_{k}}} = \operatorname{vec} (\widehat T).$
Помимо введенного правдоподобия, требуется также определить распределение для скрытых представлений $\{ {{{\mathbf{z}}}_{k}}\} _{{k = 1}}^{K}$. В силу того что правдоподобие моделируется с помощью нормального распределения, в качестве априорных распределений $p({{{\mathbf{z}}}_{k}})$ также выберем нормальное (как сопряженное к правдоподобию, см. [24]):
(19)
$p({{{\mathbf{z}}}_{k}}) = N({{{\mathbf{m}}}_{k}},\Lambda _{k}^{{ - 1}}) = N({{{\mathbf{m}}}_{k}}({\mathbf{x}},{{W}_{k}}),\Lambda _{k}^{{ - 1}}).$
Отличием от случая PPCA является использование диагональных матриц $\Lambda _{k}^{{ - 1}}$ и ненулевых средних ${{{\mathbf{m}}}_{k}}$. Как и в случае вариационных автоэнкодеров (см. [26]), для средних ${{{\mathbf{m}}}_{k}}$ была введена зависимость от входных данных и параметров разложения, в данном случае ${{{\mathbf{m}}}_{k}}({\mathbf{x}},{{W}_{k}}) = {{(W_{k}^{{\text{T}}}{{W}_{k}})}^{{ - 1}}}W_{k}^{{\text{T}}}{\mathbf{x}} = W_{k}^{ + }{\mathbf{x}}$. Такой выбор объясняется следующими соображениями: рассмотрим условное среднее $E{\kern 1pt} [{\mathbf{x}}\,{\text{|}}\,k] = {{W}_{k}}{{{\mathbf{m}}}_{k}}$, и для конкретной реализации ${\mathbf{x}}$ данное выражение есть ни что иное как проекция в пространство столбцов матрицы ${{W}_{k}}$, ${{W}_{k}}{{{\mathbf{m}}}_{k}}({\mathbf{x}},{{W}_{k}}) = {{P}_{{{{W}_{k}}}}}{\mathbf{x}}$, т.е. в итоге кластерное среднее определяется соответствующей проекцией наблюдения.

Параметры ковариационных матриц $\Lambda _{k}^{{ - 1}}$ и $L_{k}^{{ - 1}}$ полагались реализациями из лог-нормальных распределений с нулевым средним и единичной дисперсией. Наконец, укажем итоговое выражение для модели вариационного BTD разложения (VBTD) с нормальными априорным распределением и правдоподобием:

(20)
$\begin{gathered} p({\mathbf{x}},\{ {{{\mathbf{z}}}_{k}},{{\Lambda }_{k}},L_{k}^{{ - 1}}\} _{{k = 1}}^{K}) = \sum\limits_{k = 1}^K {{{p}_{k}}N(lnL_{k}^{{ - 1}}\,{\text{|}}\,0,I)N(ln\Lambda _{k}^{{ - 1}}\,{\text{|}}\,0,I)} , \\ {{p}_{k}} = {{\pi }_{k}}({\mathbf{x}},{{W}_{k}})N({\mathbf{x}}\,{\text{|}}\,{{W}_{k}}{{{\mathbf{z}}}_{k}},L_{k}^{{ - 1}})N({{{\mathbf{z}}}_{k}}\,{\text{|}}\,{{{\mathbf{m}}}_{k}}({\mathbf{x}},{{W}_{k}}),\Lambda _{k}^{{ - 1}}),\quad {{W}_{k}}{{{\mathbf{z}}}_{k}} = \operatorname{vec} ({{T}_{k}}({{{\mathbf{z}}}_{k}})). \\ \end{gathered} $
Текущая реализация предполагает, что параметры $\{ {{W}_{k}}\} _{{k = 1}}^{K}$ являются числовыми величинами, однако возможно построение более общей версии разложения, если моделировать их как случайные величины.

Аналогично детерминистической модели, введем групповую часть разложения, ${{{\mathbf{z}}}_{g}}$, а также введем обозначение ${{{\mathbf{\hat {z}}}}_{k}} = \mathop {[{{{\mathbf{z}}}_{k}}{{{\mathbf{z}}}_{g}}]}\nolimits^{\text{T}} $. Полагая независимость обеих частей скрытого представления ${{{\mathbf{z}}}_{g}} \bot {{{\mathbf{z}}}_{k}}$, получаем $p({{{\mathbf{\hat {z}}}}_{k}}) = p({{{\mathbf{z}}}_{k}},{{{\mathbf{z}}}_{g}}) = p({{{\mathbf{z}}}_{g}})p({{{\mathbf{z}}}_{k}})$, и используя свойства нормального распределения, приходим к следующему виду правдоподобия:

(21)
${{p}_{k}}({\mathbf{x}}\,{\text{|}}\,{{{\mathbf{\hat {z}}}}_{k}}) = N({\mathbf{x}}\,{\text{|}}\,{{W}_{k}}{{{\mathbf{z}}}_{k}} + {{W}_{g}}{{{\mathbf{z}}}_{g}},L_{k}^{{ - 1}} + L_{g}^{{ - 1}}),\quad k = \overline {1,K} ,$
где ${{W}_{g}}$ – параметры разложения, отвечающие общему скрытому представлению ${{{\mathbf{z}}}_{g}}$, $L_{g}^{{ - 1}}$ – ковариационная матрица для группового слагаемого.

4.2. Априорное решающее правило

Рассмотрим подробнее априорные вероятности принадлежности к кластерам, $\{ {{\pi }_{k}}\} _{{k = 1}}^{K}$. В общем виде они могут иметь сложную структуру, учитывающую дополнительные зависимости, и одной из наиболее важных является зависимость от входных данных ${\mathbf{x}}$, ${{\pi }_{k}} \equiv {{\pi }_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}})$. От подобной функциональной зависимости требуется выполнение условия нормировки, т.е. ${{\pi }_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}}) \geqslant 0$, $\sum\nolimits_{k = 1}^K \,{{\pi }_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}}) = 1$. Распространенный способ учета такого условия состоит в использовании многопеременной логистической функции (softmax) на результате отображения некоторой (дифференцируемой) функции ${{a}_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}})$ (см. [24]):

(22)
${{\pi }_{k}} \equiv {{\pi }_{k}}({\mathbf{x}},\theta ) = {\text{softma}}{{{\text{x}}}_{k}}\left( {a({\mathbf{x}},{{\theta }_{\pi }})} \right) = \frac{{exp\left( {{{a}_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}})} \right)}}{{\sum\limits_i {exp\left( {{{a}_{l}}({\mathbf{x}},{{\theta }_{{{{\pi }_{l}}}}})} \right)} }}.$

Смеси распределений с параметризованными априорными вероятностями известны как модели коллектива экспертов (см. [24]). В случае, если для нормализации используется функция ${\text{softmax}}(a({\mathbf{x}},{{\theta }_{\pi }}))$, отображение $a({\mathbf{x}},{{\theta }_{\pi }})$ должно быть задано таким образом, чтобы значения элементов вектора-результата были прямо пропорциональны близости наблюдения ${\mathbf{x}}$ к соответствующим кластерам. Выбор параметризации для предложенной модели на основе BTD был мотивирован результатами работы [19] и нашим опытом из [27], в которой аналогичное отображение использовалось для построения классификаторов на основе разложений Таккера, которые, в частности, позволили получить более устойчивые результаты при изменении оборудования сбора данных и метода химической экстракции в сравнении с рядом иных рассмотренных классификаторов.

Параметризация основана на главном угле (см. [28]) между пространствами источников (столбцы фактор-матрицы для соответствующей оси) и соответствующей матризацией наблюдения. Пусть имеется $d$-мерное тензорное наблюдение ${{X}_{i}}$, матрица развертки параметров $A_{k}^{{(m)}}$ при выбранной моде источников $m \in \{ 1,2, \ldots ,d\} $, $k = \overline {1,K} $, и их сингулярные разложения, ${{({{X}_{i}})}_{{(m)}}} = {{U}_{i}}{{\Sigma }_{i}}V_{i}^{{\text{T}}}$ и $A_{k}^{{(m)}} = {{U}_{k}}{{\Sigma }_{k}}V_{k}^{{\text{T}}}$, тогда отображение ${{a}_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}})$ может быть вычислено следующим образом:

(23)
${{a}_{k}}({{X}_{i}},A_{k}^{{(m)}}) = {{{\mathbf{\tilde {u}}}}^{T}}\Sigma _{k}^{{ - 1}}V_{k}^{{\text{T}}}{{(A_{k}^{{(m)}})}^{{\text{T}}}}{{({{X}_{i}})}_{{(m)}}}{{V}_{i}}\Sigma _{i}^{{ - 1}}{\mathbf{\tilde {v}}} = {{{\mathbf{\tilde {u}}}}^{T}}U_{k}^{{\text{T}}}{{U}_{i}}{\mathbf{\tilde {v}}} = {{\sigma }_{{max}}}(U_{k}^{{\text{T}}}{{U}_{i}}),$
где ${\mathbf{\tilde {u}}}$, ${\mathbf{\tilde {v}}}$ – левый и правый сингулярные векторы для матрицы $U_{k}^{{\text{T}}}{{U}_{i}}$, отвечающие ее максимальному сингулярному значению ${{\sigma }_{{max}}}({\kern 1pt} \cdot {\kern 1pt} )$.

5. ОБУЧЕНИЕ МОДЕЛЕЙ

5.1. Нелинейная задача наименьших квадратов

Для предложенных детерминистических моделей GLRO и GTLD вычисление параметров производится через условную оптимизацию нелинейных наименьших квадратов. Для заданного набора $d$-мерных тензоров $\{ {{X}_{i}}\} _{{i = 1}}^{N}$ одинаковых размеров, объединенных вдоль групповой оси в $(d + 1)$-мерный тензор $X$, задача поиска параметров $\theta $ принимает следующий вид:

(24)
$\begin{gathered} \mathop {min}\limits_\theta \left\| {X - T(\theta )} \right\|_{F}^{2} = \mathop {min}\limits_{{{\theta }_{g}},\{ {{\theta }_{i}}\} _{{i = 1}}^{N}} \left\| {X - \left( {T({{\theta }_{g}}) + \sum\limits_{i = 1}^N \,T({{\theta }_{i}})} \right)} \right\|_{F}^{2} \\ {\text{s}}.{\text{t}}.\;\;h({{\theta }_{g}},\{ {{\theta }_{i}}\} _{{i = 1}}^{N}) = 0, \\ \end{gathered} $
где через ${{\theta }_{g}}$, $\{ {{\theta }_{i}}\} _{{i = 1}}^{N}$ обозначены параметры общих и индивидуальных частей соответственно (согласно (6) или (7)), и в выражении $h({{\theta }_{g}},\{ {{\theta }_{i}}\} _{{i = 1}}^{N}) = 0$ собраны условия на фактор-матрицы групповой моды и условия отделимости (9).

Среди существующих методов оптимизации параметров в задаче нелинейных наименьших квадратов был использован ряд методов первого и второго порядка (выражения для градиентов и матрицы Гессе приведены в п. 10.2., 10.3.), а также метод попеременных наименьших квадратов (ALS). Для вычисления параметров различных тензорных разложений существуют готовые Матлаб пакеты (например, Tensorlab, см. [29]), однако для проведения исследований из области машинного обучения более распространено использование языка Python, на котором и был реализован вычислительный код для настоящей работы.

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

5.2. Стохастический вариационный вывод (SVI)

Вероятностные модели способны ассимилировать новые наблюдения и делать предсказания с помощью статистического вывода. Для моделей со скрытым представлением ${\mathbf{z}}$ в общем случае для этого требуется вычисление апостериорного распределения $p({\mathbf{z}}\,{\text{|}}\,{\mathbf{x}}) = p({\mathbf{x}},{\mathbf{z}}){\text{/}}p({\mathbf{x}})$. Часто проблему представляет вычисление $p({\mathbf{x}})$, и для статистического вывода используются различные приближения. Вариационный вывод представляет собой метод преобразования задачи вывода к задаче оптимизации. Основная идея заключается в замене точного апостериорного распределения на некоторого представителя из определенного семейства распределений, который как можно точнее приближает истинное апостериорное (подробнее о методе стохастического вариационного вывода см., например, [30]).

Непосредственная оптимизация параметров распределения производилась через максимизацию нижней грани правдоподобия (ELBO) (см. [31]). Для заданных моделью правдоподобия $p({\mathbf{x}}\,{\text{|}}\,{\mathbf{z}},{{\theta }_{{{\text{lh}}}}})$ и априорного распределения $p({\mathbf{z}}\,|\,{{\theta }_{z}})$, ${{p}_{\theta }}({\mathbf{x}},{\mathbf{z}}) = p({\mathbf{x}}\,|\,{\mathbf{z}},{{\theta }_{{{\text{lh}}}}})p({\mathbf{z}}\,|\,{{\theta }_{z}})$, а также для вариационного распределения $q({\mathbf{z}}\,|\,\phi ) = {{q}_{\phi }}({\mathbf{z}})$ вариационная нижняя грань возникает в уравнении $ln{{p}_{\theta }}({\mathbf{x}}) = {\text{ELBO}} + {\text{KL}}\left( {{{q}_{\phi }}({\mathbf{z}})\left. {\parallel {{p}_{\theta }}({\mathbf{z}}\,|\,{\mathbf{x}})} \right)} \right.$ и определяется как

(25)
${\text{ELBO}} \equiv {{E}_{{{{q}_{\phi }}({\mathbf{z}})}}}[ln{{p}_{\theta }}({\mathbf{x}},{\mathbf{z}}) - ln{{q}_{\phi }}({\mathbf{z}})],$
${\text{KL}}({{q}_{\phi }}({\mathbf{z}})\,||\,{{p}_{\theta }}({\mathbf{z}}\,|\,{\mathbf{x}})) = \int {{{q}_{\phi }}({\mathbf{z}})} ln\tfrac{{{{q}_{\phi }}({\mathbf{z}})}}{{{{p}_{\theta }}({\mathbf{z}}\,|\,{\mathbf{x}})}}d{\mathbf{z}}$ – дивергенция Кульбака–Лейблера, статистическая мера различий между двумя распределениями. Важно, что метод ELBO поддерживает обучение по подвыборкам и оценки Монте-Карло для стохастического градиента. Для реализации предложенной модели и оптимизации ее параметров методом стохастического вариационного вывода был использован пакет Pyro (см. [32]), в котором вычисления базируются на построении динамических вычислительных графов с помощью pytorch (см. [33]).

Для модели VBTD (20) было использовано следующее вариационное распределение:

(26)
$\begin{array}{*{20}{c}} {p(\{ {{{\mathbf{z}}}_{k}},{{\Lambda }_{k}},\sigma _{k}^{2}\} _{{k = 1}}^{K}) = \prod\limits_{k = 1}^K {\gamma _{k}^{{{{z}_{{{{\pi }_{k}}}}}}}} \cdot \prod\limits_{k = 1}^K {N({{{\mathbf{z}}}_{k}}\,|\,{{{{\mathbf{\hat {m}}}}}_{k}},\hat {\Lambda }_{k}^{{ - 1}})} ,} \\ {{{\gamma }_{k}} = \frac{{{{p}_{k}}}}{{\sum\limits_l {{{p}_{l}}} }},\quad {{z}_{{{{\pi }_{k}}}}} \in \{ 0,1\} ,\quad \sum\limits_{k = 1}^K {{{z}_{{{{\pi }_{k}}}}}} = 1,\quad p({{z}_{{{{\pi }_{k}}}}} = 1) = {{\pi }_{k}},\quad p({{{\mathbf{z}}}_{\pi }}) = \prod\limits_{k = 1}^K {\pi _{k}^{{{{z}_{{{{\pi }_{k}}}}}}}} ,} \end{array}$
где $N({{{\mathbf{z}}}_{k}}\,|\,{{\widehat {\mathbf{m}}}_{k}},\widehat \Lambda _{k}^{{ - 1}})$ – истинные апостериорные распределения для каждой отдельной модели непрерывного смешивания на основе нормального распределения.

6. КРАТКИЙ ОБЗОР ИНЫХ МОДЕЛЕЙ

В [34] представлена практическая методология использования специальных матричных разложений для группового анализа данных. Calhoun и др. строят групповой анализ данных функциональной магнитно-резонансной томографии (фМРТ) на основе метода независимых компонент (ICA). В деталях обсуждаются проблемы выбора числа компонент, моделирования зависимости данных от общих компонент, а также статистической валидации результатов. Метод ориентирован на выделение только общей активности, включает в себя двухшаговое выделение главных компонент (на уровне отдельных субъектов и на групповом уровне, для объединения найденных индивидуальных главных компонент) с последующим выделением независимых компонент. Вкратце модель выглядит следующим образом:

(27)
$\mathop {{{X}_{i}} \approx {{Y}_{i}}V_{i}^{{\text{T}}}}\limits_{({\text{PCA}},\;{\text{инд}}.)} ,\quad \mathop {Z \approx [{{Y}_{1}}, \ldots ,{{Y}_{N}}]{\kern 1pt} W}\limits_{({\text{PCA}},\;{\text{общ}}.)} ,\quad \mathop {Z \approx S{{A}^{{\text{T}}}}}\limits_{({\text{ICA}})} ,\quad {{X}_{i}} \approx S{{A}^{{\text{T}}}}W_{i}^{{\text{T}}}V_{i}^{{\text{T}}} = S{{({{V}_{i}}{{W}_{i}}A)}^{{\text{T}}}} = SM_{i}^{{\text{T}}},$
где столбцы матрицы $S$ – источники, ${{M}_{i}}$ – матрица смешивания для наблюдения ${{X}_{i}}$, $i = \overline {1,N} $.

Несколькими годами позже появился обзор [18] прикладных исследований, использующих групповой метод независимых компонент (GICA), в частности, для анализа фМРТ данных, вызванных потенциалов, данных точечного нуклеотидного полиморфизма. В частности, в обзоре перечислены исследования, в которых метод независимых компонент адаптирован к задаче группового анализа данных; подходы сгруппированы в пять категорий: (1) построение отдельных факторизаций для каждого отдельного субъекта с последующим поиском взаимосвязей между ними; (2)–(3) объединение данных по одной из осей (пространственный/временной ICA); (4) усреднение данных по групповой оси; (5) тензорные разложения. Кроме того, авторы демонстрируют возможность использования ICA для мультимодального (т.е. включающего разнородные данные) анализа через модели общего (joint) и параллельного (parallel) ICA.

Другая известная модель, тензорный вероятностный метод независимых компонент (TPICA), основана на совмещении канонического разложения и вероятностного ICA (см. [10]). В данном подходе нет строгого разделения параметров на общие и индивидуальные, а матрица смешивания имеет структуру произведения Хатри–Рао между фактор-матрицами групповой оси и оси времени.

Позднее данные методы были рассмотрены с более общей точки зрения Guo, Pagnoni (см. [35]). В работе предложен модифицированный EM-алгоритм для максимизации правдоподобия с учетом трех вариантов структуры матрицы смешивания:

(28)
${{M}_{{{\text{TPICA}}}}} = {{M}^{{{{{\text{m}}}_{{{\text{subj}}}}}}}} \otimes {{M}^{{{{{\text{m}}}_{{{\text{time}}}}}}}},\quad {{M}_{{{\text{GICA}}}}} = [\begin{array}{*{20}{c}} {{{M}_{1}}}& \ldots &{{{M}_{N}}} \end{array}],\quad {{M}_{{{\text{gTPICA}}}}} = \left[ {\begin{array}{*{20}{c}} {M_{1}^{{{{{\text{m}}}_{{{\text{subj}}}}}}} \otimes M_{1}^{{{{{\text{m}}}_{{{\text{time}}}}}}}} \\ \vdots \\ {M_{N}^{{{{{\text{m}}}_{{{\text{subj}}}}}}} \otimes M_{N}^{{{{{\text{m}}}_{{{\text{time}}}}}}}} \end{array}} \right].$
Дополнительно авторами рассмотрена методология выбора оптимальной структуры матрицы смешивания на основе критерия отношения правдоподобия, валидация которой была произведена на фМРТ данных.

Работа [36] посвящена построению группового анализа данных электрической активности мозга (ЭЭГ) на основе неотрицательных матричных разложений с учетом как групповой, так и индивидуальных составляющих: разделение индивидуальных фактор-матриц $\{ S_{i}^{{(I)}}\} _{{i = 1}}^{N}$, как и увеличение степени схожести между групповыми фактор-матрицами $\{ S_{i}^{{(C)}}\} _{{i = 1}}^{N}$, достигалось с помощью введения соответствующих слагаемых в оптимизируемый функционал:

(29)
$\begin{gathered} \mathop {min}\limits_{S_{i}^{{({\text{C}})}},S_{i}^{{({\text{I}})}},\left. {{{M}_{i}}} \right|_{{i = 1}}^{N}} \sum\limits_{i = 1}^N {\left[ {\left\| {{{X}_{i}} - [S_{i}^{{({\text{C}})}},S_{i}^{{({\text{I}})}}]M_{i}^{{\text{T}}}} \right\|_{F}^{2} + \mathop {\gamma \left\| {[S_{i}^{{({\text{C}})}},S_{i}^{{({\text{I}})}}]} \right\|_{F}^{2}}\limits_{} } \right.} + \\ \, + \left. {\sum\limits_{j = 1}^{i - 1} {\left( {\alpha \left\| {S_{i}^{{({\text{C}})}} - S_{j}^{{({\text{C}})}}} \right\|_{F}^{2} - \beta \left\| {S_{i}^{{({\text{I}})}} - S_{j}^{{({\text{I}})}}} \right\|_{F}^{2}} \right)} } \right],\quad {\text{s}}.{\text{t}}.\quad S_{i}^{{({\text{C}})}},S_{i}^{{({\text{I}})}},\left. {{{M}_{i}}} \right|_{{i = 1}}^{N} \geqslant 0. \\ \end{gathered} $

В [37] представлен вариант модели группового анализа данных на основе неотрицательного разреженного канонического разложения и концепта “связанных тензорных разложений”, в которой с помощью иерархического ALS алгоритма (вариант блочно-координатного спуска, каждая итерация которого состоит в последовательном обновлении строк/столбцов фактор-матриц, см. [38], [39]) оцениваются статистически независимые индивидуальные компоненты. Оптимизационная задача для вычисления разложения имеет вид (условия неотрицательности и разреженности опущены):

(30)
$\begin{gathered} min\sum\limits_{i = 1}^N {\left\| {{{X}_{i}} - \left| {\left[ {\Lambda ;U_{1}^{{(i)}}, \ldots ,U_{d}^{{(i)}}} \right]} \right|{\kern 1pt} } \right\|_{F}^{2}} \\ {\text{s}}.{\text{t}}.\quad \left. {U_{k}^{{(i)}}} \right|_{{i = 1}}^{N} = {{U}_{k}},\quad {{\left\| {{{U}_{k}}[:,r]} \right\|}_{F}} = 1,\quad r = \overline {1,R} ,\quad k = \overline {1,K} ,\quad K < d,\quad \Lambda = \operatorname{Diag} ({{\lambda }_{1}}, \ldots ,{{\lambda }_{R}}). \\ \end{gathered} $

В [40] была предложена модель на основе разложения Таккера с двумя фактор-матрицами для случая 3D данных. Общие компоненты моделируются как первые $C$ столбцов левой сингулярной матрицы соответствующей фактор-матрицы, а индивидуальные компоненты вычисляются как первые ${{I}_{k}}$ столбцов левой сингулярной матрицы для проекции той же фактор-матрицы на пространство, ортогональное пространству общих компонент.

В [19] был предложен метод выделения общей и индивидуальной активностей на основе матричной модели (COBE). Предполагаемая структура данных схожа с рассматриваемой в [36], но модель выделяет общие источники иным образом:

(31)
$\begin{gathered} \mathop {min}\limits_{\hat {A},{{A}_{i}}} \sum\limits_{i = 1}^N {\left\| {{{X}_{i}} - \hat {A}\hat {B}_{i}^{T} - {{A}_{i}}B_{i}^{{\text{T}}}} \right\|_{F}^{2}} \\ {\text{s}}{\text{.t}}{\text{.}}\quad {{{\hat {A}}}^{{\text{T}}}}\hat {A} = {{I}_{C}},\quad A_{i}^{{\text{T}}}{{A}_{i}} = {{I}_{{{{n}_{2}} - C}}},\quad {{{\hat {A}}}^{{\text{T}}}}{{A}_{i}} = 0,\quad i = \overline {1,N} , \\ \end{gathered} $
где $\hat {A}$ – общие источники ($C$ столбцов), ${{A}_{i}}$ – индивидуальные источники.

Из существующих расширений BTD выделим вариант разложения со слагаемыми в формате Таккера, предложенный в [41] для случая совмещенных фактор-матриц. Предложенные вычислительные алгоритмы опираются на расширенный QZ алгоритм и метод Якоби. Кроме того, важно отметить, что идея использования BTD разложения с разным типом слагаемых не является новой, например, в [42] были использованы одновременно слагаемые в формате Таккера и каноническом формате для построения приближений классических потенциалов.

7. ВЫЧИСЛИТЕЛЬНЫЕ ЭКСПЕРИМЕНТЫ

7.1. Программная реализация

Все модели и вычислительные эксперименты были реализованы на языке программирования Python с использованием дистрибутива Anaconda (см. [43]), который включает в себя различные пакеты для научных вычислений. В данном исследовании были использованы следующие пакеты: numpy [44], scipy [45], pandas [46], scikit-learn [47], matplotlib [48], seaborn [49], pyro [32]. Весь код опубликован в следующих репозиториях: https://github.com/kharyuk/gbtd; https://github.com/kharyuk/vbtd. Эксперименты были систематизированы с помощью Jupyter Notebooks (см. [50]).

Фиг. 1.

ETH80, изображения разных классов (a) и выборочные ракурсы 10 объектов класса “чашка” (б). Набор данных доступен по адресу: http://datasets.d2.mpi-inf.mpg.de/eth80/eth80-contours.tgz

7.2. Наборы данных

ETH-80. ETH80 (см. [51]) состоит из изображений различных объектов, снятых с разных ракурсов. Всего доступно восемь классов: “яблоко”, “машина”, “корова”, “чашка”, “собака”, “лошадь”, “груша”, “томат”; каждый класс включает 10 объектов (наблюдений) (фиг. 1). Полученный в результате набор данных имеет размеры $({{N}_{{{\text{samples}}}}},{{N}_{{{\text{pixels}}}}},{{N}_{{{\text{angles}}}}},{{N}_{{{\text{colors}}}}})$, где ${{N}_{{{\text{angles}}}}}{\text{ = 41}}$. Разрешение изображений было понижено до $32 \times 32$, откуда ${{N}_{{{\text{pixels}}}}} = 1024$.

SMNI EEG (ERP). Набор данных SMNI EEG состоит из измерений электрической активности мозга по ${{N}_{{{\text{electrode}}}}}{\text{ = 64}}$ электродам для двух групп субъектов: страдающих алкогольной зависимостью и контрольной группы. Согласно описанию, данные были записаны с частотой 256 Гц и нарезаны на фрагменты длиной в 1 с, содержащие момент предъявления стимула. Данные были собраны для ${{N}_{{{\text{condition}}}}} = 3$ условий: предъявление одиночного стимула (изображения), двух совпадающих и двух различных стимулов. Данные были дополнительно усреднены по повторностям, помеченные как испорченные образцы были исключены, и итоговый набор данных был представлен как тензор размерами $({{N}_{{{\text{subject}}}}},{{N}_{{{\text{time}}}}},{{N}_{{{\text{electrode}}}}},{{N}_{{{\text{condition}}}}})$, где ${{N}_{{{\text{subject}}}}} = 119$, из них 76 субъектов из первой группы, 43 из контрольной (фиг. 2).

Фиг. 2.

Усредненные данные ЭЭГ для первого (a) и третьего (б) экспериментов; на изображениях сверху сигналы изображены в 2D виде,  где более светлый цвет соответствует большим значениям магнитуды; разные цвета на изображениях снизу отвечают разным каналам. Набор данных доступен по ссылке: https://archive.ics.uci.edu/ml/datasets/eeg+database

7.3. Оценка качества методов кластеризации

Кластеризация представляет собой одну из стандартных задач машинного обучения, в которой требуется сгруппировать различные образцы по заданному критерию близости. Более формально, указывается некоторый функционал близости $\rho ({{X}_{i}},{{X}_{j}})$, и для входных данных $\{ {{X}_{i}}\} _{{i = 1}}^{N}$ необходимо ввести разбиение на непересекающиеся подмножества (кластеры) таким образом, что объекты внутри одного кластера были более близки друг к другу, чем к образцам из прочих кластеров, в смысле $\rho ({{X}_{i}},{{X}_{j}})$.

Если для данных известны истинные метки, качество кластеризации можно измерить следующим образом. Рассмотрим два разбиения, $U = \{ {{U}_{i}}\} $, $V = \{ {{V}_{j}}\} $, где $U$ отвечает истинным меткам, и введем на них четыре стандартные величины: число истинно-положительных (TP), истинно-отрицательных (TN), ложноотрицательных (FN), ложноположительных (FP) пар:

(32)
$\begin{array}{*{20}{c}} {{\text{TP}} = \left| {\{ ({{a}_{k}},{{a}_{l}})\,|\,{{a}_{k}},{{a}_{l}} \in {{U}_{i}},{{a}_{k}},{{a}_{l}} \in {{V}_{j}}\} } \right|,\quad {\text{TN}} = \left| {\{ ({{a}_{k}},{{a}_{l}})\,|\,{{a}_{k}} \in {{U}_{{{{i}_{1}}}}},{{a}_{l}} \in {{U}_{{{{i}_{2}}}}},{{a}_{k}} \in {{V}_{{{{j}_{1}}}}},{{a}_{l}} \in {{V}_{{{{j}_{2}}}}}\} } \right|,} \\ {{\text{FP}} = \left| {\{ ({{a}_{k}},{{a}_{l}})\,|\,{{a}_{k}} \in {{U}_{{{{i}_{1}}}}},{{a}_{l}} \in {{U}_{{{{i}_{2}}}}},{{a}_{k}},{{a}_{l}} \in {{V}_{j}}\} } \right|,\quad {\text{FN}} = \left| {\{ ({{a}_{k}},{{a}_{l}})\,|\,{{a}_{k}},{{a}_{l}} \in {{U}_{i}},{{a}_{k}} \in {{V}_{{{{j}_{1}}}}},{{a}_{l}} \in {{V}_{{{{j}_{2}}}}}\} } \right|.} \end{array}$

Также определим следующие вспомогательные величины:

(33)
${{P}_{{ij}}} = \left| {{{U}_{i}} \cap {{V}_{j}}} \right|{\text{/}}N,\quad {{P}_{i}} = \left| {{{U}_{i}}} \right|{\text{/}}N,\quad {{P}_{j}} = \left| {{{V}_{j}}} \right|{\text{/}}N,\quad H(W) = E[ - logP(W)],$
где $W$ – некоторое разбиение, $H({\kern 1pt} \cdot {\kern 1pt} )$ – энтропия.

Используя величины (32), (33), определим три стандартных показателя качества, индекс Фолькс–Мэллоуз (FMI), как геометрическое среднее между точностью и полнотой, скорректированный индекс Рэнда (ARI), который схож с долей верных ответов, и скорректированная взаимная информация (AMI), $C_{N}^{2} = \left( {\begin{array}{*{20}{c}} N \\ 2 \end{array}} \right)$:

(34)
$\begin{gathered} {\text{RI}} = (TP + TN){\text{/}}C_{N}^{2},\quad {\text{MI}} = \sum\limits_{i = 1}^I {\sum\limits_{j = 1}^J {{{P}_{{ij}}}log\frac{{{{P}_{{ij}}}}}{{{{P}_{i}}{{P}_{j}}}}} } , \\ {\text{ARI }} = \frac{{{\text{RI}} - E\left[ {{\text{RI}}} \right]}}{{\max ({\text{RI}}) - E[{\text{RI}}]}}{\text{,}}\quad {\text{AMI}} = \tfrac{{{\text{MI}} - E[{\text{MI}}]}}{{max(H(U),H(V)) - E[{\text{MI}}]}}, \\ {\text{FMI}} = \sqrt {\tfrac{{TP}}{{TP + FP}}\tfrac{{TP}}{{TP + FN}}} . \\ \end{gathered} $

В силу того что модели GLRO и GTLD были построены как расширение модели, предложенной в [19] (будем ссылаться на нее по названию алгоритма COBE), модель COBE была также использована для сравнения. Другой альтернативой для сравнения качества была выбрана модель группового метода независимых компонент GICA (см. [18]). Оригинальная версия модели была разработана для выделения только общей активности; в наших экспериментах была воспроизведена процедура выделения общих независимых компонент, дополненная вычитанием восстановленных общих частей из исходных данных (т.е. контрастированием). Выделение индивидуальных частей в случае COBE было проведено способом, описанным в оригинальной статье, для GLRO и GTLD были использованы соответствующие слагаемые полного разложения. Для моделей были использованы разные сочетания общих рангов, ${{r}_{C}} = \overline {1,5} $, и рангов индивидуальных, ${{r}_{I}} = \overline {1,5} $ (кроме COBE, где используются только ${{r}_{C}}$). Для GICA ${{r}_{C}}$ соответствовал одновременно числу вторичных главных компонент и числу независимых компонент. В случае GTLD все ранги Таккера полагались равными ${{r}_{C}}$.

Контрастированные и необработанные данные были кластеризованы с помощью следующих алгоритмов: иерархическая агломеративная кластеризация HAC (см. [52]), K-средних (см. [52]) и кластеризации на основе модели смеси гауссовских распределений GMM (см. [52]) (с диагональными ковариационными матрицами). В табл. 1, 2 представлены результаты, полученные для ETH80 и SMNI EEG соответственно. Результаты для всех методов контрастирования были усреднены по 10 запускам; для алгоритмов K-средних/GMM результаты дополнительно усреднены по 20 разным инициализациям для каждого запуска.

Таблица 1.  

Сравнительные результаты моделей с подобранными гиперпараметрами для набора данных ETH80

    ARI AMI FMI
HAC ${\text{Ra}}{{{\text{w}}}^{1}}$ .568 .749 .638
  ${\text{COB}}{{{\text{E}}}^{1}}$ .643 ± .023 .774 ± .016 .691 ± .019
  ${\text{COB}}{{{\text{E}}}^{2}}$ .573 ± .013 .780 ± .005 .663 ± .009
  ${\text{GIC}}{{{\text{A}}}^{1}}$ .722 ± .000 .794 ± .000 .753 ± .000
  ${\text{GIC}}{{{\text{A}}}^{2}}$ .577 ± .000 .802 ± .000 .661 ± .000
  ${\text{GLR}}{{{\text{O}}}^{1}}$ .472 ± .093 .665 ± .060 .559 ± .074
  ${\text{GLR}}{{{\text{O}}}^{2}}$ .471 ± .098 .693 ± .076 .580 ± .067
  ${\text{GTL}}{{{\text{D}}}^{1}}$ .497 ± .071 .668 ± .047 .573 ± .054
  ${\text{GTL}}{{{\text{D}}}^{2}}$ .493 ± .072 .699 ± .043 .585 ± .047
Kmeans Raw .533 ± .068 .690 ± .047 .597 ± .055
  ${\text{COB}}{{{\text{E}}}^{3}}$ .571 ± .049 .729 ± .030 .635 ± .037
  ${\text{GIC}}{{{\text{A}}}^{3}}$ .559 ± .057 .712 ± .042 .624 ± .046
  ${\text{GLR}}{{{\text{O}}}^{2}}$ .461 ± .077 .645 ± .062 .543 ± .062
  ${\text{GLR}}{{{\text{O}}}^{3}}$ .445 ± .065 .645 ± .047 .535 ± .050
  ${\text{GTL}}{{{\text{D}}}^{3}}$ .473 ± .065 .652 ± .050 .550 ± .052
GMM Raw .475 ± .072 .639 ± .064 .549 ± .060
  ${\text{COB}}{{{\text{E}}}^{3}}$ .333 ± .262 .438 ± .341 .376 ± .293
  ${\text{GIC}}{{{\text{A}}}^{4}}$ .384 ± .229 .513 ± .301 .439 ± .258
  ${\text{GLR}}{{{\text{O}}}^{2}}$ .473 ± .079 .651 ± .065 .553 ± .063
  ${\text{GLR}}{{{\text{O}}}^{3}}$ .466 ± .072 .662 ± .054 .554 ± .054
  ${\text{GTL}}{{{\text{D}}}^{3}}$ .476 ± .073 .656 ± .056 .555 ± .058

Примечание. Формат ячеек таблицы: “$x \pm y$”, где $x$ – среднее значение, $y$ – стандартное отклонение. Для HAC за ${{R}_{{{\text{average}}}}}$ и ${{R}_{{{\text{complete}}}}}$ обозначены пересчет расстояний через усреднение и взятие максимума, за ${{\rho }_{{{\text{corr}}}}}$, ${{\rho }_{{cos}}}$, ${{\rho }_{{{\text{Canb}}.}}}$ – корреляция, косинусное расстояние и расстояние Канберра; ${\text{Ra}}{{{\text{w}}}^{1}}$: ${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{average}}}}}$; ${\text{COB}}{{{\text{E}}}^{1}}$: ${{r}_{c}} = 2$, ${{\rho }_{{cos}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{COB}}{{{\text{E}}}^{2}}$: ${{r}_{c}} = 2$, ${{\rho }_{{cos}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{average}}}}}$; ${\text{COB}}{{{\text{E}}}^{3}}$: ${{r}_{c}} = 4$; ${\text{GIC}}{{{\text{A}}}^{1}}$: ${{r}_{c}} = 2$, ${{r}_{i}} = \overline {3,5} $, ${{\rho }_{{cos}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GIC}}{{{\text{A}}}^{2}}$: ${{r}_{c}} = 2$, ${{r}_{i}} = \overline {4,5} $, ${{\rho }_{{{\text{Canb}}.}}}$, ${{R}_{{{\text{average}}}}}$; ${\text{GIC}}{{{\text{A}}}^{3}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 3$; ${\text{GIC}}{{{\text{A}}}^{4}}$: ${{r}_{c}} = 2$, ${{r}_{i}} = 2$; ${\text{GLR}}{{{\text{O}}}^{1}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 4$, ${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GLR}}{{{\text{O}}}^{2}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 3$, для HAC – ${{\rho }_{{{\text{Canb}}.}}}$, ${{R}_{{{\text{average}}}}}$; ${\text{GLR}}{{{\text{O}}}^{3}}$: ${{r}_{c}} = 5$, ${{r}_{i}} = 2$; ${\text{GTL}}{{{\text{D}}}^{1}}$: ${{r}_{c}} = 5$, ${{r}_{i}} = 3$, ${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GTL}}{{{\text{D}}}^{2}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 4$, ${{\rho }_{{{\text{Canb}}.}}}$, ${{R}_{{{\text{average}}}}}$; ${\text{GTL}}{{{\text{D}}}^{3}}$: ${{r}_{c}} = 5$, ${{r}_{i}} = 4$.

Таблица 2.  

Сравнительные результаты моделей с подобранными гиперпараметрами для набора данных SMNI EEG

    ARI AMI FMI
HAC ${\text{Ra}}{{{\text{w}}}^{1}}$ .161 .104 .601
  ${\text{COB}}{{{\text{E}}}^{1}}$ .092 ± .060 .057 ± .044 .577 ± .029
  ${\text{COB}}{{{\text{E}}}^{2}}$ .077 ± .053 .057 ± .032 .591 ± .047
  ${\text{GIC}}{{{\text{A}}}^{1}}$ .253 ± .000 .198 ± .000 .715 ± .000
  ${\text{GLR}}{{{\text{O}}}^{1}}$ .048 ± .075 .036 ± .044 .621 ± .059
  ${\text{GTL}}{{{\text{D}}}^{1}}$ .060 ± .078 .046 ± .037 .597 ± .052
  ${\text{GTL}}{{{\text{D}}}^{2}}$ .009 ± .047 .058 ± .060 .568 ± .038
Kmeans Raw .196 ± .015 .123 ± .009 .624 ± .009
  ${\text{COB}}{{{\text{E}}}^{2}}$ .168 ± .059 .106 ± .037 .608 ± .033
  ${\text{GIC}}{{{\text{A}}}^{2}}$ .236 ± .004 .153 ± .005 .646 ± .006
  ${\text{GLR}}{{{\text{O}}}^{2}}$ .039 ± .054 .027 ± .037 .676 ± .080
  ${\text{GTL}}{{{\text{D}}}^{2}}$ .126 ± .094 .083 ± .063 .674 ± .051
GMM Raw .156 ± .049 .103 ± .032 .597 ± .027
  ${\text{COB}}{{{\text{E}}}^{2}}$ .151 ± .060 .101 ± .038 .597 ± .031
  ${\text{GIC}}{{{\text{A}}}^{2}}$ .196 ± .092 .127 ± .059 .630 ± .147
  ${\text{GLR}}{{{\text{O}}}^{3}}$ .037 ± .059 .018 ± .044 .709 ± .037
  ${\text{GLR}}{{{\text{O}}}^{4}}$ .033 ± .021 .021 ± .021 .720 ± .006
  ${\text{GTL}}{{{\text{D}}}^{2}}$ .052 ± .056 .038 ± .049 .712 ± .030
  ${\text{GTL}}{{{\text{D}}}^{3}}$ .058 ± .060 .034 ± .041 .704 ± .032

Примечание. Формат ячеек таблицы: “$x \pm y$”, где $x$ – среднее значение, $y$ – стандартное отклонение. Для HAC за ${{R}_{{{\text{complete}}}}}$ обозначен пересчет расстояний через взятие максимума, за ${{\rho }_{{{\text{corr}}}}}$, ${{\rho }_{{cos}}}$, ${{\rho }_{{{\text{RBF}}}}}$, ${{\rho }_{{{{l}_{2}}}}}$, ${{\rho }_{{{\text{Canb}}.}}}$ – корреляция, косинусное, RBF, ${{l}_{2}}$ расстояния и расстояние Канберра; ${\text{Ra}}{{{\text{w}}}^{1}}$: для HAC – ${{\rho }_{{{{l}_{2}}}}}$/${{\rho }_{{cos}}}$/${{\rho }_{{{\text{RBF}}}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{COB}}{{{\text{E}}}^{1}}$: ${{r}_{c}} = 3$, ${{\rho }_{{{{l}_{2}}}}}$/${{\rho }_{{cos}}}$/${{\rho }_{{{\text{RBF}}}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{COB}}{{{\text{E}}}^{2}}$: ${{r}_{c}} = 4$, для HAC – ${{\rho }_{{{{l}_{2}}}}}$/${{\rho }_{{cos}}}$/${{\rho }_{{{\text{RBF}}}}}$/${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GIC}}{{{\text{A}}}^{1}}$: ${{r}_{c}} = 2$, ${{r}_{i}} = 1$, ${{\rho }_{{{{l}_{2}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GIC}}{{{\text{A}}}^{2}}$: ${{r}_{c}} = 1$, ${{r}_{i}} = 2$; ${\text{GLR}}{{{\text{O}}}^{1}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 5$, ${{\rho }_{{{\text{Canb}}.}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GLR}}{{{\text{O}}}^{2}}$: ${{r}_{c}} = 3$, ${{r}_{i}} = 5$; ${\text{GLR}}{{{\text{O}}}^{3}}$: ${{r}_{c}} = 3$, ${{r}_{i}} = 4$; ${\text{GLR}}{{{\text{O}}}^{4}}$: ${{r}_{c}} = 1$, ${{r}_{i}} = 3$; ${\text{GTL}}{{{\text{D}}}^{1}}$: ${{r}_{c}} = 1$, ${{r}_{i}} = 3$, ${{\rho }_{{{\text{corr}}}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GTL}}{{{\text{D}}}^{2}}$: ${{r}_{c}} = 4$, ${{r}_{i}} = 3$, для HAC – ${{\rho }_{{{\text{Canb}}.}}}$, ${{R}_{{{\text{complete}}}}}$; ${\text{GTL}}{{{\text{D}}}^{3}}$: ${{r}_{c}} = 2$, ${{r}_{i}} = 1$.

На фиг. 3 приведены примеры общих и индивидуальных частей, полученных разными способами для одного из изображений из данных ETH80.

Фиг. 3.

Пример изображения из набора ETH80 и соответствующие ему индивидуальная и групповая части: (а) оригинал; (б) COBE, ${{r}_{c}} = 2$; (в) GICA, ${{r}_{c}} = 2$, ${{r}_{i}} = 4$; (г) GLRO, ${{r}_{c}} = 4$, ${{r}_{i}} = 3$; (д) GTLD, ${{r}_{c}} = 5$, ${{r}_{i}} = 3$. (б)–(д) – Верхний ряд изображений отвечает групповым частям, нижний – индивидуальным.

Кластеризация с помощью предложенной VBTD модели может быть выполнена через использование как априорного правила ${{\pi }_{k}}({\mathbf{x}},{{\theta }_{{{{\pi }_{k}}}}})$, так и через вычисление апостериорных вероятностей ${{\gamma }_{k}}$. На наборах данных были протестированы различные конфигурации VBTD моделей с фиксированными рангами (равными 3). Индивидуальные слагаемые моделировались в (Lr, 1) формате. Результаты для лучших пяти конфигураций на апостериорном правиле приведены в табл. 3 в смысле AMI и в табл. 4 в смысле ARI. Все эксперименты проведены с ограничением на число итераций ${{N}_{{{\text{it}}}}} = 100$ при фиксированном случайном состоянии.

Таблица 3.  

VBTD, топ-5 конфигураций в смысле апостериорного AMI

Данные Формат слагаемых Шум,
$L_{k}^{{ - 1}}$
Априорное правило Апостериорное правило
(инд.) (груп.) AMI ARI FMI AMI ARI FMI
ETH80 (Lr, 1) изотр. .646 (.619) .335 (.320) .512 (.498) .683 (.640) .412 (.404) .560 (.538)
  (Lr, 1) TT диаг. .585 (.529) .383 (.272) .522 (.486) .650 (.340) .392 (.116) .544 (.378)
  (Lr, 1) Таккер
(ядро)
изотр. .709 (.458) .508 (.257) .616 (.455) .646 (.420) .454 (.175) .567 (.405)
  (Lr, 1) TT изотр. .427 (.382) .169 (.139) .412 (.389) .623 (.496) .395 (.298) .527 (.429)
  (Lr, 1) диаг. .541 (.507) .255 (.235) .484 (.467) .600 (.541) .336 (.255) .531 (.484)
SMNI Таккер
(ядро)
изотр. .032 (–.005) .081 (.001) .641 (.528) .166 (.166) .215 (.215) .708 (.708)
  Таккер
(фактор)
изотр. .047 (.013) .094 (.039) .705 (.697) .139 (.000) .060 (.000) .731 (.731)
  Таккер
(ядро)
диаг. .026 (–.010) .053 (.010) .712 (.710) .128 (.000) .186 (.000) .731 (.731)
  TT диаг. .017 (.000) .046 (.026) .644 (.644) .104 (.091) .161 (.147) .601 (.596)
  (Lr, 1) изотр. .044 (.001) .098 (.020) .669 (.560) .090 (.006) .162 (.004) .660 (.585)

Примечание. Формат ячеек таблицы: “$x(y)$”, где $x$ – максимальное значение индекса, $y$ – после 100 итераций.

Таблица 4.  

VBTD, топ-5 конфигураций в смысле апостериорного ARI

Данные Формат слагаемых Шум,
$L_{k}^{{ - 1}}$
Априорное правило Апостериорное правило
(инд.) (груп.) AMI ARI FMI AMI ARI FMI
ETH80 (Lr, 1) Таккер (ядро) изотр. .709 (.458) .508 (.257) .616 (.455) .646 (.420) .454 (.175) .567 (.405)
  TT диаг. .310 (.238) .188 (.091) .394 (.344) .555 (.463) .433 (.271) .501 (.410)
  (Lr, 1) изотр. .646 (.619) .335 (.320) .512 (.498) .683 (.640) .412 (.404) .560 (.538)
  (Lr, 1) TT изотр. .427 (.382) .169 (.139) .412 (.389) .623 (.496) .395 (.298) .527 (.429)
  (Lr, 1) TT диаг. .585 (.529) .383 (.272) .522 (.486) .650 (.340) .392 (.116) .544 (.378)
SMNI Таккер (ядро) изотр. .032 (–.005) .081 (.001) .641 (.528) .166 (.166) .215 (.215) .708 (.708)
  Таккер (ядро) диаг. .026 (–.010) .053 (.010) .712 (.710) .128 (.000) .186 (.000) .731 (.731)
  (Lr, 1) изотр. .044 (.001) .098 (.020) .669 (.560) .090 (.006) .162 (.004) .660 (.585)
  TT диаг. .017 (.000) .046 (.026) .644 (.644) .104 (.091) .161 (.147) .601 (.596)
  (Lr, 1) (Lr, 1) диаг. .035 (.000) .074 (–.021) .720 (.689) .071 (.057) .136 (.117) .731 (.640)

Примечание. Формат ячеек таблицы: “$x(y)$”, где $x$ – максимальное значение индекса, $y$ – после 100 итераций.

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

Фиг. 4.

Примеры источников, построенных с помощью VBTD. Конфигурация модели: (Lr, 1) формат индивидуальных слагаемых (ранг 3, $P = 2$), TT формат для общей части ($r = (1,3,3,3)$), диагональная ковариационная матрица шума. Столбцы отвечают кластерам (последний из них – общей части), строки – отдельным восстановленным источникам (без упорядочивания).

8. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Если ориентироваться на индекс FMI, все алгоритмы показывают достаточно оптимистичные результаты. Однако индекс не является устойчивым к предсказаниям с константными значениями, что делает его особенно ненадежным в случае несбалансированных наборов данных с малым числом кластеров (SMNI EEG). В то же время показатели AMI и ARI отражают качество кластеризации более надежным образом. Следует отметить особенность набора данных SMNI EEG: в целом, все использованные алгоритмы показали невысокое качество кластеризации; предполагаем, что результаты можно улучшить с помощью дополнительной предварительной обработки данных, которая выделит более тонкие различия между группами субъектов, но этот вопрос оставлен за рамками работы.

GLRO и GTLD модели показали несколько худшие результаты в сравнении с другими рассмотренными методами, что может быть объяснено особенной структурой контрастированных данных. В случае GICA или COBE предположение о малоранговости применялось только к общей части разложения, которая затем вычиталась из данных. В моделях GLRO и GTLD малоранговое представление строится также и для индивидуальных частей. Предварительно настроенные VBTD модели показали качество кластеризации, сравнимое с GICA и COBE. Вместе с тем в текущей реализации модели мы столкнулись с нестабильностями процесса обучения, которые в ряде случаев проявляют себя заметным зазором между лучшими значениями показателей и значениями на последней итерации. Хотя рассмотренные матричные модели имеют высокие показатели качества, данные модели опираются на избыточное представление данных. В отличие от них, предложенная VBTD модель не требует вычисления параметров для каждого отдельного образца, что делает ее предпочтительнее с точки зрения затрат на хранение параметров. Кроме того, в модели может быть учтено масштабирование по числу размерностей (например, с помощью формата квантизированного тензорного поезда QTT, см. [53]). Используемый метод для вычисления параметров, стохастический вариационный вывод через оптимизацию вариационной нижней грани правдоподобия, поддерживает обучение по подвыборкам, что делает модель масштабируемой в смысле числа данных.

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

9. ВЫВОДЫ И ДАЛЬНЕЙШАЯ РАБОТА

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

Как и матричные модели GICA и COBE, предложенные модели на основе BTD показали способность к приближению групповой и индивидуальной активности, что было проверено методом кластеризации контрастированных (т.е. сохранивших только индивидуальную часть) данных. Заметим, что результаты кластеризации зависят не только от метода контрастирования, но и от метода кластеризации. Вместе с тем, если общие и индивидуальные компоненты разложения потребуется использовать в смысле, отличном от признакового представления для кластеризации, то для выбора модели потребуется использование иного способа тестирования.

Одним из недостатков предложенных моделей GLRO и GTLD является большая вычислительная сложность, которая может быть компенсирована возможностью выделения факторов, специфичных для разных мод, в случае малого числа образцов. Для выделения зависимостей из больших наборов данных более перспективным выбором выглядит VBTD модель.

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

ПРИЛОЖЕНИЕ

1. Использованные в работе обозначения

В настоящей статье использовались следующие обозначения.

Определение. Многомерный массив с вещественными элементами $X \in {{\mathbb{R}}^{{{{n}_{1}} \times \ldots \times {{n}_{d}}}}}$ называется (вещественным) тензором. Элемент тензора $X$, отвечающий мульти-индексу $({{i}_{1}}, \ldots ,{{i}_{d}})$, обозначается как $X[{{i}_{1}}, \ldots ,{{i}_{d}}]$. Его $k$-я ось называется также модой размера ${{n}_{k}}$. Любой тензор может быть записан как вектор-столбец с помощью операции векторизации. Соответствие между индексами оригинала и векторизированного результата задается правилом

$f({{i}_{1}}, \ldots ,{{i}_{d}}) = 1 + \sum\nolimits_{k = 1}^d \,\left[ {({{i}_{k}} - 1)\prod\nolimits_{l = 1}^{k - 1} \,{{n}_{l}}} \right],\quad {\text{vec}}\left( X \right)[f({{i}_{1}}, \ldots ,{{i}_{d}})] = X[{{i}_{1}}, \ldots ,{{i}_{d}}],\quad {{i}_{k}} = \overline {1,{{n}_{k}}} ,\quad k = \overline {1,d} .$

Для преобразования тензора в вектор-строку использовано схожее обозначение: $\widehat {vec}({\kern 1pt} \cdot {\kern 1pt} )$. Для тензоров размерности $d$ также вводится операция матризации, для чего в случае $d > 2$ объединяют его оси. Операцию можно вводить различными способами, здесь рассмотрим подход с сохранением в неизменном виде одной из осей с номером $k$. Будем называть такую операцию разверткой по индексу $k$: ${\text{unfol}}{{{\text{d}}}_{k}}({{X}^{{{{n}_{1}} \times \ldots \times {{n}_{k}} \times \ldots {{n}_{d}}}}}) = {{X}_{{(k)}}} = {{\hat {X}}^{{{{n}_{k}} \times \prod\nolimits_{l \ne k} {{n}_{l}}}}}$. Кроме того, если у тензора есть мода с некоторым номером $l$, размер ${{n}_{l}}$ которой может быть разложен на множители, ${{n}_{l}} = \prod\nolimits_{p = 1}^{{{P}_{l}}} \,{{n}_{{l,p}}}$, к такому тензору можно применить операцию тензоризации (искусственного повышения размерности). Использовано обозначение $\operatorname{Tens} ({\kern 1pt} \cdot {\kern 1pt} )$. Операция требует принять соглашение о пересчете индексов, в данной статье применялся стандартный порядок “по столбцам” (так называемый фортран-порядок).

Пусть имеется два тензора $X \in {{\mathbb{R}}^{{n_{1}^{{(x)}} \times \ldots \times {{n}_{k}} \times \ldots n_{{{{d}_{1}}}}^{{(x)}}}}}$, $Y \in {{\mathbb{R}}^{{n_{1}^{{(y)}} \times \ldots \times {{n}_{k}} \times \ldots n_{{{{d}_{2}}}}^{{(y)}}}}}$ размерностями ${{d}_{1}}$ и ${{d}_{2}}$, и пусть для них имеется общая мода с номером $k$ размером ${{n}_{k}}$. В этом случае определим операцию умножения по $k$-й моде, $X{{ \times }_{k}}Y = Z$, между $X$ и $Y$ по следующему поэлементному правилу:

$Z{\kern 1pt} [ \ldots ,{{i}_{{k - 1}}},{{i}_{{k + 1}}}, \ldots ,{{j}_{{k - 1}}},{{j}_{{k + 1}}}, \ldots ] = \sum\nolimits_{l = 1}^{{{n}_{k}}} \,X[ \ldots ,{{i}_{{k - 1}}},l,{{i}_{{k + 1}}}, \ldots ]Y[ \ldots ,{{j}_{{k - 1}}},l,{{j}_{{k + 1}}}, \ldots ],$

где ${{i}_{p}} = 1, \ldots ,n_{p}^{{(x)}}$, ${{j}_{q}} = 1, \ldots ,n_{q}^{{(y)}}$, $p = \overline {1,{{d}_{1}}} $, $q = \overline {1,{{d}_{2}}} $. В литературе данную операцию также называют свертыванием (см. [2]).

Фробениусова норма для тензора $X$ определяется выражением вида $\left\| X \right\|_{F}^{2} = \left\langle {\operatorname{vec} (X),\operatorname{vec} (X)} \right\rangle $, где $\left\langle { \cdot , \cdot } \right\rangle $ – скалярное (внутреннее) произведение между двумя векторами. Операцию можно обобщить на случай тензоров без дополнительной их векторизации: результат эквивалентен последовательному свертыванию аргументов-тензоров одинакового размера по всем их осям. Если же имеется некоторое подмножество осей, по которым перемножение выполнять не следует, то будем обозначать такую ситуацию следующим образом, используя скобки внутреннего произведения: ${{\left\langle {X,Y} \right\rangle }_{{\{ \alpha \in \Omega \} }}} = X{{ \times }_{{k \in \{ 1, \ldots ,d\} {{\backslash }}\Omega }}}Y$. Иными словами, подстрочные скобки для операции означают, что операция выполняется для всех объектов (по всем индексам), кроме указанных внутри этих скобок. В дополнение к внутреннему, внешнее произведение $X \circ Y = Z$ между двумя тензорами $X$, $Y$ определяется как $Z[{{i}_{1}}, \ldots ,{{i}_{{{{d}_{1}}}}},{{j}_{1}}, \ldots ,{{j}_{{{{d}_{2}}}}}] = X[{{i}_{1}}, \ldots ,{{i}_{{{{d}_{1}}}}}] \cdot Y[{{j}_{1}}, \ldots ,{{j}_{{{{d}_{2}}}}}]$.

Отметим еще два важных произведения. Кронекерово произведение двух матриц $A$ и $B$ размерами $({{n}_{1}},{{n}_{2}})$ и $({{m}_{1}},{{m}_{2}})$ определяется как $\left( {A \otimes B} \right)[{{m}_{1}}(i - 1) + k,{{m}_{2}}(j - 1) + l] = A[i,j]B[k,l]$. Если обе матрицы $A$, $B$ имеют одинаковое количество столбцов (т.е. выполнено ${{n}_{2}} = {{m}_{2}}$), то для них определено произведение Хатри–Рао, которое удобно использовать в контексте развертки тензоров, представленных в каноническом формате: $\left( {A \odot B} \right)[{{m}_{1}}(i - 1) + k,j] = A[i,j]B[k,j]$.

Для удобства ряд обозначений, использованных в работе, приведен в табл. 5. Более подробный разбор базовых операций и тензорных разложений можно найти в [1].

Таблица 5.  

Использованные в работе обозначения

$\boxed1$ Массив из единиц ${{X}_{{(k)}}}$ Матризация тензора $X$ по оси $k$
$i = \overline {1,I} $ $i \in \{ 1, \ldots ,I\} $ ${{A}^{ + }}$ Псевдообращение Мура–Пенроуза
${{( \ldots )}_{{\{ k\} }}}$ Операция без $k$-го составляющего $X{{ \times }_{k}}Y$ Перемножение тензоров по $k$-й оси
$\operatorname{diag} (A)$ Главная диагональ матрицы $\operatorname{diag} ({\mathbf{v}})$ Диагональная матрица, ${\mathbf{v}}$ – диагональ
${\text{off}}\;\operatorname{diag} (A)$ Матрица $A$ с нулями на диагонали $\operatorname{Diag} ({\mathbf{p}})$ Диагональный тензор, ${\mathbf{p}}$ – диагональ
$\widehat {{\text{vec}}}{\kern 1pt} (X)$ Векторизация тензора $X$ (строка) $\operatorname{vec} (X)$ Векторизация тензора $X$ (столбец)
$\left\langle {X,Y} \right\rangle $ Внутреннее произведение $X \circ Y$ Внешнее произведение двух тензоров
$\left| {\left[ {{{C}_{1}}, \ldots ,{{C}_{d}}} \right]} \right|$ Канонический формат $A \odot B$ Матричное произведение Хатри–Рао
$\left| {\left[ {G;{{A}_{1}}, \ldots ,{{A}_{d}}} \right]} \right|$ Формат Таккера $A \otimes B$ Матричное кронекерово произведение

2. Структура матрицы Гессе для BTD разложения с тремя типами слагаемых

Пусть требуется приблизить данный тензор $T \in {{\mathbb{R}}^{{{{n}_{1}} \times \ldots \times {{n}_{d}}}}}$ в формате (8). Обозначим через ${\mathbf{x}}$ вектор-столбец, составленный из векторизации всех параметров модели:

(35)
${\mathbf{x}} = \mathop {\left[ { \ldots \widehat {\operatorname{vec} }{\kern 1pt} (C_{k}^{{[s]}}) \ldots \widehat {\operatorname{vec} }{\kern 1pt} (A_{k}^{{[m]}}) \ldots \widehat {\operatorname{vec} }{\kern 1pt} ({{G}^{{[m]}}}) \ldots } \right]}\nolimits^{\text{T}} ,$

где $C_{k}^{{[s]}}$ есть $k$-я фактор-матрица для $s$-го слагаемого в каноническом формате, $A_{k}^{{[m]}}$ есть $k$-я фактор-матрица для $m$-го слагаемого в формате Таккера, ${{G}^{{[m]}}}$ – ядро Таккера для $m$-го слагаемого в формате Таккера, и пусть $F = F({\mathbf{x}})$ – восстановленный тензор с помощью текущего приближения параметров:

$F({\mathbf{x}}) = \sum\nolimits_{m = 1}^M \,\left| {\left[ {{{G}^{{[m]}}};A_{1}^{{[m]}}, \ldots ,A_{d}^{{[m]}}} \right]} \right| + \left| {\left[ {{{C}_{1}}, \ldots ,{{C}_{d}}} \right]} \right|$
(как в (8)). Для дальнейших целей введем дополнительное обозначение $Z({\mathbf{x}}) = F({\mathbf{x}}) - T$. Соответствующая оптимизационная задача может быть сформулирована в виде задачи нелинейных наименьших квадратов, градиент и матрица Гессе которой имеют специальный вид:
(36)
$\begin{gathered} \mathop {min}\limits_{\mathbf{x}} f({\mathbf{x}}) = \frac{1}{2}\mathop {\left. {\parallel F({\mathbf{x}}) - T} \right\|}\nolimits_F^2 ,\quad H[f]({\mathbf{x}}) = {{J}^{{\text{T}}}}({\mathbf{x}})J({\mathbf{x}}) + Q({\mathbf{x}}), \\ \nabla f({\mathbf{x}}) = {{J}^{{\text{T}}}}({\mathbf{x}}){\text{vec}}(F({\mathbf{x}}) - T),\quad Q({\mathbf{x}}) = \sum\limits_i {\operatorname{vec} (Z({\mathbf{x}}))[i] \cdot {{\nabla }^{2}}(\operatorname{vec} (Z({\mathbf{x}}))[i])} . \\ \end{gathered} $
Здесь $J({\mathbf{x}})$ – матрица Якоби для $f({\mathbf{x}})$:
(37)
$J({\mathbf{x}}) = \left[ {(V_{{\{ k\} }}^{{[s]}} \otimes {{I}_{{{{n}_{k}}}}}) \ldots (G_{{(k)}}^{{[m]}}V_{{\{ k\} }}^{{[m]{\text{T}}}} \otimes {{I}_{{{{n}_{k}}}}}) \ldots {{V}^{{[m]}}} \ldots } \right],$
$k = \overline {1,d} ,$ $m = \overline {1,M} ,$ $s = \overline {1,S} $ и $V_{{\{ k\} }}^{{[s]}} = C_{d}^{{[s]}} \odot \ldots \odot C_{{k + 1}}^{{[s]}} \odot C_{{k - 1}}^{{[s]}} \odot \ldots \odot C_{1}^{{[s]}}$, $V_{{\{ k\} }}^{{[m]}} = A_{d}^{{[m]}} \otimes \ldots \otimes A_{{k + 1}}^{{[m]}} \otimes $ $ \otimes \;A_{{k - 1}}^{{[m]}} \otimes \ldots \otimes A_{1}^{{[m]}}$, ${{V}^{{[m]}}} = A_{d}^{{[m]}} \otimes \ldots \otimes A_{1}^{{[m]}}$. Главная часть матрицы Гессе имеет следующую структуру:
(38)
${{J}^{{\text{T}}}}J = \left[ {\begin{array}{*{20}{c}} {{{{({\text{G}}{{{\text{r}}}^{{{\text{CC}}}}})}}_{{{{k}_{1}},{{k}_{2}};{{s}_{1}},{{s}_{2}}; \cdot , \cdot }}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{CA}}}}})}}_{{{{k}_{1}},{{k}_{2}};{{s}_{1}}, \cdot ; \cdot ,{{m}_{2}}}}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{CG}}}}})}}_{{{{k}_{1}}, \cdot ;{{s}_{1}}, \cdot ; \cdot ,{{m}_{2}}}}}} \\ {{{{({\text{G}}{{{\text{r}}}^{{{\text{AC}}}}})}}_{{{{k}_{1}},{{k}_{2}}; \cdot ,{{s}_{2}};{{m}_{1}}, \cdot }}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{AA}}}}})}}_{{{{k}_{1}},{{k}_{2}}; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{AG}}}}})}}_{{{{k}_{1}}, \cdot ; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}}} \\ {{{{({\text{G}}{{{\text{r}}}^{{{\text{GC}}}}})}}_{{{{k}_{1}},{{k}_{2}}; \cdot ,{{s}_{2}};{{m}_{1}}, \cdot }}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{GA}}}}})}}_{{ \cdot ,{{k}_{2}}; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}}}&{{{{({\text{G}}{{{\text{r}}}^{{{\text{GG}}}}})}}_{{ \cdot , \cdot ; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}}} \end{array}} \right],$
где ${{k}_{1}},{{k}_{2}} = \overline {1,d} ,{{m}_{1}},{{m}_{2}} = \overline {1,M} ,{{s}_{1}},{{s}_{2}} = \overline {1,S} $, и соответствующие блоки имеют вид
$\begin{gathered} {\text{Gr}}_{{{{k}_{1}},{{k}_{2}};{{s}_{1}},{{s}_{2}}; \cdot , \cdot }}^{{{\text{CC}}}} = (V_{{\{ {{k}_{1}}\} }}^{{[{{s}_{1}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{1}}}}}}}})P_{{1,{{k}_{2}}}}^{{{{k}_{1}},1}}(V_{{\{ {{k}_{2}}\} }}^{{[{{s}_{2}}]}} \otimes {{I}_{{{{n}_{{{{k}_{2}}}}}}}}), \\ {\text{Gr}}_{{{{k}_{1}},{{k}_{2}};{{s}_{1}}, \cdot ; \cdot ,{{m}_{2}}}}^{{{\text{CA}}}} = (V_{{\{ {{k}_{1}}\} }}^{{[{{s}_{1}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{1}}}}}}}})P_{{1,{{k}_{2}}}}^{{{{k}_{1}},1}}(V_{{\{ {{k}_{2}}\} }}^{{[{{m}_{2}}]}}G_{{({{k}_{2}})}}^{{[{{m}_{2}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{2}}}}}}}}), \\ \end{gathered} $
(39)
$\begin{gathered} {\text{Gr}}_{{{{k}_{1}},{{k}_{2}};{{s}_{1}}, \cdot ; \cdot ,{{m}_{2}}}}^{{{\text{CG}}}} = (V_{{\{ {{k}_{1}}\} }}^{{[{{s}_{1}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{1}}}}}}}})P_{1}^{{{{k}_{1}}}}({{V}^{{[{{m}_{2}}]}}}), \\ {\text{Gr}}_{{{{k}_{1}},{{k}_{2}};{{m}_{1}},{{m}_{2}}}}^{{{\text{AA}}}} = (G_{{({{k}_{1}})}}^{{[{{m}_{1}}]}}V_{{\{ {{k}_{1}}\} }}^{{[{{m}_{1}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{1}}}}}}}})P_{{1,{{k}_{2}}}}^{{{{k}_{1}},1}}(V_{{\{ {{k}_{2}}\} }}^{{[{{m}_{2}}]}}G_{{({{k}_{2}})}}^{{[{{m}_{2}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{2}}}}}}}}), \\ \end{gathered} $
$\begin{gathered} {\text{Gr}}_{{{{k}_{1}}, \cdot ; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}^{{{\text{AG}}}} = (G_{{({{k}_{1}})}}^{{[{{m}_{1}}]}}V_{{\{ {{k}_{1}}\} }}^{{[{{m}_{1}}]{\text{T}}}} \otimes {{I}_{{{{n}_{{{{k}_{1}}}}}}}})P_{1}^{{{{k}_{1}}}}({{V}^{{[{{m}_{2}}]}}}), \\ {\text{Gr}}_{{ \cdot , \cdot ; \cdot , \cdot ;{{m}_{1}},{{m}_{2}}}}^{{{\text{GG}}}} = ({{V}^{{[{{m}_{1}}]{\text{T}}}}})({{V}^{{[{{m}_{2}}]}}}), \\ \end{gathered} $
здесь $P_{ \cdot }^{ \cdot }$ – коммутационные матрицы (“матрицы перестановки осей”), верхние и нижние индексы обозначают исходные и новые номера мод соответственно. Например, векторизация транспонирования матрицы $A$ может быть записана в виде $P_{{2,1}}^{{1,2}}\operatorname{vec} (A) = \operatorname{vec} ({{A}^{{\text{T}}}})$. В тензорном случае все прочие индексы при этом упорядочиваются согласно их естественной нумерации.

На практике слагаемое $Q({\mathbf{x}})$ часто исключается по нескольким причинам: требует дополнительных вычислительных ресурсов, конфликтует с требованием оптимизационного метода на симметричность и положительную определенность приближения гессиана, а также по причине того, что в окрестности локального минимума его влиянием можно пренебречь. Вместе с тем приводим выражения для вычисления данной квадратичной части гессиана, используя тензорное представление:

(40)
${\text{Tens}}\left( {\frac{{{{\partial }^{2}}{\mathbf{z}}}}{{\partial \widehat {{\text{vec}}}{\kern 1pt} (C_{k}^{{[{{s}_{1}}]}})\partial \widehat {{\text{vec}}}{\kern 1pt} (C_{l}^{{[{{s}_{2}}]}})}} \cdot {\mathbf{z}}} \right) = \left\{ \begin{gathered} {{\left\langle {{{{\left. {\left| {\left[ {C_{1}^{{[{{s}_{1}}]}} \ldots ,C_{d}^{{[{{s}_{1}}]}}} \right]} \right|{\kern 1pt} } \right|}}_{{C_{k}^{{[{{s}_{1}}]}},C_{l}^{{[{{s}_{1}}]}} = {{I}_{{{{L}_{{{{s}_{1}}}}}}}}}}},\quad Z} \right\rangle }_{{\{ k,l\} }}},\quad k \ne l,\quad {{s}_{1}} = {{s}_{2}}, \hfill \\ 0,\quad {\text{иначе}}, \hfill \\ \end{gathered} \right.$
(41)
${\text{Tens}}\left( {\tfrac{{{{\partial }^{2}}{\mathbf{z}}}}{{\partial \widehat {{\text{vec}}}{\kern 1pt} (A_{k}^{{[{{m}_{1}}]}})\partial \widehat {{\text{vec}}}{\kern 1pt} (A_{l}^{{[{{m}_{2}}]}})}} \cdot {\mathbf{z}}} \right) = \left\{ \begin{gathered} {{\langle \underbrace {\left| {\left[ {{{G}^{{[{{m}_{1}}]}}};A_{1}^{{[{{m}_{1}}]}}, \ldots ,A_{d}^{{[{{m}_{1}}]}}} \right]} \right|}_{A_{k}^{{[{{m}_{1}}]}} = {{I}_{{r_{k}^{{[{{m}_{1}}]}}}}},A_{l}^{{[{{m}_{1}}]}} = {{I}_{{r_{l}^{{[{{m}_{1}}]}}}}}},Z\rangle }_{{\{ k,l\} }}},\quad k \ne l,\quad {{m}_{1}} = {{m}_{2}}, \hfill \\ 0,\quad {\text{иначе}}, \hfill \\ \end{gathered} \right.$
(42)
${\text{Tens}}\left( {\tfrac{{{{\partial }^{2}}{\mathbf{z}}}}{{\partial \widehat {{\text{vec}}}{\kern 1pt} (A_{k}^{{[{{m}_{1}}]}})\partial \widehat {{\text{vec}}}{\kern 1pt} ({{G}^{{[{{m}_{2}}]}}})}} \cdot {\mathbf{z}}} \right) = \left\{ \begin{gathered} {{\left. {\left| {\left[ {Z;A_{1}^{{[{{m}_{1}}]}}, \ldots ,A_{d}^{{[{{m}_{1}}]}}} \right]} \right|{\kern 1pt} } \right|}_{{A_{k}^{{[{{m}_{1}}]}} = {{I}_{{r_{k}^{{[{{m}_{1}}]}}}}}}}},\quad {{m}_{1}} = {{m}_{2}}, \hfill \\ 0,\quad {{m}_{1}} \ne {{m}_{2}}, \hfill \\ \end{gathered} \right.$
${{m}_{1}},{{m}_{2}} = \overline {1,M} $, $k,l = \overline {1,d} $, ${\mathbf{z}} = \operatorname{vec} (Z({\mathbf{x}}))$, $Z = Z({\mathbf{x}})$, и все остальные вторые производные равны нулю. Важно отметить, что использование явного представления матрицы Гессе требует больших затрат по памяти. Для оптимизации расходования памяти имеет смысл использовать не полное представление, а использовать результат умножения матрицы на вектор, опираясь на структуру параметров, в том или ином итеративном методе, что и было сделано в используемой реализации.

3. Встраивание групповых условий в оптимизационную задачу

В данном приложении описано, как встроить сформулированные условия в оптимизационную задачу. Предложенная детерминистическая модель для группового анализа данных включает в себя условия на фактор-матрицу групповой оси, а также условие отделимости вида (9). Суммируем их для обоих вариантов модели, полной (Lr, 1), GLRO (6), и модели GTLD (7):

$\begin{gathered} {{C}_{{d + 1}}} = \left[ {\begin{array}{*{20}{c}} {{{I}_{N}}}&{\mathbf{p}} \end{array}} \right],\quad \quad \quad \quad \quad \quad \quad \quad {{A}_{{d + 1}}} = \operatorname{diag} ({\mathbf{p}}), \\ {{p}_{1}} + \ldots + {{p}_{N}} = {{p}_{{{\text{cum}}}}},\quad {{p}_{i}} \geqslant {{p}_{{{\text{min}}}}},\quad \quad \quad {{p}_{1}} + \ldots + {{p}_{N}} = {{p}_{{{\text{cum}}}}},\quad {{p}_{i}} \geqslant {{p}_{{{\text{min}}}}}, \\ {{(C_{\gamma }^{{[N + 1]}})}^{{\text{T}}}}C_{\gamma }^{{[i]}} = 0,\quad \gamma \in \Omega ,\quad i = \overline {1,N} ,\quad \quad \quad A_{\gamma }^{{\text{T}}}C_{\gamma }^{{[i]}} = 0,\quad \gamma \in \Omega ,\quad i = \overline {1,N} , \\ \end{gathered} $
где ${{p}_{{{\text{min}}}}}$ и ${{p}_{{{\text{cum}}}}}$ – фиксированные константы, $\Omega $ – набор мод для условия отделимости, $N$ – общее число образцов. Существуют разные способы учета этих условий в задаче (36), обратимся к двум методам: проекционному и методу множителей Лагранжа. В проекционном методе каждая итерация состоит из шага обновления решения безусловной задачи оптимизации с последующим его проектированием на область допустимости, заданной через условия. Первая часть проектора, отвечающая  условию отделимости, является проектированием столбцов индивидуальных фактор-матриц на пространство, ортогональное столбцам общей фактор-матрицы: $P_{Y}^{ \bot }(X) = (I - Y{{({{Y}^{{\text{T}}}}Y)}^{{ - 1}}}{{Y}^{{\text{T}}}})X$. Оставшаяся часть проектора является проекцией ${\mathbf{p}}$ на внутренность ${{l}_{1}}$ шара с вырезом в окрестности нуля и положительными значениями координат (см. [54]), $\{ {\mathbf{y}} \in {{\mathbb{R}}^{N}}\,{\text{|}}\,{{y}_{1}} + \ldots + {{y}_{N}} = {{p}_{{{\text{cum}}}}},{{y}_{i}} \geqslant {{p}_{{{\text{min}}}}} > 0\} $.

Второй метод состоит в использовании множителей Лагранжа, сводящих условную задачу оптимизации к безусловной. Запишем соответствующее слагаемое целевого функционала $G({\mathbf{x}})$:

(44)
$\begin{gathered} G\left( {\mathbf{x}} \right){\text{ }} = \sum\limits_{\gamma \in \Omega } {\sum\limits_{i = 1}^N {\left\langle {\operatorname{vec} (U_{\gamma }^{{\text{T}}}C_{\gamma }^{{[i]}}),{{{\hat {\mu }}}_{{\gamma ,i}}}} \right\rangle } } - \kappa \cdot \left( {{{p}_{{{\text{cum}}}}} - \left\langle {{\mathbf{p}},\boxed1} \right\rangle } \right) - \left\langle {\zeta ,{\mathbf{p}} - {{p}_{{{\text{min}}}}} \cdot \boxed1} \right\rangle , \\ {{U}_{\gamma }} = \left\{ \begin{gathered} C_{\gamma }^{{[N + 1])}},\quad {\text{для}}\;{\text{модели}}\;(6), \hfill \\ {{A}_{\gamma }},\quad {\text{для}}\;{\text{модели}}\;(7), \hfill \\ \end{gathered} \right.\quad \left\{ \begin{gathered} {{C}_{{d + 1}}} = \left[ {{{I}_{N}}{\mathbf{p}}} \right],\quad {\text{для}}\;{\text{модели}}\;(6) \hfill \\ {{C}_{{d + 1}}} = {{I}_{N}},{{A}_{{d + 1}}} = \operatorname{diag} ({\mathbf{p}}),\quad {\text{для}}\;{\text{модели}}\;(7). \hfill \\ \end{gathered} \right. \\ \end{gathered} $

Число дополнительных параметров можно уменьшить, ослабив условия на ортогональность групповых и индивидуальных фактор-матриц:

(45)
$\begin{array}{*{20}{c}} {G({\mathbf{x}}) = \sum\limits_{\gamma \in \Omega } {\frac{{{{\mu }_{\gamma }}}}{2}} {\kern 1pt} \sum\limits_{i = 1}^N {\left\| {U_{\gamma }^{{\text{T}}}C_{\gamma }^{{[i]}}} \right\|_{F}^{2}} - \kappa \cdot \left( {{{p}_{{{\text{cum}}}}} - \left\langle {{\mathbf{p}},\boxed1} \right\rangle } \right) - \left\langle {\zeta ,{\mathbf{p}} - {{p}_{{{\text{min}}}}}\boxed1} \right\rangle ,} \end{array}$
с теми же ${{U}_{\gamma }}$, что и в (44). Обозначим вектор производных по вспомогательным переменным для $G({\mathbf{x}})$ как $g({\mathbf{x}})$. Градиент для $g({\mathbf{x}})$ в случае модели (7) имеет следующий вид:
(46)
$\begin{gathered} \begin{array}{*{20}{c}} {\nabla g({\mathbf{x}}) = [\begin{array}{*{20}{c}} {\nabla {{g}_{{{{\mu }_{{{{\gamma }_{1}}}}}}}}({\mathbf{x}})}& \ldots &{\nabla {{g}_{{{{\mu }_{{{{\gamma }_{{|\Omega |}}}}}}}}}({\mathbf{x}})}&{\nabla {{g}_{\kappa }}({\mathbf{x}})}&{\nabla {{g}_{{{{\zeta }_{1}}}}}({\mathbf{x}})}& \ldots &{\nabla {{g}_{{{{\zeta }_{N}}}}}({\mathbf{x}})} \end{array}],} \end{array} \\ \nabla {{g}_{{{{\mu }_{{{{\gamma }_{j}}}}}}}}{{({\mathbf{x}})}^{{\text{T}}}} = [\begin{array}{*{20}{c}} \ldots &0&{\widehat {{\text{vec}}}{\kern 1pt} ({{A}_{{{{\gamma }_{j}}}}}A_{{{{\gamma }_{j}}}}^{{\text{T}}}{{C}_{{{{\gamma }_{j}}}}})}&0& \ldots &0&{\widehat {{\text{vec}}}{\kern 1pt} ({{C}_{{{{\gamma }_{j}}}}}C_{{{{\gamma }_{j}}}}^{{\text{T}}}{{A}_{{{{\gamma }_{j}}}}})}&0& \ldots \end{array}{\kern 1pt} ], \\ \nabla {{g}_{\kappa }}{{({\mathbf{x}})}^{{\text{T}}}} = [\begin{array}{*{20}{c}} \ldots &0&{{{{\boxed1}}_{{1 \times N}}}}&0& \ldots \end{array}],\quad \nabla {{g}_{{{{\zeta }_{i}}}}}{{({\mathbf{x}})}^{{\text{T}}}} = \left[ {\begin{array}{*{20}{c}} \ldots &0&{\underbrace { - 1}_{i{\text{ - я}}\;{\text{позиция}}\;{\mathbf{p}}}}&0& \ldots \end{array}} \right]. \\ \end{gathered} $
Наконец, каждая итерация требует решения системы с расширенной матрицей Гессе, $\tau = \mathop {[{{\mu }_{{{{\gamma }_{1}}}}}, \ldots ,{{\mu }_{{{{\gamma }_{{|\Omega |}}}}}},\kappa ,{{\zeta }_{1}}, \ldots ,{{\zeta }_{N}}]}\nolimits^{\text{T}} $:

(47)
$\left[ {\begin{array}{*{20}{c}} {H[f]({\mathbf{x}})}&{\nabla g({\mathbf{x}})} \\ {\nabla g{{{({\mathbf{x}})}}^{{\text{T}}}}}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta {\mathbf{x}}} \\ {\Delta \tau } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - \nabla f({\mathbf{x}})} \\ { - g({\mathbf{x}})} \end{array}} \right].$

Список литературы

  1. Kolda T.G., Bader B.W. Tensor decompositions and applications // SIAM Rev. 2009. V. 51 Iss. 3. P. 455–500.

  2. Cichocki A., Lee N., Oseledets I., Phan A.-H., Zhao Q., Mandic D.P. Tensor networks for dimensionality reduction and large-scale optimization: Part 1, low-rank tensor decompositions // Found. Trends Mach. Learn. 2016. V. 9. Iss. 4–5. P. 249–429.

  3. Sidiropoulos N.D., De Lathauwer L., Fu X., Huang K., Papalexakis E.E., Faloutsos C. Tensor decomposition for signal processing and machine learning // IEEE Trans. Signal Proc. 2017. V. 65. Iss. 13. P. 3551–3582.

  4. Phan A.-H., Cichocki A. Tensor decompositions for feature extraction and classification of high dimensional datasets // IEICE Nonlin. Theory and its App. 2010. V. 1. Iss. 1. P. 37–68.

  5. De Lathauwer L. Decompositions of a higher-order tensor in block terms – part I: Lemmas for partitioned matrices // SIAM J. Matrix Anal. Appl. 2008. V. 30. Iss. 3. P. 1022–1032.

  6. De Lathauwer L. Decompositions of a higher-order tensor in block terms – part II: Definitions and uniqueness // SIAM J. Matrix Anal. Appl. 2008. V. 30. Iss. 3. P. 1033–1066.

  7. De Lathauwer L., Nion D. Decompositions of a higher-order tensor in block terms – part III: Alternating least squares algorithms // SIAM J. Matrix Anal. Appl. 2008. V. 30. Iss. 3. P. 1067–1083.

  8. Prasad G., Jahanshad N., Aganj I., Lenglet C., Sapiro G., Toga A.W., Thompson P.M. Atlas-based fiber clustering for multi-subject analysis of high angular resolution diffusion imaging tractography // Proc. IEEE Int. Symp. Biomed. Imaging. 2011. P. 276–280.

  9. Calhoun V.D., Adali T. Multisubject independent component analysis of fMRI: a decade of intrinsic networks, default mode, and neurodiagnostic discovery // IEEE Rev. Biomed. Eng. 2012. V. 5. P. 60–73.

  10. Beckmann C.F., Smith S.M. Tensorial extensions of independent component analysis for multisubject fMRI analysis. // Neuroimage. 2005. V. 25. Iss. 1. P. 294–311.

  11. Nazarenko D.V., Kharyuk P.V., Oseledets I.V., Rodin I.A., Shpigun O.A. Machine learning for LC–MS medicinal plants identification // Chemom. Intell. Lab. Syst. 2016. V. 156. P. 174–180.

  12. Xia P., Bai Z., Liang T., Yang D., Liang Z., Yan X., Liu Y. High-performance liquid chromatography based chemical fingerprint analysis and chemometric approaches for the identification and distinction of three endangered panax plants in Southeast Asia // J. Sep. Sci. 2016. V. 39. Iss. 20. P. 3880–3888.

  13. Smilde A., Bro R., Geladi P. Multi-way analysis: applications in the chemical sciences // UK: John Wiley & Sons, 2005.

  14. Cichocki A., Mandic D., De Lathauwer L., Zhou G., Zhao Q., Caiafa C., Phan A.H. Tensor decompositions for signal processing applications: From two-way to multiway component analysis // IEEE Signal Process. Mag. 2015. V. 32. Iss. 2. P. 145–163.

  15. Zhou G., Zhao Q., Zhang Y., Adalı T., Xie Sh., Cichocki A. Linked component analysis from matrices to high-order tensors: Applications to biomedical data // Proc. IEEE. 2016. V. 104. Iss. 2. P. 310–331.

  16. Харюк П.В. Групповой анализ данных на основе блочного канонического разложения // Сб. тезисов 59-й научной конференции МФТИ. M.: МФТИ, 2016.

  17. Харюк П.В. Классификация сигналов с помощью блочного тензорного разложения в задаче группового анализа данных // Сб. тезисов XXIV Международной научной конференции Ломоносов-2017. М.: ООО “МАКС Пресс”, 2017. С. 152–153.

  18. Calhoun V.D., Liu J., Adalı T. A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data. Neuroimage, 2009. V. 45. Iss. 1. P. S163–S172.

  19. Zhou G., Cichocki A., Zhang Y., Mandic D.P. Group component analysis for multiblock data: Common and individual feature extraction. IEEE Trans. Neural Netw. Learn. Syst. 2016. V. 27. Iss. 11. P. 2426–2439.

  20. Lock E.F., Hoadley K.A., Marron J.S., Nobel A.B. Joint and individual variation explained (JIVE) for integrated analysis of multiple data types // Ann. Appl. Stat. 2013. V. 7. Iss. 1. P. 523.

  21. Sorber L., Van Barel M., De Lathauwer L. Optimization-based algorithms for tensor decompositions: Canonical polyadic decomposition, decomposition in rank-(Lr, Lr, 1) terms, and a new generalization // SIAM J. Optim. 2013. V. 23. Iss. 2. P. 695–720.

  22. Tipping M.E., Bishop C.M. Mixtures of probabilistic principal component analyzers // Neural Comput. 1999. V. 11. Iss. 2. P. 443–482.

  23. Tipping M.E., Bishop C.M. Probabilistic principal component analysis // J. R. Stat. Soc. Series B Stat. Methodol. 1999. V. 61. Iss. 3. P. 611–622.

  24. Bishop C.M. Pattern recognition and machine learning // NY: Springer, 2006.

  25. Oseledets I.V. Tensor-train decomposition // SIAM J. Sci. Comput. 2011. V. 33. Iss. 5. P. 2295–2317.

  26. Kingma D.P., Welling M. Auto-encoding variational bayes // arXiv:1312.6114, 2013.

  27. Kharyuk P., Nazarenko D., Oseledets I., Rodin I., Shpigun O., Tsitsilin A., Lavrentyev M. Employing fingerprinting of medicinal plants by means of LC-MS and machine learning for species identification task // Sci. Rep. 2018. V. 8. Iss. 1. N. 17053.

  28. Björck A., Golub G.H. Numerical methods for computing angles between linear subspaces // Math. of Computat. 1973. V. 27. Iss. 123. P. 579–594.

  29. Vervliet N., Debals O., Sorber L., Van Barel M., De Lathauwer L. Tensorlab user guide // 2016.

  30. Blei D., Ranganath R., Mohamed S. Variational inference: Foundations and modern methods // NIPS Tutorial, 2016.

  31. Wingate D., Weber T. Automated variational inference in probabilistic programming arXiv:1301.1299, 2013.

  32. Bingham E., Chen J.P., Jankowiak M., Obermeyer F., Pradhan N., Karaletsos T., Singh R., Szerlip P., Horsfall P., Goodman N.D. Pyro: Deep universal probabilistic programming // J. Mach. Learn. Res. 2019. V. 20. Iss. 1. P. 973–978.

  33. Paszke A., Gross S., Chintala S., Chanan G., Yang E., DeVito Z., Lin Z., Desmaison A., Antiga L., Lerer A. Automatic differentiation in PyTorch // NIPS Autodiff Workshop, 2017.

  34. Calhoun V.D., Adali T., Pearlson G.D., Pekar J.J. A method for making group inferences from functional MRI data using independent component analysis // Hum. Brain Mapp. 2001. V. 14. Iss. 3. P. 140–151.

  35. Guo Y., Pagnoni G. A unified framework for group independent component analysis for multi-subject fMRI data // NeuroImage. 2008. V. 42. Iss. 3. P. 1078–1093.

  36. Lee H., Choi S. Group nonnegative matrix factorization for EEG classification // Artif. Int. Stat. 2009. P. 320–327.

  37. Yokota T., Cichocki A., Yamashita Y. Linked PARAFAC/CP tensor decomposition and its fast implementation for multi-block tensor analysis // Springer Int. Conf. Neural Inform. Proc. 2012. P. 84–91.

  38. Cichocki A., Zdunek R., Amari S.-I. Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization // Springer Int. Conf. Indep. Comp. Anal. Sig. Sep. 2007. P. 169–176.

  39. Gillis N., Glineur F. Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization // Neural Comput. 2012. V. 24. Iss. 4. P. 1085–1105.

  40. Yokota T., Cichocki A. Linked Tucker2 decomposition for flexible multi-block data analysis // Springer Int. Conf. Neural Inform. Proc. 2014. P. 111–118.

  41. Gong X.F., Lin Q.-H., Debals O., Vervliet N., De Lathauwer L. Coupled rank-(Lm, Ln, 1) block term decomposition by coupled block simultaneous generalized Schur decomposition // IEEE ICASSP. 2016. P. 2554–2558.

  42. Khoromskij B., Khoromskaia V. Low rank tucker-type tensor approximation to classical potentials // Open Math. 2007. V. 5. Iss. 3. P. 523–550.

  43. Continuum Analytics. Anaconda software distribution // [Электронный ресурс] vers. 2-2.4.0, 11.2015. Дата обращения: 1.11.2019.

  44. Oliphant T.E. A guide to NumPy // USA: Trelgol Publ. USA. V. 1. 2006.

  45. Jones E., Oliphant T., Peterson P. SciPy: Open source scientific tools for Python // [Электронный ресурс] 2001. Дата обращения: 1.11.2019.

  46. McKinney W. Data structures for statistical computing in python // Proceed. of the 9th Python in Sci. Conf. V. 2010. 445, P. 51–56.

  47. Pedregosa F., Varoquaux G., Gramfort A., Michel V., Thirion B., Grisel O., Blondel M., Prettenhofer P., Weiss R., Dubourg V., Vanderplas J., Passos A., Cournapeau D., Brucher M., Perrot M., DuchesnayE. Scikit-learn: machine learning in Python // J. of Mach. Learn. Res. 2011. V. 12. P. 2825–2830.

  48. Hunter J.D. Matplotlib: A 2D graphics environment // Comput. Sci. Eng. 2007. V. 9. Iss. 3. P. 90–95.

  49. Waskom M., Botvinnik O., O’Kane D., Hobson P., Lukauskas S., Gemperline D.C., Augspurger T., Halchenko Y., Cole J.B., Warmenhoven J., de Ruiter J., Pye C., Hoyer S., Vanderplas J., Villalba S., Kunter G., Quintero E., Bachant P., Martin M., Meyer K., Miles A., Ram Y., Yarkoni T., Williams M.L., Evans C., Fitzgerald C., Brian, Fonnesbeck C., Lee A., Qalieh A. Seaborn: statistical data visualization // [Электронный ресурс], v.0.8.1, 09.2017. Дата обращения: 1.11.2019.

  50. Kluyver T., Ragan-Kelley B., Pérez F., Granger B., Bussonnier M., Frederic J., Kelley K., Hamrick J., Grout J., Corlay S., Ivanov P., Avila D., Abdalla S., Willing C. Jupyter notebooks – a publishing format for reproducible computational workflows // ELPUB 2016. IOS Press, P. 87–90.

  51. Leibe B., Schiele B. Analyzing appearance and contour based methods for object categorization // IEEE CVPR. 2003. V. 2.

  52. Hastie T., Tibshirani R., Friedman J., Franklin J. The elements of statistical learning: data mining, inference and prediction // Springer Series in Statistics, 2009.

  53. Khoromskij B.N., Oseledets I.V. QTT approximation of elliptic solution operators in higher dimensions // Russ. J. Numer. Anal. Math. Model. 2011. V. 26. Iss. 3. P. 303–322.

  54. Gupta M.D., Kumar S., Xiao J. L1 projections with box constraints. // arXiv:1010.0141, 2010.

Дополнительные материалы отсутствуют.