Журнал вычислительной математики и математической физики, 2022, T. 62, № 11, стр. 1804-1821

Поиск разреженных решений для сверхбольших систем, обладающих тензорной структурой

Д. А. Желтков^1, *, Н. Л. Замарашкин^1, **, С. В. Морозов^1, ***

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

^* E-mail: dmitry.zheltkov@gmail.com
^** E-mail: nikolai.zamarashkin@gmail.com
^*** E-mail: stanis-morozov@yandex.ru

Поступила в редакцию 30.12.2021
После доработки 06.06.2022
Принята к публикации 07.07.2022

EDN: GANDKY
DOI: 10.31857/S0044466922110151

Аннотация

Задача поиска разреженного решения для больших систем линейных уравнений возникает во многих приложениях, связанных с обработкой сигналов. В некоторых случаях размеры возникающих систем оказываются столь велики, что известные методы становятся неэффективными. Решение таких систем возможно только при наличии в них дополнительной структуры. В настоящей работе предлагается эффективный метод поиска разреженных решений сверхбольших систем линейных уравнений, обладающих тензорной структурой определенного вида. Приведенный теоретический анализ и экспериментальные результаты позволяют судить об эффективности предложенного метода. Библ. 14. Фиг. 7.

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

1. ВВЕДЕНИЕ

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

$y = X\beta + \varepsilon ,$

где $\beta \in {{\mathbb{R}}^{n}}$ – сигнал, $y \in {{\mathbb{R}}^{m}}$ – вектор наблюдений, $X \in {{\mathbb{R}}^{{m \times n}}}$ – матрица измерений, а $\varepsilon \in {{\mathbb{R}}^{m}}$ – вектор шума. В некоторых прикладных задачах интерес представляет ситуация, когда размерность сигнала $n$ велика и существенно превосходит число измерений $m$ (т.е. $n \gg m$), а вектор сигнала $\beta $ разреженный (имеет малое число ненулевых компонент).

Обозначим множество ненулевых компонент вектора $\beta \in {{\mathbb{R}}^{n}}$ через

$\operatorname{supp} \beta = \{ i:{{\beta }_{i}} \ne 0\} $

и будем искать $k$-разреженное решение $\beta $ задачи наименьших квадратов

(1)

${{\left\| {y - X\beta } \right\|}_{2}} \to \min ,\quad {{\left\| \beta \right\|}_{0}} = \left| {\operatorname{supp} \beta } \right| \leqslant k.$

Среди методов, используемых для решения (1), наиболее популярным является алгоритм OMP (Orthogonal Matching Pursuit [1]). Сложность OMP оценивается как $\mathcal{O}(mnk).$ Это приемлемо при $n$, сравнимых с ${{10}^{6}}$, но уже недостаточно эффективно при $n$ порядка ${{10}^{9}}$.

В данной работе основное внимание уделяется решению сверхбольших систем линейных уравнений, в которых $n$ имеет порядок ${{10}^{{13}}}$ и более. Очевидно, что при таких размерах линейной системы требуются дополнительные предположения о структуре задачи. В работе таких предположений три:

1) пространство параметров тензоризовано ${{\mathbb{R}}^{n}} = {{\mathbb{R}}^{{{{n}_{1}}}}} \times {{\mathbb{R}}^{{{{n}_{2}}}}} \times \cdots {{\mathbb{R}}^{{{{n}_{d}}}}},$ $n = {{n}_{1}} \times {{n}_{2}} \times \cdots \times {{n}_{d}},$ а действие матрицы $X$ на пространстве ${{\mathbb{R}}^{n}}$ соответствует действию линейного оператора $\mathcal{A}$ на тензоризованном пространстве. Относительно структуры оператора $\mathcal{A}$ будем предполагать, что он допускает возможность эффективного применения к тензорам ранга $1$;

2) нормальное решение задачи наименьших квадратов

$v = \mathop {\arg \min }\limits_u {{\left\| {\mathcal{A}u - y} \right\|}_{2}}$

может быть приближено в ${{\mathbb{R}}^{{{{n}_{1}}}}} \times \cdots \times {{\mathbb{R}}^{{{{n}_{d}}}}}$ тензорами малого тензорного ранга;

3) строки матрицы измерений задают отображение ограниченной изометрии из пространства размерности $m$ в пространство размерности $n$

(2)

$\left( {1 - \delta } \right)\left( {h,l} \right) < \left( {{{\mathcal{A}}^{{\text{т}}}}h,{{\mathcal{A}}^{{\text{т}}}}l} \right) < \left( {1 + \delta } \right)\left( {h,l} \right),\quad {\text{где}}\quad h,l \in {{\mathbb{R}}^{m}}.$

В работе изучается метод, позволяющий при изложенных выше предположениях находить разреженное решение для систем со сверхбольшой размерностью пространства параметров сигнала $n$. Новый метод является адаптацией метода OMP, в которой учитывается структура задачи. По сути, предлагается новый алгоритм выбора существенных (ненулевых) компонент разреженного решения. Реализованная идея выглядит несколько парадоксально. Поиск существенных компонент разреженного решения (1) осуществляется с помощью анализа компонент нормального решения. На первый взгляд может показаться, что в алгоритме происходит замена задачи меньшей вычислительной сложности на задачу большей вычислительной сложности. Это, как будет показано далее, не так, если допустить существование тензорной структуры в нормальном решении и свойства квазиизометрии оператора измерений $\mathcal{A}$.

Оставшаяся часть текста организована следующим образом. В разд. 2 рассмотрен OMP алгоритм поиска разреженных решений. В разд. 3 дается описание предлагаемого метода, который мы называем тензоризованным OMP. В разд. 4 приведен теоретический анализ различных этапов нового метода. Разд. 5 содержит результаты численных экспериментов, сравнение работы классического и тензоризованного алгоритмов, а также ряд дополнительных эвристик, позволяющих повысить качество работы предлагаемого метода. Разд. 6 содержит заключительные замечания.

2. АЛГОРИТМ OMP

Одним из наиболее эффективных методов поиска разреженных решений для задачи наименьших квадратов является Orthogonal Matching Pursuit (OMP) алгоритм [1]. Это жадная итерационная процедура [2], получившая значительное внимание благодаря простоте записи алгоритма и его эффективности в большом числе реальных приложений [3]–[6]. На каждом шаге OMP увеличивает список выбранных компонент на одну так, чтобы добиться наибольшего падения невязки. Стандарное описание OMP метода приводится ниже в алгоритме 1.

Алгоритм 1. Orthogonal Matching Pursuit

Вход: матрица $A \in {{\mathbb{R}}^{{m \times n}}}$, вектор правой части $b \in {{\mathbb{R}}^{m}}$, порог ошибки ${{\varepsilon }_{0}}$.

Инициализация:

• $k = 0$;

• ${{x}^{0}} = 0;$

• ${{r}^{0}} = b - A{{x}^{0}} = b$;

• множество выбранных столбцов ${{S}^{0}} = \not {0};$

while ${{\left\| {{{r}^{k}}} \right\|}_{2}} \geqslant {{\varepsilon }_{0}}$ do

1. Вычислить ошибки для каждого из столбцов $\epsilon (j) = \mathop {\min }\limits_{{{z}_{j}}} \left\| {{{a}_{j}}{{z}_{j}} - {{r}^{k}}} \right\|_{2}^{2}$ для всех $j$. Решением задачи минимизации является $z_{j}^{*} = \frac{{a_{j}^{{\text{т}}}{{r}^{k}}}}{{\left\| {{{a}_{j}}} \right\|_{2}^{2}}}$.

2. Найти столбец ${{j}_{0}}$, сильнее всего уменьшающий невязку, а именно ${{j}_{0}}:\forall j \notin {{S}^{k}},$ $\epsilon ({{j}_{0}}) \leqslant \epsilon (j)$. Затем обновляем множество выбранных элементов: ${{S}^{{k + 1}}} = {{S}^{k}} \cup \{ {{j}_{0}}\} $.

3. Найти оптимальные коэффициенты при столбцах из ${{S}^{{k + 1}}}$, т.е. вычислить ${{x}^{{k + 1}}}$, минимизирующий $\left\| {Ax - b} \right\|_{2}^{2}$ при условии ${\text{supp}}\,{{x}^{{k + 1}}} = {{S}^{{k + 1}}}$.

4. Вычислить невязку ${{r}^{{k + 1}}} = b - A{{x}^{{k + 1}}}$.

5. $k = k + 1$.

end while

return ${{S}^{k}}$

Сложность итерации OMP оценивается через $\mathcal{O}(mn)$ арифметических операций. При размерах $n,$ по порядку больших ${{10}^{9}},$ алгоритм становится неэффективным. Цель данной работы адаптировать OMP алгоритм на случай больших $n$, при условии, что пространство ${{\mathbb{R}}^{n}}$ тензоризовано, а в операторе $\mathcal{A}$ системы и нормальном решении задачи наименьших квадратов присутствует тензорная структура.

3. ОБЩАЯ СХЕМА ТЕНЗОРИЗОВАННОГО OMP

Рассмотрим линейный оператор

$\mathcal{A}:{{\mathbb{R}}^{{{{n}_{1}}}}} \times {{\mathbb{R}}^{{{{n}_{2}}}}} \times \ldots \times {{\mathbb{R}}^{{{{n}_{d}}}}} \to {{\mathbb{R}}^{m}}.$

Будем считать, что столбцы матрицы оператора $\mathcal{A}$ имеют одинаковую длину. От оператора потребуем возможности быстрого умножения на векторы, представимые тензорами ранга $1,$ т.е. на векторы вида ${{v}_{1}} \otimes {{v}_{2}} \otimes \ldots \otimes {{v}_{d}}$, где ${{v}_{i}} \in {{\mathbb{R}}^{{{{n}_{i}}}}}$ (здесь $ \otimes $ обозначает кронекерово произведение).

Кроме того, предположим, что нормальное решение задачи наименьших квадратов

$v = \mathop {\arg \min }\limits_u {{\left\| {\mathcal{A}u - y} \right\|}_{2}}$

приближается (возможно довольно грубо) тензором малого ранга. Воспользуемся этим предположением, чтобы снизить сложность стандартного OMP.

Начнем со следующего простого наблюдения: наиболее трудоемкий этап OMP относится к вычислению невязок

(3)

$\epsilon (j) = \mathop {\min }\limits_{{{z}_{j}}} {{\left\| {{{a}_{j}}{{z}_{j}} - r} \right\|}_{2}}$

для всех столбцов $\mathcal{A}$. Действительно, для $j$-го столбца оптимальное значение $z_{j}^{*}$ дается выражением $z_{j}^{*} = \frac{{a_{j}^{{\text{т}}}r}}{{\left\| {{{a}_{j}}} \right\|_{2}^{2}}}$, где $r$ обозначает текущий вектор невязки. Прямое вычисление $z_{j}^{*}$ имеет сложность $\mathcal{O}(m).$ Всего таких вычислений $n.$ Таким образом, только прямое вычисление всех $z_{j}^{*}$ уже имеет сложность $\mathcal{O}(mn)$.

Поскольку ${{a}_{j}}z_{j}^{*}$ и ${{a}_{j}}z_{j}^{*} - r$ ортогональны, то решение задачи (3) эквивалентно поиску максимума среди величин

${{c}_{j}} = \frac{{({{a}_{j}},y)}}{{{{{\left\| {{{a}_{j}}} \right\|}}_{2}}}}$

с последующим выбором столбца с максимальным ${{c}_{j}}$. Учитывая, что все длины столбцов ${{\left\| {{{a}_{j}}} \right\|}_{2}}$ совпадают, наилучшая компонента на текущей итерации OMP дается максимумом модуля скалярного произведения ${\text{|}}{\kern 1pt} ({{a}_{j}},y){\kern 1pt} {\text{|}}$. Таким образом, снизить алгоритмическую сложность OMP возможно, предъявив алгоритм малой сложности для вычисления величин ${\text{|}}{\kern 1pt} ({{a}_{j}},y){\kern 1pt} {\text{|}}$.

Сделаем элементарное, но важное наблюдение. Для оператора ${{\mathcal{A}}^{{\text{т}}}},$ обладающего свойством ограниченной изометрии (см. предположение 3), оценки на величины ${\text{|}}{\kern 1pt} ({{a}_{j}},y){\kern 1pt} {\text{|}}$ можно получить, рассматривая нормальное решение задачи наименьших квадратов. Действительно, нормальное решение $\tilde {u}$ представляется в виде

(4)

$\tilde {u} = {{\mathcal{A}}^{\dag }}y = {{\mathcal{A}}^{{\text{т}}}}{{\left( {\mathcal{A}{{\mathcal{A}}^{{\text{т}}}}} \right)}^{{ - 1}}}y.$

Тогда для компоненты ${{\tilde {u}}_{i}}$

$\left| {{{{\tilde {u}}}_{i}}} \right| = \left| {a_{i}^{{\text{т}}}{{{(\mathcal{A}{{\mathcal{A}}^{{\text{т}}}})}}^{{ - 1}}}y} \right| = \left| {\left( {{{a}_{i}},(\mathcal{A}{{\mathcal{A}}^{{\text{т}}}}{{)}^{{ - 1}}}y} \right)} \right| \leqslant \frac{1}{{1 - \delta }}\left| {\left( {{{\mathcal{A}}^{{\text{т}}}}{{a}_{i}},{{\mathcal{A}}^{{\text{т}}}}{{{(\mathcal{A}{{\mathcal{A}}^{{\text{т}}}})}}^{{ - 1}}}y} \right)} \right| = \frac{1}{{1 - \delta }}\left| {\left( {{{a}_{i}},y} \right)} \right|.$

Аналогично получается оценка снизу

$\left| {{{{\tilde {u}}}_{i}}} \right| \geqslant \frac{1}{{1 + \delta }}\left| {\left( {{{a}_{i}},y} \right)} \right|.$

Имея в виду полученные оценки, будем выбирать существенные компоненты в OMP на основе модуля компонент нормального решения задачи наименьших квадратов.

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

Остается еще вопрос о быстром поиске больших компонент в полученном приближении к нормальному решению. Для тензорных представлений ответ хорошо известен, а соответствующие алгоритмы описаны, например в [7], [8]. Однако в наших численных экспериментах тензорный ранг $T$ всегда выбирался равным $1$. В этом случае выбор $C$ максимальных элементов тензора может быть еще упрощен. А именно, пусть после решения задачи наименьших квадратов для тензора ${{v}_{T}}$ получается представление в виде

${{v}_{T}} = {{u}_{1}} \otimes {{u}_{2}} \otimes \ldots \otimes {{u}_{d}}.$

Выберем из ${{u}_{1}}$ не более $C$ максимальных по модулю элементов и обозначим получившееся множество через ${{G}_{1}}$. Затем перемножим попарно все элементы из ${{G}_{1}}$ с элементами из ${{u}_{2}}$ и выберем не более $C$ максимальных по модулю элементов, а полученное множество обозначим через ${{G}_{2}}$. Продолжая процедуру, в конце получим множество ${{G}_{d}}$, содержащее $C$ максимальных по модулю элементов тензора ${{v}_{T}}$. Сложность такого поиска при простейшей реализации составляет $\mathcal{O}(C({{n}_{1}} + \ldots + {{n}_{d}} + d\log C))$.

Основываясь на сказанном, предлагается следующий порядок вычислений в тензоризованном алгоритме 2.

Algorithm 2. Тензоризованный OMP

Вход: тензоризованный линейный оператор $\mathcal{A}$, вектор правой части $b$, количество компонент разреженного решения $K$, ранг тензора для малоранговой аппроксимации $T$, количество кандидатов $C$.

Инициализация:

• ${{x}^{0}} = 0;$

• ${{r}^{0}} = b - \mathcal{A}{{x}^{0}} = b$;

• множество выбранных столбцов ${{S}^{0}} = \not {0};$

• множество столбцов для ортогонализации ${{Q}^{0}} = \not {0};$

for $k = 0, \ldots ,K - 1$ do

1. Вычислить решение ранга $T$ для решения задачи наименьших квадратов

Для простоты мы обозначили ${{Q}^{k}}$ и множество столбцов, и матрицу ортогонализации к этим столбцам.

2. Найти $C$ элементов с максимальными по модулю элементами в тензоре ${{v}_{T}}$. Обозначим полученное множество $\hat {S}$.

3. Построим матрицу $\hat {A}$, состоящую из столбцов с номерами из $\hat {S}$; $i$-й столбец может быть получен действием оператора $\mathcal{A}$ на единичный вектор.

4. Применим один шаг метода OMP к задаче с матрицей ${{Q}^{k}}\hat {A}$ и правой частью ${{Q}^{k}}b$. В результате получим некоторый номер столбца ${{j}_{0}}$ для матрицы исходного оператора $\mathcal{A}$.

5. Добавим столбец ${{j}_{0}}$ в множество выбранных столбцов: ${{S}^{{k + 1}}} = {{S}^{k}} \cup \{ {{j}_{0}}\} $, а также столбец ${{Q}^{k}}{{a}_{{{{j}_{0}}}}}$ в множество столбцов для ортогонализации: ${{Q}^{{k + 1}}} = {{Q}^{k}} \cup \{ {{Q}^{k}}{{a}_{{{{j}_{0}}}}}\} $.

end for

return ${{S}^{K}}$

Остается вопрос об эффективном поиске решения малого тензорного ранга для ортогонализованной задачи наименьших квадратов

${{v}_{T}} = \mathop {\arg {\text{min}}}\limits_{_{{\substack{ u \in {{\mathbb{R}}^{{{{n}_{1}} \times \ldots \times {{n}_{d}}}}}: \\ \operatorname{rank} u\, \leqslant \,T } }}} {{\left\| {{{Q}^{k}}\mathcal{A}u - {{Q}^{k}}b} \right\|}_{2}}.$

Для этого предлагается использовать простую модификацию алгоритма ALS. Действительно, для решаемой задачи

(5)

${{\left\| {Q\mathcal{A}\left( {\sum\limits_{t = 1}^{\text{т}} \,{{u}_{{1t}}} \otimes {{u}_{{2t}}} \otimes \ldots \otimes {{u}_{{dt}}}} \right) - Qb} \right\|}_{2}} \to \min .$

выберем некоторое $i$ от $1$ до $d$ и фиксируем все ${{u}_{{jt}}}$ при $j \ne i.$ Решим задачу относительно ${{u}_{{it}}}$. Пусть ${{U}_{i}} = {{\left[ {\begin{array}{*{20}{c}} {u_{{i1}}^{{\text{т}}}}&{u_{{i2}}^{{\text{т}}}}& \ldots &{u_{{iT}}^{{\text{т}}}} \end{array}} \right]}^{{\text{т}}}}$, перепишем (5) в виде

${{\left\| {Q{{M}_{i}}{{U}_{i}} - Qb} \right\|}_{2}} \to \min ,$

где ${{M}_{i}} \in {{\mathbb{R}}^{{m \times T{{n}_{i}}}}}$ размерность новой задачи наименьших квадратов. Для получения решения алгоритм ALS многократно итерирует по всем $i$ от $1$ до $d$.

Заметим также, что наличие матрицы ортогонализации $Q$ не существенно влияет на оценку сложности. Действительно, $Q$ имеет вид

$Q = (I - {{w}_{1}}w_{1}^{{\text{т}}})(I - {{w}_{2}}w_{2}^{{\text{т}}}) \ldots (I - {{w}_{k}}w_{k}^{{\text{т}}}),$

где $k$ соответствует количеству векторов, к которым происходит ортогонализация. Таким образом, умножение на $Q$ соответствует последовательности из $k$ одноранговых преобразований и имеет сложность $\mathcal{O}(kmT{{n}_{i}}),$ которая сопоставима со сложностью решения задачи наименьших квадратов $\mathcal{O}(m{{T}^{2}}n_{i}^{2})$. Таким образом, решение ортогонализованной задачи оказывается ненамного сложнее классического алгоритма ALS для нахождения малорангового решения.

Пусть алгоритм совершает ${{N}_{{ALS}}}$ полных итераций (т.е. при которых $i$ пробегает от $1$ до $d$). Тогда полная сложность поиска решения малого ранга равна

$\mathcal{O}\left( {m{{N}_{{ALS}}}\left( {kT\sum\limits_{i = 1}^d \,{{n}_{i}} + {{T}^{2}}{\kern 1pt} \sum\limits_{i = 1}^d \,n_{i}^{2}} \right)} \right).$

Из приведенных рассуждений легко видеть, что алгоритм имеет полиномиальную сложность относительно параметров $m,{{n}_{1}}, \ldots ,{{n}_{d}}$. Более того, сложность линейна по $m$.

4. ТЕОРЕТИЧЕСКИЙ АНАЛИЗ

В этом разделе мы не ставим целью доказать строгие результаты о сходимости метода тензоризованного OMP. Мы лишь предполагаем дать более глубокое понимание различных шагов алгоритма.

Тензоризованное решение задачи наименьших квадратов малого тензорного ранга представляет собой довольно сложный для анализа объект. Кроме того, согласно одному из базовых предположений, нормальное решение задачи наименьших квадратов приближается тензором малого ранга. Поэтому вместо анализа решений малого ранга начнем с анализа нормального решения общей задачи наименьших квадратов для случайных матриц. Для достаточно богатого класса случайных матриц, в частности, матриц из нормального распределения, можно показать наличие свойства ограниченной изометрии [9]. Как было показано выше, для матриц, обладающих свойством ограниченной изометрии, нормальное решение задачи наименьших квадратов

$v = \mathop {\arg {\text{min}}}\limits_u {{\left\| {\mathcal{A}u - y} \right\|}_{2}}$

является хорошей оценкой на величину

${{c}_{j}} = \frac{{({{a}_{j}},y)}}{{{{{\left\| {{{a}_{j}}} \right\|}}_{2}}}}.$

В качестве подтверждения этому был проведен численный эксперимент. Были сгенерированы случайная матрица $A$ из нормального распределения размера $m \times {{2}^{m}}$ при m = 20 и случайный вектор правой части $y$. Для каждого $j = 1, \ldots {{,2}^{m}}$ были вычислены значение ${{{v}}_{j}}$ компоненты решения задачи наименьших квадратов и величина $({{a}_{j}},y)$. На фиг. 1 изображены пары точек $(({{a}_{j}},y),{\text{|}}{{v}_{j}}{\kern 1pt} {\text{|}}{\kern 1pt} )$. Отсюда легко видеть, что компоненты $\left| {{{v}_{j}}} \right|$ являются хорошими оценками на величины ${\text{|}}{{c}_{j}}{\kern 1pt} {\text{|}}$ при малых значениях $\delta $.

Фиг. 1.

График зависимости модуля компонент решения ${\text{|}}{{v}_{j}}{\text{|}}$ от скалярного произведения $({{a}_{j}},y)$ для случайной матрицы из нормального распределения. Размер матрицы $20 \times {{2}^{{20}}}$.

4.1. Влияние ортогонализации

Опишем, что изменится, если вместо задачи наименьших квадратов для $Ax = y$ решается задача

${{\left\| {QAx - Qy} \right\|}_{2}} \to \min .$

Рассуждения этого пункта верны для любых матриц. Пусть имеются матрица $A \in {{\mathbb{R}}^{{m \times n}}}$ и вектор правой части $y$, и на первом шаге метода мы некоторым образом нашли вектор ${{a}^{{(1)}}}$ и ортогонализуем к нему матрицу. Матрица ортогонализации может быть представлена в виде

${{Q}_{1}} = I - {{a}^{{(1)}}}{{({{a}^{{(1)}}})}^{{\text{т}}}}{\text{/}}\left\| {{{a}^{{(1)}}}} \right\|_{2}^{2}.$

Таким образом, решается система ${{Q}_{1}}A = {{Q}_{1}}y$. Пусть далее мы нашли столбец ${{a}^{{(2)}}}$ и ортогонализуем систему к нему. Поскольку ${{a}^{{(2)}}}$ выбирался из столбцов матрицы ${{Q}_{1}}A$, векторы ${{a}^{{(1)}}}$ и ${{a}^{{(2)}}}$ будут ортогональны. Матрица ${{Q}_{2}}$ строится аналогичным образом:

${{Q}_{2}} = I - {{a}^{{(2)}}}{{({{a}^{{(2)}}})}^{{\text{т}}}}{\text{/}}\left\| {{{a}^{{(2)}}}} \right\|_{2}^{2}$

и решается система ${{Q}_{2}}{{Q}_{1}}A = {{Q}_{2}}{{Q}_{1}}y$. Верна следующая

Теорема 1. Пусть $A \in {{\mathbb{R}}^{{m \times n}}}$ и $y \in {{\mathbb{R}}^{m}}$. Пусть ${{Q}_{1}},{{Q}_{2}}, \ldots ,{{Q}_{k}}$ – набор матриц ортогонализации, порожденных попарно ортогональными векторами. В этом случае решения задач наименьших квадратов систем

(6)

${{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}A = {{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y$

(7)

$A = {{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y$

совпадают.

Доказательство. Докажем, что псевдообратная матрица к матрице ортогонализации ${{Q}_{i}}$ совпадает с ней самой, т.е. $Q_{i}^{\dag } = {{Q}_{i}}$. Легко проверить, что

$Q_{i}^{2} = {{Q}_{i}},\quad Q_{i}^{{\text{т}}} = {{Q}_{i}},$

откуда следуют все аксиомы псевдообратной матрицы. Далее покажем, что матрицы ${{Q}_{i}}$ и ${{Q}_{j}}$ коммутируют. Действительно,

$\begin{gathered} {{Q}_{i}}{{Q}_{j}} = \left( {I - {{a}^{{(i)}}}{{{\left( {{{a}^{{(i)}}}} \right)}}^{{\text{т}}}}{\kern 1pt} {\text{/}}\left\| {{{a}^{{(i)}}}} \right\|_{2}^{2}} \right)\left( {I - {{a}^{{(j)}}}{{{\left( {{{a}^{{(j)}}}} \right)}}^{{\text{т}}}}{\kern 1pt} {\text{/}}\left\| {{{a}^{{(j)}}}} \right\|_{2}^{2}} \right) = \\ = I - {{a}^{{(i)}}}{{\left( {{{a}^{{(i)}}}} \right)}^{{\text{т}}}}{\kern 1pt} {\text{/}}\left\| {{{a}^{{(i)}}}} \right\|_{2}^{2} - {{a}^{{(j)}}}{{\left( {{{a}^{{(j)}}}} \right)}^{{\text{т}}}}{\kern 1pt} {\text{/}}\left\| {{{a}^{{(j)}}}} \right\|_{2}^{2} + \frac{1}{{\left\| {{{a}^{{(i)}}}} \right\|_{2}^{2}\left\| {{{a}^{{(j)}}}} \right\|_{2}^{2}}}{{a}^{{(i)}}}{{\left( {{{a}^{{(i)}}}} \right)}^{{\text{т}}}}{{a}^{{(j)}}}{{\left( {{{a}^{{(j)}}}} \right)}^{{\text{т}}}}. \\ \end{gathered} $

Но последнее слагаемое равно $0$ в силу того, что ${{({{a}^{{(i)}}})}^{{\text{т}}}}{{a}^{{(j)}}} = 0$, откуда видно, что матрицы ${{Q}_{i}}$ и ${{Q}_{j}}$ коммутируют.

Решение системы (6) представляется в виде

(8)

$\begin{gathered} {{\left( {{{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}A} \right)}^{\dag }}{{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y = {{A}^{\dag }}Q_{1}^{\dag }Q_{2}^{\dag } \ldots Q_{n}^{\dag }{{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y = \\ = {{A}^{\dag }}{{Q}_{1}}{{Q}_{2}} \ldots {{Q}_{n}}{{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y = {{A}^{\dag }}Q_{n}^{2}Q_{{n - 1}}^{2} \ldots Q_{1}^{2}y = {{A}^{\dag }}{{Q}_{n}}{{Q}_{{n - 1}}} \ldots {{Q}_{1}}y, \\ \end{gathered} $

что в точности совпадает с решением задачи наименьших квадратов (7). Теорема доказана.

Из этого видно, что решения задач

${{\left\| {QAx - Qy} \right\|}_{2}} \to \min $

${{\left\| {Ax - Qy} \right\|}_{2}} \to \min $

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

Заметим, что это утверждение верно лишь для нормального решения задачи наименьших квадратов и может оказаться неверным для решения задачи наименьших квадратов малого ранга.

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

${{c}_{j}} = \frac{{({{a}_{j}},r)}}{{{{{\left\| {{{a}_{j}}} \right\|}}_{2}}}}.$

Однако на практике мы не можем найти точное решение задачи наименьших квадратов для больших систем, поэтому алгоритм тензоризованного OMP выбирает $C$ кандидатов на основе компонент малорангового решения. С учетом предположения, что точное решение задачи наименьших квадратов может быть хорошо приближено тензором малого ранга, это дает возможность надеяться, что в множество из $C$ кандидатов попадет столбец, который бы выбрал алгоритм OMP.

4.2. О падении невязки на различных шагах метода

Изучим вопрос о скорости падения невязки для матрицы $A$, имеющей вид

(9)

$A = \left[ {\begin{array}{*{20}{c}} 1&1&1&1& \ldots &{ - 1}&{ - 1} \\ 1&1&1&1& \ldots &{ - 1}&{ - 1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 1&1&1&1& \ldots &{ - 1}&{ - 1} \\ 1&1&{ - 1}&{ - 1}& \ldots &{ - 1}&{ - 1} \\ 1&{ - 1}&1&{ - 1}& \ldots &1&{ - 1} \end{array}} \right].$

Эта матрица содержит всевозможные столбцы из $ \pm 1$, взятые в лексикографическом порядке. Легко показать, что эта матрица обладает свойством ограниченной изометрии с $\delta = 0$. Пусть правая часть задана в виде суммы $k$ столбцов матрицы $A$

$b = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}{{a}_{{{{i}_{j}}}}} + \varepsilon .$

Пусть произошла ортогонализация к вектору $p$, где ${{\left\| p \right\|}_{2}} = 1$. Разложим столбцы матрицы $A$, входящие в правую часть, на компоненты, коллинеарные $p$ и ортогональные ему:

${{a}_{{{{i}_{j}}}}} = a_{{{{i}_{j}}}}^{\parallel } + a_{{{{i}_{j}}}}^{ \bot },$

откуда $a_{{{{i}_{j}}}}^{\parallel } = {{\beta }_{j}}p$. Тогда

$b = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{\parallel } + {{\varepsilon }^{\parallel }} + \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{ \bot } + {{\varepsilon }^{ \bot }}.$

Тогда после ортогонализации

${{b}_{1}} = \left( {I - p{{p}^{{\text{т}}}}} \right)b = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{\parallel } + {{\varepsilon }^{\parallel }} + \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{ \bot } + {{\varepsilon }^{ \bot }} - \left( {\sum\limits_{j = 1}^k \,{{\alpha }_{j}}(p,a_{{{{i}_{j}}}}^{\parallel }) + (p,{{\varepsilon }^{\parallel }})} \right)p.$

Откуда, используя $a_{{{{i}_{j}}}}^{\parallel } = {{\beta }_{j}}p$ и ${{\varepsilon }^{\parallel }} = \gamma p$, получаем

${{b}_{1}} = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{ \bot } + {{\varepsilon }^{ \bot }}.$

Добавим и вычтем коллинеарную составляющую

${{b}_{1}} = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{ \bot } + {{\varepsilon }^{ \bot }} + \sum\limits_{j = 1}^k \,{{\alpha }_{j}}a_{{{{i}_{j}}}}^{\parallel } + {{\varepsilon }^{\parallel }} - \sum\limits_{j = 1}^k \,{{\alpha }_{j}}{{\beta }_{j}}p - \gamma p = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}{{a}_{{{{i}_{j}}}}} + \tau p + \varepsilon {\kern 1pt} .$

Заметим, что $p$ с точностью до нормировки является одним из столбцов матрицы $A$. Таким образом, после шага метода в правой части становится на 1 столбец больше.

Оценим норму вектора правой части ${{b}_{1}}$ после ортогонализации. По теореме Пифагора имеем

$\left\| b \right\|_{2}^{2} = \left\| {{{b}_{1}}} \right\|_{2}^{2} + \frac{{{{{(p,b)}}^{2}}}}{{\left\| p \right\|_{2}^{2}}},$

откуда

$\frac{{\left\| {{{b}_{1}}} \right\|_{2}^{2}}}{{\left\| b \right\|_{2}^{2}}} = 1 - \frac{{{{{(p,b)}}^{2}}}}{{\left\| p \right\|_{2}^{2}\left\| b \right\|_{2}^{2}}},$

что равносильно

$\frac{{{{{\left\| {{{b}_{1}}} \right\|}}_{2}}}}{{{{{\left\| b \right\|}}_{2}}}} = \left| {\sin \angle (p,b)} \right|.$

Оценим для начала худший случай. Выберем нормированный вектор $b$ так, что $\mathop {\max }\limits_p {{(p,b)}^{2}}$ минимально. Поскольку $p$ состоит из $ \pm 1$, можно, не ограничивая общности, считать, что все компоненты вектора $b$ неотрицательны и вектор $p$ состоит из всех $1$. Тогда нужно найти вектор $b$ такой, что ${{\left\| b \right\|}_{2}} = 1$, все элементы неотрицательны и их сумма минимальна. Очевидно, это вектор $b = {{e}_{1}}$, имеющий единицу в 1 позиции и нули в остальных.

Тогда в худшем случае имеем

$\frac{{\left\| {{{b}_{1}}} \right\|_{2}^{2}}}{{\left\| b \right\|_{2}^{2}}} = 1 - \frac{1}{m}.$

Заметим, что этот случай вовсе не является нереалистичным:

$b = \frac{1}{2}\left[ {\begin{array}{*{20}{c}} 1 \\ 1 \\ \vdots \\ 1 \end{array}} \right] - \frac{1}{2}\left[ {\begin{array}{*{20}{c}} { - 1} \\ 1 \\ \vdots \\ 1 \end{array}} \right].$

Однако в среднем все не так плохо. Были проведены следующие численные эксперименты. Для правой части выбиралось $k$ случайных столбцов матрицы и $k$ коэффициентов из нормального распределения. Находилось падение квадрата нормы правой части и результат усреднялся по ${{10}^{6}}$ экспериментов. На фиг. 2 приведен график зависимости $1 - \left\| {{{b}_{1}}} \right\|_{2}^{2}{\text{/}}\left\| b \right\|_{2}^{2}$ от количества столбцов $k$. Легко видеть, что среднее падение всегда больше $0.5$ и очень слабо зависит от числа строк матрицы $m$.

Фиг. 2.

График зависимости квадрата падения нормы невязки от количества компонент в правой части.

Оценим величину

(10)

$\mathop {\max }\limits_p \frac{{{{{(p,b)}}^{2}}}}{{{{{\left\| p \right\|}}^{2}}{{{\left\| b \right\|}}^{2}}}}.$

Во-первых, ясно, что $\left\| p \right\|_{2}^{2} = m$. Пусть

$b = \sum\limits_{j = 1}^k \,{{\alpha }_{j}}{{a}_{{{{i}_{j}}}}}.$

Тогда все компоненты имеют вид

${{b}_{i}} = \pm {{\alpha }_{1}} \pm {{\alpha }_{2}} \pm \ldots \pm {{\alpha }_{k}},$

где ${{\alpha }_{j}} \sim \mathcal{N}(0,1)$, и знаки + и – выбираются с равной вероятностью. Легко видеть, что $ \pm {{\alpha }_{j}} \sim \mathcal{N}(0,1)$, следовательно,

${{b}_{i}} \sim \mathcal{N}(0,k){\kern 1pt} .$

Однако компоненты вектора $b$ не являются независимыми.

Далее, из (10) легко посчитать максимум

$\mathop {\max }\limits_p \frac{{{{{(p,b)}}^{2}}}}{{{{{\left\| p \right\|}}^{2}}{{{\left\| b \right\|}}^{2}}}} = \frac{{\left\| b \right\|_{1}^{2}}}{{m\left\| b \right\|_{2}^{2}}}.$

Нетрудно понять, что, чем больше слагаемых в $b$ (т.е. чем больше $k$), тем более независимы между собой будут компоненты ${{b}_{i}}$, и наоборот – чем меньше $k$, тем более зависимы компоненты. Так, если $k = 1$, то будет всего один вариант ${\text{|}}{{\alpha }_{1}}{\text{|}}$, если $k = 2$, то $2$ варианта: ${\text{|}}{{\alpha }_{1}} + {{\alpha }_{2}}{\kern 1pt} {\text{|}}$ и ${\text{|}}{{\alpha }_{1}} - {{\alpha }_{2}}{\kern 1pt} {\text{|}}$. В общем случае ${{2}^{{k - 1}}}$ вариантов. Величина $\frac{{\left\| b \right\|_{1}^{2}}}{{\left\| b \right\|_{2}^{2}}}$ будет максимальна, если все компоненты равны, т.е. наиболее зависимы. Из этого наблюдения и из фиг. 2 можно видеть, что с ростом $k$ (по оси $x$) величина (10) убывает. Отсюда вытекает гипотеза, что оценкой снизу на (10) будет ситуация, когда все компоненты независимы. Доказательством этой гипотезы мы не обладаем, однако, численные эксперименты хорошо согласуются с ней. В случае если все компоненты ${{b}_{i}}$ независимые и одинаково распределенные с ${{b}_{i}} \sim \mathcal{N}(0,1)$ (в случае ${{b}_{i}} \sim \mathcal{N}(0,k)$ можно отнормировать числитель и знаменатель), имеем

$\mathbb{E}\left( {\frac{{\left\| b \right\|_{1}^{2}}}{{\left\| b \right\|_{2}^{2}}}} \right) = 1 + \frac{2}{\pi }(m - 1).$

Для сравнения эти оценки изображены на фиг. 2 штриховыми линиями. Эту оценку можно доказать аналитически.

Действительно,

$\left\| b \right\|_{1}^{2} = ({\kern 1pt} {\text{|}}{{b}_{1}}{\kern 1pt} {\text{|}}\; + \;{\text{|}}{{b}_{2}}{\kern 1pt} {\text{|}}\; + \ldots + \;{\text{|}}{{b}_{m}}{\kern 1pt} {\text{|}}{\kern 1pt} {{)}^{2}} = \sum\limits_{j = 1}^m \,b_{j}^{2} + \sum\limits_{i,j = 1,i \ne j}^k {\kern 1pt} {\text{|}}{{b}_{i}}{{b}_{j}}{\kern 1pt} {\text{|}}{\kern 1pt} {\kern 1pt} {\text{.}}$

Тогда

$\frac{{\left\| b \right\|_{1}^{2}}}{{\left\| b \right\|_{2}^{2}}} = 1 + \sum\limits_{i,j = 1,i \ne j}^k \frac{{{\text{|}}{{b}_{i}}{{b}_{j}}{\text{|}}{\kern 1pt} }}{{b_{1}^{2} + \ldots + b_{m}^{2}}},$

откуда имеем, что

$\mathbb{E}\left( {\frac{{\left\| b \right\|_{1}^{2}}}{{\left\| b \right\|_{2}^{2}}}} \right) = 1 + m(m - 1)\mathbb{E}\left( {\frac{{{\text{|}}{{b}_{1}}{{b}_{2}}{\kern 1pt} {\text{|}}}}{{b_{1}^{2} + \ldots + b_{m}^{2}}}} \right).$

Последнее матожидание равно

$\begin{gathered} \mathbb{E}\left( {\frac{{{\text{|}}{{b}_{1}}{{b}_{2}}{\kern 1pt} {\text{|}}}}{{b_{1}^{2} + \ldots + b_{m}^{2}}}} \right) = \frac{1}{{{{{(2\pi )}}^{{m/2}}}}}\int\limits_{ - \infty }^\infty \ldots \int\limits_{ - \infty }^\infty \frac{{{\text{|}}{{x}_{1}}{{x}_{2}}{\kern 1pt} {\text{|}}}}{{x_{1}^{2} + \ldots + x_{m}^{2}}}\exp \left( { - \frac{{x_{1}^{2} + \ldots + x_{m}^{2}}}{2}} \right)d{{x}_{1}} \ldots d{{x}_{m}} = \\ = \frac{{{{2}^{{m/2}}}}}{{{{\pi }^{{m/2}}}}}\int\limits_0^\infty \ldots \int\limits_0^\infty \frac{{{{x}_{1}}{{x}_{2}}}}{{x_{1}^{2} + \ldots + x_{m}^{2}}}\exp \left( { - \frac{{x_{1}^{2} + \ldots + x_{m}^{2}}}{2}} \right)d{{x}_{1}} \ldots d{{x}_{m}}. \\ \end{gathered} $

Вычислим при $\varepsilon > 0$ интеграл

$\int\limits_0^\infty \ldots \int\limits_0^\infty {\kern 1pt} \int\limits_\varepsilon ^\infty {\kern 1pt} \int\limits_\varepsilon ^\infty \frac{{{{x}_{1}}{{x}_{2}}}}{{x_{1}^{2} + \ldots + x_{m}^{2}}}\exp \left( { - \frac{{x_{1}^{2} + \ldots + x_{m}^{2}}}{2}} \right)d{{x}_{1}} \ldots d{{x}_{m}}.$

Сделаем замену переменных

$\begin{array}{*{20}{c}} {{{y}_{1}} = x_{1}^{2},} \\ {{{y}_{2}} = x_{2}^{2},} \\ {{{y}_{3}} = {{x}_{3}},} \\ \vdots \\ {{{y}_{m}} = {{x}_{m}}.} \end{array}$

Легко видеть, что якобиан такого преобразования равен $\frac{1}{{4{{x}_{1}}{{x}_{2}}}}$, откуда последний интеграл равен

$\frac{1}{4}\int\limits_0^\infty \ldots \int\limits_0^\infty {\kern 1pt} \int\limits_{{{\varepsilon }^{2}}}^\infty {\kern 1pt} \int\limits_{{{\varepsilon }^{2}}}^\infty \frac{{\exp \left( { - \frac{{{{y}_{1}} + {{y}_{2}} + y_{3}^{2} + \ldots + y_{m}^{2}}}{2}} \right)}}{{{{y}_{1}} + {{y}_{2}} + y_{3}^{2} + \ldots + y_{m}^{2}}}d{{y}_{1}} \ldots d{{y}_{m}}.$

Сделаем еще одну замену

$\begin{array}{*{20}{c}} {{{t}_{1}} = \frac{{{{y}_{1}} + {{y}_{2}} + y_{3}^{2} + \ldots + y_{m}^{2}}}{2},} \\ {{{t}_{2}} = {{y}_{2}},} \\ {{{t}_{3}} = {{y}_{3}},} \\ \vdots \\ {{{t}_{m}} = {{y}_{m}}.} \end{array}$

В этом случае якобиан равен $2$ и интеграл принимает вид

(11)

$\frac{1}{4}\int\limits_0^\infty \ldots \int\limits_0^\infty {\kern 1pt} \int\limits_{{{\varepsilon }^{2}}}^\infty \int\limits_{\frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}}^\infty \frac{{{{e}^{{ - {{t}_{1}}}}}}}{{{{t}_{1}}}}d{{t}_{1}} \ldots d{{t}_{m}}.$

Введем обозначение для экспоненциального интеграла

${{E}_{n}}(x) = \int\limits_1^\infty \frac{{{{e}^{{ - xt}}}}}{{{{t}^{n}}}}dt.$

Заметим, что

${{E}_{1}}(x) = \int\limits_1^\infty \frac{{{{e}^{{ - xt}}}}}{t}dt$

и сделаем замену $y = xt$ при $x > 0$. Тогда получим

${{E}_{1}}(x) = \int\limits_x^\infty \frac{{{{e}^{{ - y}}}}}{y}dy.$

В этих обозначениях интеграл (11) равен

$\frac{1}{4}\int\limits_0^\infty \ldots \int\limits_0^\infty {\kern 1pt} \int\limits_{{{\varepsilon }^{2}}}^\infty {{E}_{1}}\left( {\frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right)d{{t}_{2}} \ldots d{{t}_{m}}.$

Вычислим производную

$\left( {{{E}_{n}}\left( {\frac{x}{2} + c} \right)} \right)_{x}^{'} = \int\limits_1^\infty \frac{{{{e}^{{ - \left( {\frac{x}{2} + c} \right)t}}}\left( { - \frac{t}{2}} \right)}}{{{{t}^{n}}}}dt = - \frac{1}{2}{{E}_{{n - 1}}}\left( {\frac{x}{2} + c} \right).$

${{E}_{n}}\left( {\frac{x}{2} + c} \right) = - 2\left( {{{E}_{{n + 1}}}\left( {\frac{x}{2} + c} \right)} \right)_{x}^{'}.$

Посчитаем неопределенный интеграл

(12)

$\begin{gathered} \int \,{{E}_{1}}\left( {\frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right)d{{t}_{2}} = - 2\int \left( {{{E}_{2}}\left( {\frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right)} \right)_{{{{t}_{2}}}}^{'}d{{t}_{2}} = \\ = - 2{\kern 1pt} \int\limits_1^\infty \frac{{{{e}^{{ - \frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}{{t}_{1}}}}}}}{{t_{1}^{2}}}d{{t}_{1}} + C. \\ \end{gathered} $

Откуда имеем, что

$\int\limits_{{{\varepsilon }^{2}}}^\infty \,{{E}_{1}}\left( {\frac{{{{\varepsilon }^{2}} + {{t}_{2}} + t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right)d{{t}_{2}} = 2\int\limits_1^\infty \frac{{{{e}^{{ - \left( {{{\varepsilon }^{2}} + \frac{{t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right){{t}_{1}}}}}}}{{t_{1}^{2}}}d{{t}_{1}}.$

Переходя к пределу при $\varepsilon \to 0$, получаем, что

$\begin{gathered} \int\limits_0^\infty \ldots \int\limits_0^\infty \frac{{{{x}_{1}}{{x}_{2}}}}{{x_{1}^{2} + \ldots + x_{m}^{2}}}\exp \left( { - \frac{{x_{1}^{2} + \ldots + x_{m}^{2}}}{2}} \right)d{{x}_{1}} \ldots d{{x}_{m}} = \\ \, = \frac{1}{2}\int\limits_0^\infty \ldots \int\limits_0^\infty {\kern 1pt} \int\limits_1^\infty \frac{{{{e}^{{ - \left( {\frac{{t_{3}^{2} + \ldots + t_{m}^{2}}}{2}} \right){{t}_{1}}}}}}}{{t_{1}^{2}}}d{{t}_{1}}d{{t}_{3}} \ldots d{{t}_{m}} = \frac{1}{2}\int\limits_1^\infty \frac{1}{{t_{1}^{2}}}\left( {\int\limits_0^\infty \,{{e}^{{ - \frac{{t_{3}^{2}{{t}_{1}}}}{2}}}}d{{t}_{3}}} \right) \ldots \left( {\int\limits_0^\infty \,{{e}^{{ - \frac{{t_{m}^{2}{{t}_{1}}}}{2}}}}d{{t}_{m}}} \right)d{{t}_{1}}. \\ \end{gathered} $

Простой заменой получаем, что

$\int\limits_0^\infty \,{{e}^{{ - \frac{{t_{3}^{2}{{t}_{1}}}}{2}}}}d{{t}_{3}} = \sqrt {\frac{\pi }{{2{{t}_{1}}}}} .$

Тогда имеем

$\int\limits_0^\infty \ldots \int\limits_0^\infty \frac{{{{x}_{1}}{{x}_{2}}}}{{x_{1}^{2} + \ldots + x_{m}^{2}}}\exp \left( { - \frac{{x_{1}^{2} + \ldots + x_{m}^{2}}}{2}} \right)d{{x}_{1}} \ldots d{{x}_{m}} = \frac{{{{\pi }^{{\frac{{m - 2}}{2}}}}}}{{{{2}^{{\frac{m}{2}}}}}}\int\limits_1^\infty \frac{{d{{t}_{1}}}}{{{{t}^{{\frac{{m + 2}}{2}}}}}} = \frac{{{{\pi }^{{\frac{{m - 2}}{2}}}}}}{{{{2}^{{\frac{{m - 2}}{2}}}}m}}.$

И в итоге получаем

$\mathbb{E}\frac{{{\text{|}}{{b}_{1}}{{b}_{2}}{\text{|}}}}{{b_{1}^{2} + \ldots + b_{m}^{2}}} = \frac{{{{2}^{{\frac{m}{2}}}}}}{{{{\pi }^{{\frac{m}{2}}}}}}\frac{{{{\pi }^{{\frac{{m - 2}}{2}}}}}}{{{{2}^{{\frac{{m - 2}}{2}}}}m}} = \frac{2}{{\pi m}},$

$\mathbb{E}\left( {\frac{{\left\| b \right\|_{1}^{2}}}{{\left\| b \right\|_{2}^{2}}}} \right) = 1 + m(m - 1)\frac{2}{{\pi m}} = 1 + \frac{2}{\pi }(m - 1),$

$\mathbb{E}\left( {\mathop {\max }\limits_p \frac{{{{{(p,b)}}^{2}}}}{{{{{\left\| b \right\|}}^{2}}{{{\left\| b \right\|}}^{2}}}}} \right) = \frac{2}{\pi } + \left( {1 - \frac{2}{\pi }} \right)\frac{1}{m}.$

Тогда матожидание падения нормы правой части оказывается равно

$\mathbb{E}\frac{{\left\| {{{b}_{1}}} \right\|_{2}^{2}}}{{\left\| b \right\|_{2}^{2}}} = \left( {1 - \frac{2}{\pi }} \right)\left( {1 - \frac{1}{m}} \right),$

где ${{b}_{1}}$ – правая часть после первого шага метода.

Таким образом, для матриц вида (9) на первом шаге невязка будет падать в среднем более чем в 2 раза. Это позволяет надеяться на высокую скорость сходимости метода. Однако такое быстрое падение невязки является особенностью того, что столбцы матрицы (9) очень плотно покрывают все пространство. Для случайных матриц, не обладающих такой особенностью, падение невязки ожидается более медленным.

5. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

В этом разделе приведены результаты численных экспериментов для предложенного метода.

5.1. Описание модели

Модель, рассматриваемая в этом разделе, носит название модель EMP-Вольтерра. Введем обозначения. Пусть $X,Y \in {{\mathbb{C}}^{N}}$ и обозначают входы и выходы модели соответственно. Здесь $N$ обозначает длину сигнала, $P = \left\lceil {\frac{{\max {{X}_{i}}}}{{{{2}^{K}}}}} \right\rceil $, где $K$ является нормировочным параметром модели. Матрица $F \in {{\mathbb{C}}^{{N \times (P + 1)}}}$ является некоторой функцией от входа $X$ и считается известной, $e$ обозначает вектор из единиц. Параметрами модели, подлежащими определению, являются ${{u}_{{rd}}} \in {{\mathbb{C}}^{{{{S}_{X}}}}}$, ${{v}_{{rd}}} \in {{\mathbb{C}}^{{{{S}_{X}}}}}$, ${{w}_{{rd}}} \in {{\mathbb{C}}^{{{{S}_{L}}}}}$, ${{a}_{{rd}}} \in {{\mathbb{C}}^{{P + 1}}}$, ${{\alpha }_{{rd}}},{{\beta }_{{rd}}} \in \mathbb{C}$. Тогда модель EMP-Вольтерра может быть записана следующим образом:

$Y \approx \sum\limits_{r = 1}^R {\left( {{\kern 1pt} \mathop \odot \limits_{d = 1}^{\left\lceil {\frac{{{{D}_{X}}}}{2}} \right\rceil } \left( {X \circledast {{u}_{{rd}}} + {{\alpha }_{{rd}}}e} \right)} \right)} \odot \left( {{\kern 1pt} \mathop \odot \limits_{d = 1}^{\left\lfloor {\frac{{{{D}_{X}}}}{2}} \right\rfloor } \left( {\bar {X} \circledast {{{v}}_{{rd}}} + {{\beta }_{{rd}}}e} \right)} \right) \odot \left( {{\kern 1pt} \mathop \odot \limits_{d = 1}^{{{D}_{L}}} \left( {\left( {F{{\alpha }_{{rd}}}} \right) \circledast {{w}_{{rd}}}} \right)} \right).$

Здесь символом $ \odot $ обозначено адамарово произведение, а символом $ \circledast $ свертка между двумя векторами.

Несложно видеть, что эту модель можно записать как тензоризованный линейный оператор относительно $[{{u}_{{rd}}},{{\alpha }_{{rd}}}]$, $[{{v}_{{rd}}},{{\beta }_{{rd}}}]$, ${{w}_{{rd}}}$ и ${{a}_{{rd}}}$. Число столбцов матрицы такого оператора может быть вычислено по формуле

${{N}_{{COL}}} = ({{S}_{X}} + {{1)}^{{{{D}_{X}}}}} \times S_{L}^{{{{D}_{L}}}} \times {{(P + 1)}^{{{{D}_{L}}}}}.$

В качестве входного сигнала $X$ и выходного сигнала $Y$ были использованы реальные данные для решения задачи цифрового предыскажения (Digital Predistortion, DPD). Длина сигнала $N \sim 5 \times {{10}^{4}}$. Наименьшие размеры модели, которая обеспечивает достаточное для практического применения качество цифрового предыскажения для данного сигнала, равны ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 1$. В этом случае число столбцов получающейся матрицы равно 39,424. Чуть более сложная модель получается выбором ${{S}_{X}} = 7,\;{{S}_{L}} = 7,\;{{D}_{X}} = 3,\;{{D}_{L}} = 3$, что соответствует 233,744,896 столбцам. Максимальный размер задачи, для которой проводились эксперименты, соответствует параметрам ${{S}_{X}} = 13,$ ${{S}_{L}} = 13,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 3$, где матрица содержит приблизительно $8 \times {{10}^{9}}$ столбцов и $5 \times {{10}^{4}}$ строчек.

В качестве метрики качества работы была взята классическая для теории обработки сигналов единица измерения – децибелы. Таким образом, точность вычисляется по формуле

$ - 20{{\lg }_{{10}}}\frac{{{{{\left\| {\mathcal{A}u - Y} \right\|}}_{2}}}}{{{{{\left\| Y \right\|}}_{2}}}},$

где $u$ – найденное разреженное решение.

5.2. Результаты

Ниже приведены результаты применения метода тензоризованного OMP к описанной выше модели. На фиг. 3 приведен график сравнения алгоритма OMP и метода тензоризованного OMP для модели с параметрами ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 1$. Здесь и далее алгоритм тензоризованного OMP отбирает $C = 400$ столбцов-кандидатов для множества $\hat {S}$. В качестве метода поиска малорангового решения вместо алгоритма ALS использовался метод Левенберга-Марквардта [10], [11]. Как отмечалось выше, на каждой итерации метода тензоризованного OMP в пункте 4 Алгоритма 2 можно выбирать не 1 столбец, а несколько. Поскольку нахождение малорангового решения занимает большую часть времени работы, такой прием позволяет существенно ускорить время вычислений. Для сравнения на фиг. 3 приведен также график для метода тензоризованного OMP с выбором 5 столбцов на каждой итерации. Из графиков, во-первых, можно видеть, что выбор 5 столбцов вместо одного практически не ухудшает точности построенного решения, но позволяет почти в 5 раз ускорить работу программы. Поэтому во всех дальнейших экспериментах в методе тензоризованного OMP выбирает 5 столбцов за итерацию. Во-вторых, из графиков видно, что метод тензоризованного OMP работает ненамного хуже классического OMP, хотя имеет полиномиальную сложность, в отличие от OMP, имеющего экспоненциальную сложность. На данном этапе качество и не могло получиться лучше, чем у OMP, поскольку метод тензоризованного OMP в некотором смысле является его аппроксимацией, предназначенной для работы со сверхбольшими тензоризованными матрицами.

Фиг. 3.

График сравнения методов OMP и тензоризованного OMP с различным числом шагов в пункте 4 Алгоритма 2. Модель ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 1$. Для метода тензоризованного OMP отбиралось $C = $ = 400 столбцов.

Метод тензоризованного OMP может быть дополнительно улучшен с помощью следующего подхода. На шагах 3–5 Алгоритма 2 будем строить матрицу $\hat {A}$ содержащей не только столбцы из $\hat {S}$, но также и все столбцы из множества уже отобранных ранее столбцов $S$. Пусть на некоторой итерации Алгоритма 2 уже выбрано множество $S$. Получим множество новых столбцов-кандидатов $\hat {S}$. Построим матрицу $\hat {A}$, содержащую столбцы и $S \cup \hat {S}$, и отберем из нее (после ортогонализации) ${\text{|}}S{\kern 1pt} {\text{|}} + 5$ столбцов с помощью классического алгоритма OMP. Такой подход гарантированно не ухудшит качества работы, но может помочь найти хорошее решение с меньшим числом разреженных компонент, что можно видеть из фиг. 4.

Фиг. 4.

График сравнения методов простого тензоризованного OMP и тензоризованного OMP с переранжированием на каждом шаге. Модель ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 1$. Для метода тензоризованного OMP отбиралось $C = 400$ столбцов.

Заметим, что матрица ${{Q}^{k}}\hat {A}$, для которой применяется OMP в пункте 4 Алгоритма 2 имеет всего $C = 400$ столбцов, что позволяет применять здесь не только жадный алгоритм OMP, но и более продвинутые стратегии оптимизации. В качестве примера были использованы метод имитации отжига и алгоритм K-SVD [12], [13]. Метод имитации отжига представляет собой метод глобальной оптимизации, который на каждом шаге переходит либо в состояние, уменьшающее значение функционала, либо, с некоторой вероятностью, в состояние, ухудшающее значение функционала, причем вероятность последнего зависит от температурного коэффициента, значение которого убывает в процессе работы алгоритма [14]. Алгоритм K-SVD отличается от OMP своей возможностью не только жадно добавлять столбцы, но и выбрасывать ненужные. На фиг. 5 приведены результаты сравнения методов. Алгоритм K-SVD по построению всегда пересматривает множество ранее выбранных столбцов, поэтому для OMP и метода имитации отжига было использовано переранжирование всех столбцов на каждой итерации. Из графиков можно видеть, что на первых итерациях метод имитации отжига выбирает столбцы лучше, поскольку является менее жадным методом, однако затем лидерство переходит к K-SVD, поскольку задача оптимизации становится слишком сложной для метода имитации отжига.

Фиг. 5.

График сравнения различных методов отбора столбцов в пункте 4 Алгоритма 2. Использованные методы: OMP (с переранжированием после каждого шага), метод имитации отжига (с переранжированием после каждого шага) и K-SVD [12], [13]. Модель ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3,$ ${{D}_{L}} = 1$. Для метода тензоризованного OMP отбиралось $C = 400$ столбцов.

Одной из проблем метода тензоризованного OMP является борьба с влиянием уже выбранных столбцов. Если столбец $j$ уже был выбран ранее, то после ортогонализации он будет нулевым, поэтому коэффициент на $j$-й позиции может быть выбран любым, и это не скажется на невязке. В случае тензоризованного OMP ситуация несколько менее критична, поскольку малоранговость решения представляет собой своего рода регуляризацию. Однако возможность ставить произвольные значения в некоторые позиции может позволить находить лучшее с точки зрения невязки решение, что в корне противоречит тому, как дальше будут использованы компоненты решения. Большая по модулю компонента должна соответствовать более полезному для решения столбцу. Для того чтобы добиться этого, можно добавить простейший ${{l}_{2}}$–регуляризатор с коэффициентом $\lambda $. Минимизируемая функция в этом случае записывается как

(13)

$\left\| {Q\mathcal{A}\left( {\sum\limits_{t = 1}^{\text{т}} \,{{u}_{{1t}}} \otimes {{u}_{{2t}}} \otimes \ldots \otimes {{u}_{{dt}}}} \right) - Qb{\kern 1pt} } \right\|_{2}^{2} + \lambda \sum\limits_{k,t} \left\| {{{u}_{{kt}}}} \right\|_{2}^{2} \to \min .$

На фиг. 6 приведен график сравнения метода с нерегуляризованной моделью при ${{D}_{L}} = 1$ и ${{D}_{L}} = 3$ и регуляризованной модели с ${{D}_{L}} = 3$. Из графика можно видеть, что увеличение размера модели закономерно приводит к улучшения качества. Из уравнения (13) несложно заметить, что коэффициент регуляризации $\lambda $ необходимо выбирать адаптивно. Действительно, на первых итерациях метода невязка будет большой и необходима большая $\lambda $ для получения эффекта от регуляризации. В то же время слишком большая $\lambda $ будет приводить к нахождению нулевого решения в качестве оптимального. Поэтому $\lambda $ выбирается адаптивно: изначально выбирается достаточно большое значение коэффициента регуляризации, но когда поиск малорангового решения начинает давать большую невязку, коэффициент $\lambda $ уменьшается в несколько раз. Результат такого уменьшения можно видеть на графике в виде скачка.

Фиг. 6.

График сравнения качества работы на различных моделях, а также эффект от применения регуляризации. Во всех моделях ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3$. Для метода тензоризованного OMP отбиралось $C = 400$ столбцов.

Классический метод OMP может быть применен к модели с ${{D}_{L}} = 1$, но его использование для модели с ${{D}_{L}} = 3$ уже затруднительно. На фиг. 7 приведен график сравнения OMP для модели с ${{D}_{L}} = 1$ и тензоризованного OMP с ${{D}_{L}} = 3$. Можно видеть, что при большом числе компонент тензоризованный OMP начинает обходить классический OMP, что доказывает его эффективность для поиска разреженных решений сверхбольших систем линейных уравнений.

Фиг. 7.

График сравнения методов OMP и тензоризованного OMP. Во всех моделях ${{S}_{X}} = 7,$ ${{S}_{L}} = 7,$ ${{D}_{X}} = 3$. Для метода тензоризованного OMP отбиралось $C = 400$ столбцов.

6. ЗАКЛЮЧЕНИЕ

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

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

Mallat S.G., Zhang Z. Matching pursuits with time-frequency dictionaries // IEEE Transactions on Signal Proc. 1993. V. 41. № 12. P. 3397–3415.
Elad M. Sparse and redundant representations: from theory to applications in signal and image processing. Springer Science & Business Media. 2010.
Cai T.T., Wang L. Orthogonal matching pursuit for sparse signal recovery with noise // IEEE Transactions on Information Theory. 2011. V. 57. № 7. P. 4680–4688.
Candes E.J., Tao T. Decoding by linear programming // IEEE Transactions on Information Theory. 2005. V. 51. № 12. P. 4203–4215.
Davenport M.A., Wakin M.B. Analysis of orthogonal matching pursuit using the restricted isometry property // IEEE Transactions on Information Theory. 2010. V. 56. № 9. P. 4395–4401.
Tropp J.A. Greed is good: Algorithmic results for sparse approximation // IEEE Transactions on Information theory. 2004. V. 50. № 10. P. 2231–2242.
Lebedeva O. Tensor conjugate-gradient-type method for rayleigh quotient minimizationin block qtt-format // Russian Journal of Numerical Analysis and Math. Modelling. 2011. V. 26. № 5. P. 465–489.
Zheltkov D.A., Osinsky A. Global optimization algorithms using tensor trains. Internat. Conference on Large-Scale Scientific Computing. Springer. 2019. P. 197–202.
Vershynin R. Introduction to the non-asymptotic analysis of random matrices //arXiv preprint arXiv:1011.3027. 2010.
Levenberg K.A method for the solution of certain non-linear problems in least squares // Quarterly of Applied Math. 1944. V. 2. № 2. P. 164–168.
Marquardt D.W. An algorithm for least-squares estimation of nonlinear parameters // J. of the Society for In-dustrial and Applied Math. 1963. V. 11. № 2. P. 431–441.
Aharon M., Elad M., Bruckstein A. K-svd: An algorithm for designing overcomplete dictionaries for sparse representation // IEEE Transactions on Signal Proc. 2006. V. 54. № 11. P. 4311–4322.
Rubinstein R., Bruckstein A.M., Elad M. Dictionaries for sparse representation modeling // Proc. of the IEEE. 2010. V. 98. № 6. P. 1045–1057.
Kirkpatrick S., Gelatt Jr, C., Vecchi M. Optimization by Simulated Annealing // Science. 1983. V. 220. № 4598. P. 671–680.

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

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал