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

1. ВВЕДЕНИЕ

Пусть $X \in {{\mathbb{R}}^{{m \times n}}}$ – неизвестная матрица ранга $k,$ причем $k \ll min(m,n)$, а $\Omega $ – случайный набор пар индексов $(i,j)$, где $i \in 1,2, \ldots ,m$ и $j \in 1,2, \ldots ,n.$ Задача восстановления всей матрицы $X$ лишь по элементам ${{X}_{{ij}}}$ для $(i,j) \in \Omega $ называется задачей матричного восполнения (Matrix completion problem).

Современные приложения задачи матричного восполнения включают рекомендательные системы (см. [1]), обработку коррелирующих сигналов (см. [2]–[4]), машинное обучение (см. [5], [6]) и многое другое.

Одним из эффективных методов, применяемых для решения задачи восполнения, является алгоритм проекции на главные сингулярные пространства (Singular Value Projection method), предложенный в [7], [8] (далее для краткости мы будем называть этот метод SVP).

По своей структуре алгоритм SVP относится к классу методов проективного градиента. Каждая его итерация состоит из двух шагов: шага градиентного спуска и последующей проекции полученного приближения на многообразие матриц ранга $k.$ Cложность одной итерации SVP практически полностью определяется сложностью построения наилучшего приближения с помощью сингулярного разложения (SVD), применяемого к m × n-матрицам общего вида.

В работе рассматриваются способы ускорения метода SVP за счет замены медленного алгоритма SVD наилучшего приближения матрицами предписанного ранга $k$ на известные быстрые методы построения малоранговых приближений. При этом не предполагается высокой точности получаемых приближений. Напротив, достаточно такой точности, при которой алгоритм восполнения сохраняет свойство глобальной сходимости к искомой матрице $X,$ а возможное увеличение общего числа итераций компенсируется скоростью построения приближениями малого ранга на каждой итерации. Как будет показано, применение такого подхода во многих случаях приводит к существенному сокращению времени работы всего алгоритма. Семейство алгоритмов, использующих алгоритмы приближенного проектирования будем называть методом ASVP (Approximate Singular Value Projection).

В работе рассматриваются два быстрых алгоритма малоранговых приближений:

• метод проектирования на случайные пространства (см. [9], [10]);

• метод псевдо-крестовой аппроксимации, основанной на принципе большого проективного объема (см. [11]–[13]).

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

Как было отмечено в [4], алгоритм SVP может проявлять нерегулярное поведение. По ходу работы алгоритма гипотезы, используемые в доказательстве его сходимости могут начать нарушаться, даже если они верны для искомой матрицы и исходного приближения; в таких случаях на практике алгоритм расходится. Для восстановления сходимости метода в [4] предложено ограничивать длину шага в градиентном спуске. Последнее замедляет сходимость и увеличивает сложность SVP. В случае метода ASVP нерегулярное поведение имеет даже большее значение. Более того, на практике нерегулярное поведение является типичным! В работе дается анализ причин нерегулярного поведения алгоритмов SVP и ASVP и предлагаются изменения, которые сохраняют сходимость алгоритма без уменьшения скорости сходимости.

2. TEОРИЯ СХОДИМОСТИ АЛГОРИТМА SVP

Следуя [7], мы рассматриваем теорию алгоритма SVP в несколько более общей постановке задачи восполнения.

А именно, пусть на множестве $m \times n$-матриц задан аффинный оператор $\mathcal{A}:{{\mathbb{R}}^{{m \times n}}} \to {{\mathbb{R}}^{M}},$ где $M$ обозначает число аффинных соотношений на элементах $m \times n$-матриц. Тогда для произвольного вектора ${\mathbf{B}} \in {{\mathbb{R}}^{M}}$ задачу малорангового восполнения можно поставить как задачу оптимизации

с квадратичным функционалом ${{\psi }_{\mathcal{A}}}(X)$. Следуя [7], заметим, что градиент $\nabla {{\psi }_{\mathcal{A}}}(Y)$ функционала (1) в произвольной точке $Y \in {{\mathbb{R}}^{{m \times n}}}$ может быть записан в виде где оператор является транспонированным к $\mathcal{A}$ (везде предполагается, что в пространствах ${{\mathbb{R}}^{{m \times n}}}$ и ${{\mathbb{R}}^{M}}$ задано стандартное скалярное произведение).

Алгоритм SVP состоит из двух шагов: шага градиентного спуска и проекции нового приближения на многообразие матриц ранга не выше $k$. Другими словами, если ${{X}_{t}}$ – приближение, полученное на шаге $t$ алгоритма, то следующее приближение ${{X}_{{t + 1}}}$ получается в виде

где $\tau \in \mathbb{R}$ – некоторый шаг градиентного спуска, который мы определим позже, а – наилучший во фробениусовой норме проектор на множество матриц ранга не выше $k$. Для сходимости SVP алгоритма в [7] доказана следующая теорема.

Теорема 1 (см. [7]). Пусть ${{X}_{*}}$ – решение (1), для которого ${{\psi }_{\mathcal{A}}}({{X}_{*}}) = 0$, а оператор $\mathcal{A}$ удовлетворяет условию ограниченной изометрии вида

c параметром $0 < \sigma < 1$ для всех матриц $X$ ранга не выше $2k.$ Если дополнительното алгоритм SVP с постоянным шагом $\tau = 1{\text{/}}(1 + \sigma )$ сходится к решению ${{X}_{*}}$, а скорость сходимости определяется соотношением

Из (2) следует, что алгоритмическая сложность SVP определяется сложностью вычисления наилучшей проекции ${{\mathcal{P}}_{k}}$ на многообразие матриц ранга не выше $k.$ Если $m \approx n,$ то сложность такого вычисления, основанного на сингулярном разложении, будет порядка $\mathcal{O}({{n}^{3}})$. Цель настоящей работы состоит в исследовании возможности ускорения SVP путем замены оптимального проектора на быстро вычисляемый приближенный ${{\hat {\mathcal{P}}}_{k}}$.

3. ТЕОРИЯ СХОДИМОСТИ АЛГОРИТМА SVP С ПРИБЛИЖЕННЫМ ВЫЧИСЛЕНИЕМ ПРОЕКЦИЙ (ASVP)

Формализуем понятие приближенного проектирования на многообразие матриц ранга не выше $k$. Пусть оператор ${{\hat {\mathcal{P}}}_{k}}:{{\mathbb{R}}^{{m \times n}}} \to {{\mathbb{R}}^{{m \times n}}}$, заданный произвольным образом, всегда возвращает матрицу ранга не выше $k$. Будем называть ${{\hat {\mathcal{P}}}_{k}}$ оператором приближенного проектирования, если для любой матрицы $Y$ выполнено

с некоторой константой $\varepsilon > 0$.

По аналогии с алгоритмом SVP определим итерацию приближенного алгоритма ASVP равенством

Покажем, что для такого алгоритма ASVP справедлива теорема, аналогичная теореме 1.

Теорема 2. Пусть на матрице ${{X}_{*}}$ ранга не выше $k$ достигается минимум функционала

с оператором $\mathcal{A}$, для которого условие ограниченной изометриивыполняется на всех матрицах ранга не выше $2k$.

В этом случае для алгоритма ASVP с постоянным шагом $1{\text{/}}(1 + \sigma )$ и оператором ${{\hat {\mathcal{P}}}_{k}}$ приближенного проектирования, удовлетворяющим условию (4) с константой $\varepsilon $, справедливо неравенство

Доказательство. Используя арифметическое тождество

и символьную подстановку для произвольной матрицы $X$ можно записать где скалярное произведение $(U,V)$ для матриц $U,V \in {{\mathbb{R}}^{{m \times n}}}$ определяется соотношением $(U,V) = \operatorname{tr} ({{V}^{{\text{T}}}}U).$ Для произвольной матрицы $X$ ранга не выше $k$ воспользуемся во втором слагаемом (7) свойством ограниченной изометрии оператора $\mathcal{A}$ и установим верхнюю оценку для разности: Правую часть (8) обозначим через ${{f}_{t}}(X).$ Таким образом, и для всех матриц $X,$ ранг которых не превосходит $k.$

Значения ${{f}_{t}}(X)$ дают оценку на падение функционала невязки на всем множестве матриц ранга не выше $k,$ а вид ${{f}_{t}}(X)$ позволяет найти в явном виде матрицу ранга не выше $k,$ на которой падение невязки является значительным.

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

вместе с соответствием После замен ${{f}_{t}}(X)$ принимает вид Определяя ${{Y}_{{t + 1}}}$ равенством представим ${{f}_{t}}(X)$ в виде суммы: Так как второе слагаемое в (11) не зависит от $X,$ то минимум ${{f}_{t}}(X)$ на множестве всех матриц ранга не выше $k$ достигается на матрице ${{Z}_{{t + 1}}} = {{\mathcal{P}}_{k}}({{Y}_{{t + 1}}})$. При этом сама матрица ${{Y}_{{t + 1}}}$ формально совпадает с матрицей, полученной градиентным спуском из точки ${{X}_{t}}$ с величиной шага $\tau = 1{\text{/}}(1 + \sigma ),$ не зависящей от его номера $t$.

Видно, что ${{Z}_{{t + 1}}}$ совпадает с приближением на шаге $t + 1$ в алгоритме SVP.

Так как ${{Z}_{{t + 1}}}$ дает минимум ${{f}_{t}}(X)$, то значение ${{f}_{t}}({{Z}_{{t + 1}}})$ можно оценить сверху значением ${{f}_{t}}({{X}_{*}})$:

Используя свойство ограниченной изометрии оператора $\mathcal{A}$ на матрице $({{X}_{*}} - {{X}_{t}})$, ранг которой очевидным образом не превосходит $2k,$ преобразуем последнее выражение к виду Теперь воспользуемся формулой (7) для первого слагаемого, чтобы получить Положим ${{X}_{{t + 1}}} = {{\hat {\mathcal{P}}}_{k}}({{Y}_{{t + 1}}})$. В силу определения ${{\hat {\mathcal{P}}}_{k}}$ выполняется неравенство Рассмотрим разность ${{f}_{t}}({{X}_{{t + 1}}}) - {{f}_{t}}({{Z}_{{t + 1}}}).$ Подставляя для ${{f}_{t}}({{X}_{{t + 1}}})$ и ${{f}_{t}}({{Z}_{{t + 1}}})$ их представление из (11), сокращая постоянное слагаемое и учитывая (13), получаем Еще раз применим (11) теперь к правой части (14). Тогда где $\left\| {{{\mathcal{A}}^{{\text{T}}}}} \right\|$ – операторная норма ${{\mathcal{A}}^{{\text{T}}}}$, подчиненная второй норме в пространстве ${{\mathbb{R}}^{M}}$ и фробениусовой норме в пространстве ${{\mathbb{R}}^{{m \times n}}}$. Теперь мы готовы оценить разность $\psi ({{X}_{{t + 1}}}) - \psi ({{X}_{t}}).$ Действительно, как следует из (10), Заменяя ${{f}_{t}}({{Z}_{{t + 1}}})$ выражением (12), преобразуем последнее соотношение к окончательному виду Откуда после сокращения слагаемого $\psi ({{X}_{t}})$ в левой и правой частях неравенства заканчиваем доказательство.

Следствие 1. Пусть существует матрица ${{X}_{*}}$ ранга не выше $k$, для которой ${{\psi }_{\mathcal{A}}}({{X}_{*}}) = 0$ и выполняется условие ограниченной изометрии (3). Алгоритм ASVP с постоянным шагом $\tau = \tfrac{1}{{1 + \sigma }}$ сходится к решению ${{X}_{*}}$, а скорость сходимости определяется соотношением

Доказательство. Из предыдущей теоремы и в силу соотношений

имеем

В предположении $\sigma < 1{\text{/}}3,$ аналогичном тому, что используется для доказательства сходимости SVP в [7], видим, что $\varepsilon $ можно выбрать достаточно малой константой, чтобы из неравенства (16) следовала геометрическая сходимость ASVP.

Как следует из формулы (15), использование приближенного проектора ${{\hat {\mathcal{P}}}_{k}}$ в алгоритме ASVP ухудшает коэффициент линейной сходимости, что, конечно, не является неожиданным. Однако, как будет показано далее, выигрыш, получаемый за счет упрощения , превосходит проигрыш, связанный с увеличением числа итераций в ASVP.

В следствии 1 предполагается существование решения ${{X}_{*}}$, для которого ${{\psi }_{\mathcal{A}}}({{X}_{*}}) = 0$. Однако на практике для заданного вектора ${\mathbf{B}} \in {{\mathbb{R}}^{M}}$ минимум функционала может быть не равен нулю, например, в случае ${\mathbf{B}} = \mathcal{A}(X) + \Delta $, где $X$ – матрица малого ранга, а $\Delta $ имеет малую норму. Оценки для такого случая представлены в следствии 2.

Следствие 2. Пусть на матрице ${{X}_{*}}$ достигается минимум функционала ${{\psi }_{\mathcal{A}}}(X)$ среди всех матриц ранга не выше $k$, а $\Delta = {\mathbf{B}} - \mathcal{A}({{X}_{*}})$ имеет норму $\delta = \mathop {\left\| \Delta \right\|}\nolimits_2 $.

Если для приближения ${{X}_{t}}$ выполняется неравенство ${{\psi }_{\mathcal{A}}}({{X}_{t}}) \geqslant {{C}^{2}}\tfrac{{{{\delta }^{2}}}}{2}$ с константой $C > 0,$ то для следующего приближения, получаемого в методе ASVP, справедливо соотношение

с константой $\kappa $, удовлетворяющей неравенству где

Доказательство. Из теоремы 2

Учитывая, что ${{\delta }^{2}} \leqslant \tfrac{2}{{{{C}^{2}}}}\psi ({{X}_{t}})$, и раскрывая второе слагаемое, получаем Полагая

окончательно получаем

4. ПРИМЕНЕНИЕ АЛГОРИТМОВ SVP И ASVP К ЗАДАЧЕ ВОСПОЛНЕНИЯ МАТРИЦ МАЛОГО РАНГА

Рассмотрим применение алгоритма ASVP к задаче восполнения матриц малого ранга по случайному набору элементов (на случайном шаблоне). С точки зрения анализа сходимости алгоритмов SVP и ASVP задача восполнения сводится к специальному способу выбора аффинного преобразования $\mathcal{A}$.

Пусть ${{X}_{*}} \in {{\mathbb{R}}^{{m \times n}}}$ – неизвестная матрица ранга $k,$ причем $k \ll min(m,n),$ а $\Omega $ – случайный набор пар индексов $(i,j),$ где $i \in 1,2, \ldots ,m$ и $j \in 1,2, \ldots ,n$, причем любая пара $(i,j)$ может быть выбрана равновероятно среди всех $mn$ пар с вероятностью $q,$ которую мы определим позже. Обозначим мощность множества $\Omega $ как $M = {\text{|}}\Omega {\kern 1pt} {\text{|}}$, будем считать, что пары $(i,j) \in \Omega $ линейно упорядочены; будем также для простоты считать, что выполнено в точности $M = qmn$, а $q$ будем также называть плотностью известных элементов матрицы.

Определим линейный оператор $\mathcal{S}:{{\mathbb{R}}^{{m \times n}}} \to {{\mathbb{R}}^{M}}$ по правилу

где ${{(\mathcal{S}(X))}_{k}}$ – компонента вектора $\mathcal{S}(X)$, соответствующая паре $(i,j)$ в выбранном линейном упорядочении. Определим $B$ формулой $B = \mathcal{S}({{X}_{*}}) \in {{\mathbb{R}}^{M}}$ и зададим функционал ${{\psi }_{\mathcal{S}}}(X) = \tfrac{1}{2}\mathop {\left\| {\mathcal{S}(X) - {\mathbf{B}}} \right\|}\nolimits_2^2 $. Задача восполнения матрицы ${{X}_{*}}$ по части ее элементов записывается в виде Для упрощения записи в дальнейшем индекс $\mathcal{S}$ в обозначении функционала $\psi $ опускается.

В качестве оператора ${{\hat {\mathcal{P}}}_{k}}$ при этом можно использовать произвольный алгоритм построения малоранговой аппроксимации $\tilde {A} = {{\hat {\mathcal{P}}}_{k}}A$, достаточно близкой к проекции с использованием сингулярного разложения

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

При $\sigma \approx 0$ из следствия 2 подстановкой $\mathcal{A} = \mathcal{S}$, используя тривиальную оценку $\mathop {\left\| {{{\mathcal{S}}^{{\text{T}}}}} \right\|}\nolimits^2 \leqslant \tfrac{1}{q}$, мы получаем условие

а потому для сходимости метода восполнения матриц достаточно $\varepsilon = O(q)$.

Отметим, что для задачи восполнения матрицы малого ранга по значениям ее элементов на разреженном шаблоне $\Omega $ напрямую воспользоваться результатами разд. 2 и 3 нельзя.

Так, например, в случае произвольного шаблона $\Omega $ оператор $\mathcal{S}(X)$, определенный ранее, не удовлетворяет условию ограниченной изометрии, которое является ключевым при доказательстве сходимости метода. Тем не менее, как показано в [7], если $\Omega $ имеет случайный характер, то условие ограниченной изометрии для оператора $\mathcal{S}$ выполняется с вероятностью, неотличимой от $1,$ на почти всем множестве матриц ранга, не превосходящего $2k.$ Чтобы сформулировать соответствующее утверждение, нам потребуется понятие $\mu $-некогерентных матриц.

Определение 1. Пусть для матрицы $X \in {{\mathbb{R}}^{{m \times n}}}$ задано ее сингулярное разложение $X = U\Sigma {{V}^{{\text{T}}}}.$ Матрица $X$ называется $\mu $-некогерентной, если для матриц сингулярных векторов $U$ и $V$ справедливы неравенства

Используя определение 1, сформулируем следующий результат о выполнении свойства ограниченной изометрии для оператора $\mathcal{S}$.

Теорема 3 [Теорема 4.2 из [7]]. Существует константа $C \geqslant 0$ такая, что для любого $0 < \sigma < 1$, любых $\mu \geqslant 1$ и $n \geqslant m \geqslant 3$, и для любого шаблона разреженности $\Omega $, пары индексов которого выбираются случайно в соответствии с моделью Бернулли с параметром

определяющим вероятность для пары индексов $(i,j)$ принадлежать маске $\Omega ,$ оператор $\mathcal{S}$ удовлетворяет на множестве всех $\mu $-некогерентных матриц свойству ограниченной изометрии видас вероятностью не меньше $1 - exp( - nlog(n))$.

Из теорем 2 и 3 следует сходимость алгоритмов SVP и ASVP в случае, если на всех итерациях (т.е. для всех $t$) матрицы ${{X}_{{t + 1}}} - {{X}_{t}}$ и ${{X}_{t}} - {{X}_{ * }}$ обладают ограниченной константой некогерентности $\mu $. Отметим, что для разности ${{X}_{{t + 1}}} - {{X}_{t}}$ формальная проверка $\mu $-некогерентности является возможной и вычислительно недорогой процедурой. В то же время для разности ${{X}_{t}} - {{X}_{*}}$ на практике оценить константу некогерентности невозможно.

Из сказанного выше следует, что для задачи восполнения матриц малого ранга, заданных на случайном шаблоне, обоснование алгоритмов SVP и ASVP не является безусловным. Формальные реализации этих методов могут и проявляют нерегулярное поведение, т.е. расходятся при выборе теоретически обоснованного шага градиентного спуска $\tau $ (см., например, [4]) или в том случае, когда сингулярные числа решения ${{X}_{*}}$ быстро убывают.

На практике получение (быстро) сходящихся реализаций алгоритмов возможно только с использованием дополнительных приемов, среди которых мы выделяем два: последовательный набор ранга приближения и выбор шага $\tau $ в градиентном методе.

5. ВОЗМОЖНЫЕ МЕТОДЫ ПРИБЛИЖЕННОГО ПРОЕКТИРОВАНИЯ ASVP

6. СРАВНИТЕЛЬНЫЙ АНАЛИЗ ВЫЧИСЛИТЕЛЬНОЙ СЛОЖНОСТИ АЛГОРИТМОВ ASVP И SVP

Рассмотрим соотношения сложностей алгоритма SVP с использованием полного сингулярного разложения матрицы и алгоритма ASVP с использованием приближенного проектирования матрицы. Для анализа отметим следующие особенности алгоритмов.

1. Доказательство сходимости алгоритмов SVP и ASVP базируется на использовании свойства ограниченной изометрии (3) рассматриваемого оператора $A$.

2. В случае, если $A$ – оператор задачи матричного восполнения, существует асимптотическая оценка на минимальное количество известных элементов матрицы, которого достаточно для полного восполнения. Эта оценка определяется размерами неизвестной матрицы, рангом неизвестной матрицы, и “некогерентностью” неизвестной матрицы, численной величиной, характеризующей разреженность сингулярных факторов матрицы.

3. Не обязательно использовать именно такой порядок числа известных элементов матрицы, допустимо использовать и большие порядки.

4. В случае алгоритма ASVP гарантируется геометрическая сходимость алгоритма с порядком не менее $O\left( {\tfrac{\epsilon }{q}} \right)$, где $q$ – плотность известных элементов матрицы, а $\varepsilon $ – относительная ошибка приближенного проектирования. Таким образом, при уменьшении числа известных элементов матрицы требуется большая точность приближенного проектирования.

5. И в случае ASVP с использованием проектирования на случайные подпространства, и в случае ASVP с использованием псевдоскелетного метода, вводится дополнительный параметр $p$, определяющий некую вспомогательную размерность метода, причем с увеличением $p$ увеличивается и точность, и вычислительная сложность приближенного проектирования, и наоборот.

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

Проведем формальное сравнение алгоритмов. Сложность SVP и ASVP удобно выразить через параметры $m,\;n,\;k,\;q,\;p$. Кроме того, $q \geqslant {{q}_{{{\text{min}}}}}(m,n,k,\mu )$; также и для псевдоскелетного метода, и для случайного проектирования выполнено $\epsilon = O\left( {\tfrac{k}{p}} \right)$; считая множитель геометрической сходимости константой, чтобы обеспечить константный порядок числа итераций алгоритма, имеем $\epsilon = O(q)$, откуда $p = O\left( {\tfrac{k}{q}} \right)$.

Алгоритм SVP состоит из двух основных операций, оценим их сложность.

1. Вычисление невязки на маске известных элементов покомпонентно: $O(k)$ операций на каждый элемент, $O(mnkq)$ операций всего.

2. Вычисление сингулярного разложения для получения приближения следующей итерации: $O(mnmin(m,n))$.

Так как $q < 1$, $k < min(m,n)$, имеем сложность одной итерации алгоритма $O(mnmin(m,n))$, кубическую в случае квадратной матрицы.

Рассмотрим полученные выше оценки на сложность итерации алгоритма ASVP.

• $O((m + n){{p}^{2}}) + O(mnpq)$ для случайного проектирования,

• $O((m + n)kp) + O({{p}^{3}})$ для псевдоскелетного метода.

С учетом $q = O(\epsilon ) = O\left( {\tfrac{k}{p}} \right)$ для метода случайного проектирования имеем, что асимптотика сложности одной итерации псевдоскелетного метода в $O\left( {min\left( {m,n} \right){\text{/}}p} \right)$ раз меньше сложности одной итерации метода случайного проектирования.

Теперь рассмотрим сложности ASVP в зависимости от различных $q$. Пусть $n \geqslant m$. В случае минимально возможного $q = {{q}_{{{\text{min}}}}} = O\left( {\tfrac{{{{k}^{2}}{{\mu }^{2}}logn}}{m}} \right)$ имеем сложность порядка

• $O\left( {\tfrac{{(m + n){{m}^{2}}}}{{{{k}^{2}}{{\mu }^{4}}lo{{g}^{2}}n}}} \right) + O(mnk)$ для случайного проектирования;

• $O\left( {\tfrac{{{{m}^{3}}}}{{{{\mu }^{6}}{{k}^{3}}{{{\log }}^{3}}n}}} \right) + O\left( {\tfrac{{(m + n)m}}{{{{\mu }^{2}}logn}}} \right)$ для псевдоскелетного метода.

Покажем, что при больших порядках $q$ можно достичь и более низкой сложности ASVP. Пусть $q = {{m}^{{ - 1/2}}} \gg \tfrac{{{{k}^{2}}{{\mu }^{2}}logn}}{m}$ (при малых порядках $\mu $, близких к константе или к логарифму от размеров матрицы); тогда подстановкой имеем сложность

• $O((m + n)m{{k}^{2}})$ для случайного проектирования;

• $O((m + n){{m}^{{1/2}}}{{k}^{2}}) + O({{m}^{{3/2}}}{{k}^{3}})$ для псевдоскелетного метода.

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

Проведем экспериментальное сравнение алгоритмов SVP и ASVP c двумя рассмотренными вариантами методов приближенного проектирования.

• Рассматриваются квадратные матрицы размера $m = n = 1000$. Рассматриваются ранги матриц $k$ от $5$ до $25$ и плотности распределения известных элементов $q$ от $0.1$ до $0.4$.

• Матрицы формируются случайно следующим образом.

– Формируются случайные матрицы ${{\hat {U}}_{k}} \in {{\mathbb{R}}^{{m \times k}}}$ и ${{\hat {V}}_{k}} \in {{\mathbb{R}}^{{n \times k}}}$, каждый элемент которых распределен случайно по стандартному нормальному распределению.

– Вычисляются ортогональные базисы ${{U}_{k}}$ и ${{V}_{k}}$ в пространствах столбцов ${{\hat {U}}_{k}}$ и ${{\hat {V}}_{k}}$ соответственно, например, с помощью $QR$-разложения. Экспериментально установлено, что некогерентность факторов ${{U}_{k}},\;{{V}_{k}}$, выбранных таким образом, растет логарифмически от $m,\;n$.

– Искомая матрица ${{X}_{*}}$ задается в виде ${{X}_{*}} = {{U}_{k}}{{\Sigma }_{k}}V_{k}^{{\text{T}}}$ с матрицей ${{\Sigma }_{k}}$ сингулярных значений ${{\sigma }_{1}} \geqslant {{\sigma }_{2}} \geqslant {{\sigma }_{3}} \geqslant \ldots $ . Экспериментально установлено, что сходимость методов SVP, ASVP зависит от характера убывания сингулярных чисел ${{\sigma }_{j}}$; в экспериментах рассмотрены следующие законы ${{\sigma }_{j}}$: ${{\sigma }_{i}} = 1$ (фиг. 1), ${{\sigma }_{i}} = 1{\text{/}}i$ (фиг. 2), ${{\sigma }_{i}} = 1{\text{/}}{{i}^{2}}$ (фиг. 3), ${{\sigma }_{i}} = 1{\text{/}}{{2}^{i}}$ (фиг. 4).

• Сравнение проводится с ограничением по времени практической работы алгоритмов. Ограничение задается константой, которая определяется временем выполнения десяти или пятидесяти итераций алгоритма SVP с полным SVD; предполагается, что за это время приближенный алгоритм ASVP выполнит большее число итераций.

Ниже приведены графики зависимостей невязки от $q$ и $k$ при фиксированном размере матрицы в случае алгоритма SVP и двух рассмотренных вариантов алгоритма ASVP при различных законах падения сингулярных чисел неизвестной матрицы. Значения полных относительных ошибок по всей матрице на практике во всех экспериментах отличаются от относительных невязок на маске известных элементов не более, чем в два раза.

Графики сходимости в случае полного SVP приведены в двух вариантах: с шагом $q$ (фиг. 1д–4д) и c шагом $3{\text{/}}4$ (фиг. 1е–4е). Так как рассматривалось ограничение в 50 итераций SVP, вариант SVP с большим шагом на большой части графиков не достиг сходимости, так как такого числа итераций оказалось недостаточно для завершения последовательного набора ранга.

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

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