Report.tex

\documentclass[oneside,final,14pt]{extreport} % class of document
\usepackage[utf8]{inputenc} % codir and russian
\usepackage[russian]{babel} % russian
\usepackage{vmargin} % format of document:
\setpapersize{A4} % size
\setmarginsrb{2cm}{1.5cm}{1cm}{1.5cm}{0pt}{0mm}{0pt}{13mm} % fields
\usepackage{indentfirst} % gap before first paragraphe
\sloppy % not to go through fields!
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{mathrsfs}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[ruled,vlined]{algorithm2e}
\begin{document}

\begin{titlepage}
	\centering{Московский государственный университет \\}
	\centering{Факультет вычислительной математики\\}
	\centering{и кибернетики}
	\vskip10cm
	\centering{\bf \Large Методы оптимизации для низкоранговых матриц}
	\vskip5cm
	\begin{flushright}
		{\bf Автор:} Фролов Антон
	\end{flushright}
	\vfill
	\centering{Москва\\}
	\centering{2021}
\end{titlepage}
\setcounter{page}{2} 

\tableofcontents
\chapter{Введение}
\section{О задаче}

Оптимизация — это задача нахождения экстремума (минимума или максимума) целевой функции в некоторой области конечномерного пространства, ограниченной набором линейных и/или нелинейных равенств и/или неравенств.

В нашем случае, мы имеем дело с задачей оптимизации с ограничением на ранг или задачей RCOP (Rank Constrained Optimization Problem). Она заключается в оптимизации выпуклой целевой функции от прямоугольной матрицы заданного размера с учетом ограничения на её ранг. 

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

Существующие подходы к решению RCOP в основном сосредоточены на методах на основе попеременного проецирования и комбинированных алгоритмах линеаризации и факторизации с применением факторного анализа. Однако эти итерационные подходы зависят от первоначального предположения, и их быстрая сходимость не может быть гарантирована. Кроме того, был предложен метод, подобный методу Ньютона, применяемый в задаче стабилизации с обратной связью. Метод оптимизации с помощью римановых многообразий был применен для решения матричных уравнений Ляпунова большой размерности путем нахождения приближения низкого ранга.

\section{Формулировка задачи}

Задача оптимизации (миннимизации) с ограничением на ранг может быть сформулирована так

\begin{equation}
	\label{new:PF}
	\begin{aligned}
		&\min_X f(X) \\ &s.t. \ X \in \mathbb{R}^{m \times n}, \ rank(X) \le r
	\end{aligned}
\end{equation}

где $f(X)$ выпуклая функция и $\mathbb R^{m \times n}$ -- множество вещественных матриц размера $m$ на $n$. Не ограничивая общности, будем считать, что  $m \le n$. 

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


\textbf{Лемма 1.} 
{\it Пусть дана матрица $X \in \mathbb R^{m \times n}$. Тогда неравенство $rank(X) \le r$ выполнено тогда и только тогда, когда существуют матрицы $Y =Y^T \in R^{m \times m}$ и $Z=Z^T \in R^{n \times n}$ такие, что

\begin{equation}
 rank(Y) + rank(Z) \le 2r, \begin{bmatrix} Y && X \\ X^T && Z \end{bmatrix} \succeq 0  
\end{equation}
}

\textbf{Утверждение 2.}
{\it Равенство $Z = X^T X$ эквивалентно неравенству $rank\left( \begin{bmatrix} I_m && X \\ X^T && Z \end{bmatrix} \right) \le m$, где $Z \in \mathbb S^n, X \in \mathbb R^{m \times n}$ и $I_m \in \mathbb R^{m \times m}$ единичная матрица.}

\textit{Доказательство:} С учетом того, что ранг симметричной блочной матрцы равен рангу блока, стоящего на диагонали, плюс ранг его дополнения Шура (Guttman rank additivity formula), имеем следующее отношение, $rank\left( \begin{bmatrix} Y && X \\ X^T && Z \end{bmatrix} \right) \le m \Leftrightarrow rank(I_m) + rank(Z - X^TX) \le m \Leftrightarrow m + rank(Z-X^TX) \le m \Leftrightarrow rank (Z - X^TX) = 0 \Leftrightarrow Z = X^TX$.


%\textbf{Замечание}: Когда $Z=X^T X$, это говорит о том, что $\begin{bmatrix} I_m && X \\ X^T && Z \end{bmatrix} \succeq 0$ верно для $I_m \succeq 0$ и ее дополнение Шура есть нулевая матрица.

Далее приводится расширенная лемма о полупопределенном вложении.

\textbf{Лемма 3.}
{\it Пусть дана матрица $X \in \mathbb{R}^{m \times n}$. Тогда неравенство $rank(X) \le r$ выполнено тогда и только тогда, когда существует матрица $Z \in \mathbb S^n$ такая, что
\begin{equation}
 rank(Z) \le r \ \& \ rank \left( \begin{bmatrix} I_m && X \\ X^T && Z \end{bmatrix} \right) \le m  
\end{equation}
}

\textit{Доаказательство:} Из утверждения (2), $rank \left( \begin{bmatrix} I_m && X \\ X^T && Z \end{bmatrix} \right) \le m \Leftrightarrow Z = X^T X$. Известно, что $rank(Z) = rank(X)$, когда $Z = X^T X$. Таким образом, доказательство завершено.

%\textbf{Замечание}: Лемму 3 можно интерпретировать, как один из случаев Леммы 1, полагая $Y$ в Лемме 1 равным $I_m$ и $rank(Z) = rank(X)$ и заменяя $2r$ на $m+r$. Лемма 3 в этом случае переводит неравенство для ранга прямоугольной матрицы в два неравенства для рангов полуопределенных матриц, в то время как лемма 1 эквивалентно переводит его в одно неравенство для ранга и одно условие на полуопределенность.

Следовательно, основываясь на Лемме 3 задача \ref{new:PF} с ограничением на ранг прямоугольной матрицы может быть эквивалентным образом переведена в следующую формулировку с ограничениями на ранги положительно полуопределенных матриц

\begin{equation}
	\begin{aligned}
		&\min_{X,Z} f(X) \\ &X \in \mathbb{R}^{m \times n}, rank(Z) \le r \\ &rank \left( \begin{bmatrix} I_m && X \\ X^T && Z \end{bmatrix} \right) \le m
	\end{aligned}
\end{equation}
где $Z \in \mathbb S^n$ - дополнительно введенная матрица.
%Для простоты следующее обсуждение RCOP рассматривает ранговые ограничения только для полуопределенных матриц.

\section{Сингулярное разложение}
Сингулярным разложением или SVD (Singular Vector Decomposition) матрицы $M$ порядка $m \times n$ называется разложение следующего вида
$$M=U\Sigma V^{*}$$
где $\Sigma$ -- матрица размера $m \times n$ с неотрицательными элементами, у которой элементы, лежащие на главной диагонали — это сингулярные числа (а все элементы, не лежащие на главной диагонали, являются нулевыми), а матрицы $U$ (порядка $m$) и $V$ (порядка $n$) -- это две унитарные матрицы, состоящие из левых и правых сингулярных векторов соответственно (а $V^*$ — это сопряжённо-транспонированная матрица к $V$).


В некоторых практических задачах требуется приближать заданную матрицу $M$ некоторой другой матрицей $M_{k}$ с заранее заданным рангом $k$. Известна следующая теорема, которую иногда называют теоремой Эккарта — Янга.

Если потребовать, чтобы такое приближение было наилучшим в том смысле, что норма Фробениуса разности матриц $M$ и $M_k$ минимальна, при ограничении $rank(M_k) = k$, то оказывается, что наилучшая такая матрица $M_k$ получается из сингулярного разложения матрицы $M$ по формуле:

$$ M_k = U \Sigma_kV^*$$
где $\Sigma_K$ -- матрица $\Sigma$, в которой заменили нулями все диагональные элементы, кроме $k$ наибольших элементов.
Если элементы матрицы $\Sigma$  упорядочены по невозрастанию, то выражение для матрицы $M_{k}$ можно переписать в такой форме:
$$ M_k = U_k \Sigma_kV_k^*$$
где матрицы $U_k, \Sigma _{k}$ и $V_{k}$ получаются из соответствующих матриц в сингулярном разложении матрицы $M$ обрезанием до ровно $k$ первых столбцов.

Таким образом видно, что приближая матрицу $M$ матрицей меньшего ранга, мы выполняем своего рода сжатие информации, содержащейся в $M$: матрица $M$ размера $m \times n$ заменяется меньшими матрицами размеров $m \times k$ и $k \times n$ и диагональной матрицей с $k$ элементами. При этом сжатие происходит с потерями — в приближении сохраняется лишь наиболее существенная часть матрицы $M$.

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

Сингулярное разложение и теорема SVD понадобятся в дальнейшем при рассмотрении одного из методов оптимизации.

\chapter{Методы оптимизации}
\section{Попеременная оптимизация}


В подходе попеременной минимизации (AltMin -- Alternating Minimization) искомая матрица $X$ низкого ранга записывается в билинейной форме, то есть $X = UV^T$. Матрицы $U$ и $V$ имеют $r$ столбцов и $m$ и $n$ строк соответсвенно, поэтому ранг матрицы $UV^T$ не превосходит $r$.

Тогда задачу можно переписать в виде
\begin{equation}
	\begin{aligned}
		&\min_{U,V} f(UV^T) \\ &s.t. \ U \in \mathbb{R}^{m \times r}, \  V \in \mathbb{R}^{n \times r} 
	\end{aligned}
\end{equation}

Алгоритм начинается с инициализации матриц $U_0$ и $V_0$ и постепенно минимизирует целевую функцию, фиксируя одну из матриц и минимизируя функцию по второй.

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

\begin{algorithm}
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\Input{T -- количество итераций, $r$ -- ограничение на ранг}
\Output{$U \in \mathbb{R}^{m \times r}, V \in \mathbb{R}^{n \times r}$}
\textbf{begin}
\begin{enumerate}
	\item \textbf{Initialize} Инициализировать матрицу $U_0$
	\item \textbf{for} $t = 1,2,\ldots, T$:
	\item $V_t = \text{argmin}_{V \in \mathbb{R}^{n \times r}} f(U_{t-1} V^T)$
	\item $U_t = \text{argmin}_{U \in \mathbb{R}^{m \times r}} f(U V_t^T)$
	\item\textbf{end for}
\end{enumerate}
\textbf{end}
\caption{Алгоритм попеременной минимизации (AltMin)}
\end{algorithm}

\section{IRM}
Далее рассмотрим алгоритм итеративной минимизации ранга (IRM), также основанный на попеременной минимизации.

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

Сходимость метода IRM доказывается с помощью теории двойственности и условий Каруша-Куна-Таккера, доказательство опустим.

Учитывая вышесказанное, сформулируем задачу RCOP следующим образом,

\begin{equation}
	\label{new:IRM}
	\begin{aligned}
		&\min_X f(X) \\ &s.t. \  X \in \mathbb{S}_+^n, \ rank(X) \le r
	\end{aligned}
\end{equation}

Количество ненулевых собственных значений матрицы совпадает с ее рангом. Для неизвестной квадратной матрицы $U \in \mathbb S_+^n$ невозможно найти ее собственные значения пока она не определена. Обратим внимание на тот факт, что при ранге матрицы $r$ она имеет $r$ ненулевых собственных значений. Поэтому вместо ограничения ранга мы сосредоточимся на ограничении собственных значений $U$ таким образом, чтобы все $n - r$ собственных значений $U$ были нулями. Далее приведем необходимые утверждения и леммы, которые будут использоваться в дальнейшем.

\textbf{Утверждение 4.}
{\it $(r+1)$е по величине собственное значение $\lambda_{n-r}$ матрицы $U \in S_+^n$ не больше, чем $e$ тогда и только тогда, когда $eI_{n-r} - V^T U V \succeq 0$, где $I_{n-r}$ единичная матрица размера $n - r$, $V \in R^{n \times (n-r)}$ -- матрица из собственных векторов, отвечающих $n - r$ наименьшим собственным значениям $U$.
}

\textit{Доказательство:} Предположим, что собственные значения $U$ отсортированы по убыванию в виде $[ \lambda_n, \lambda_{n-1}, \ldots , \lambda_1 ]$. Так как отношение Рэлея собственного вектора это собственное значение, которому он отвечает, то $eI_{n-r} - V^T U V$ диагональная матрица с диагональными элементами $[e-\lambda_{n-r}, e-\lambda_{n-r-1}, \ldots, e-\lambda_1 ]$. Следовательно, $e \ge \lambda_{n-r}$ тогда и только тогда, когда $eI_{n-r} - V^TUV \succeq 0$.

\textbf{Следствие 5.}
{\it Если $e = 0$ и $U$ положительно полуопределенная матрица, неравенство $rank(U) \le r$ выполнено тогда и только тогда, когда $eI_{n-r} - V^TUV \succeq 0$.
%Равенство $rank(U) = r$ выполнено тогда и только тогда, когда $eI_{n-r} - V^TUV \succeq 0$ и $V_p^TUV_p \succ 0$. Матрицы $V \in R^{n \times (n-r)}$ и $V_p \in R^{n \times 1}$ состоят из собственных векторов отвечающих соответственно $n-r$ наименьшим собственным значениям и $n - r + 1$ наименьшему собственному значению $U$.
}

Однако для задачи \ref{new:IRM} до нахождения матрицы $X$ мы не можем получить матрицу $V$, т.е. собственные векторы $X$. Поэтому предлагается итерационный метод решения \ref{new:IRM} путем постепенного приближения к заданному рангу. На каждом шаге $k$ мы будем решать следующую задачу полуопределенного программирования, сформулированную ниже

\begin{equation}
	\label{new:IRM2}
	\begin{split}
		\min_{X_k,e_k} f(X_k) + \omega_k e_k \\ X_k \in \mathbb{S}_+^n \\e_kI_{n-r} - V_{k-1}^T X_kV_{k-1} \succeq 0 \\ e_k \le e_{k-1}	
	\end{split}
\end{equation}

где $w_k = w_0t^k$ - весовой коэффициент на итерации $k$. $w_k$ растет с увеличением $k$, когда заданы параметры $t > 1$ и $w_0 > 0$. $V_{k-1} \in R^{n \times (n - r)} $ -- ортонормированные собственные векторы, отвечающие $n - r$ наименьшим собственным значениям матрицы $X_{k-1}$, найденной на предыдущей итерации $k - 1$. На первой итерации, когда $k = 1$, $e_0$ является $(n - r)$-м наименьшим собственным значением матрицы $X_0$, описанной ниже. Известно, что решатель задач полуопределенного программирования (SDP -- Semi-Definite Programming), основанный на методе внутренней точки, имеет вычислительную временную сложность порядка $O(m (n ^ 2m + n ^ 3))$ для задачи SDP с $m$ линейными ограничениями и линейным матричным неравенством размерности $n$. В результате дополнительное линейное матричное неравенство $e_kI_{n-r} - V_{k-1}^T X_kV_{k-1} \succeq 0$ в подзадаче приводит к дополнительным линейным ограничениям порядка выше $O(n^2)$. Следовательно, временная сложность каждой итерации ограничена снизу величиной $O(n ^ 6)$.

На каждом шаге мы пытаемся оптимизировать исходную целевую функцию и в то же время минимизировать параметр $e$, чтобы при $e = 0$ выполнялось ограничение на ранг матрицы $X$. Весовой коэффициент, $w_k$, действует как коэффициент регуляризации, и увеличение его значения на каждом шаге приводит к постепенному уменьшению $(r + 1)$-го наибольшего собственного значения до нуля. Благодаря этому коэффициенту регуляризации наш алгоритм не просто переключается между поиском неизвестной матрицы и ее собственными векторами, он также приводит к быстрому уменьшению фиктивной переменной $e_k$ до нуля. Кроме того, дополнительное ограничение $e_k \le e_{k-1}$ гарантирует, что $e_k$ монотонно убывает. Вышеупомянутый подход повторяется до тех пор, пока не будет выполнено $e \le \varepsilon$, где $\varepsilon$ - малый порог для критерия останова. Каждая итерация решается с помощью решателя задач SDP на основе метода внутренней точки, который применим к задачам SDP малого и среднего размера. Несложно расширить \ref{new:IRM2} на задачи с многоранговыми ограничениями. Для краткости здесь описана простейшая версия с ограничением одного ранга.

Кроме того, на первой итерации $k = 1$ требуется начальная матрица $V_0$. Интуитивно понятно что, нужно использовать решение задачи \ref{new:IRM2} с ослабленными ограничениями, отбросив последнее ограничение и штрафной член целевой функции. При этом предположении $X_0$ находится через

\begin{equation}
	\label{new:IRM3}
	\begin{aligned}
	& \min_{X_0} f(X_0) \\
	& s.t. \ X_0 \in \mathbb{S}_+^n
	\end{aligned}
\end{equation}

Тогда $V_0$ -- собственные векторы, соответствующие $n - r$ наименьшим собственным значениям $X_0$.

\begin{algorithm}
\SetKwInOut{Input}{Input}\SetKwInOut{Output}{Output}
\Input{$w_0 , t, \varepsilon_1 , \varepsilon_2, k_{max}$}
\Output{матрица $X^*$, на которой достигается минимум}
\textbf{begin}
\begin{enumerate}
	\item \textbf{Initialize:} Положить $k = 0$, решить задачу \ref{new:IRM3} для получения матрицы $V$ из $X$ через спектральное разложение. Положить $k=k+1$
	\item \textbf{while} $k \le k_{max} \ \& \ e_k \ge \varepsilon_1 \ \|\  |(f(X_k)-
f(X_{k-1}))/f(X_{k-1})| \ge \varepsilon_2$
	\item Решить последовательную задачу \ref{new:IRM2} и найти $X_k, e_k$
	\itemОбновить $V_k$ из $X_k$ с помощью спектрального разложения и положить $k=k+1$
	\itemОбновить $w_k$ с помощью формулы $w_k=w_{k-1} \ast t$
	\item\textbf{end while}
\end{enumerate}
\textbf{end}
\caption{Последовательая минимизация ранга для решения задачи \ref{new:IRM}}
\end{algorithm}


\section{Алгоритм GECO}

В этом разделе рассмотрим жадный алгоритм покомпонентной оптимизации или GECO (Greedy Efficient Component Optimization).

Далее будем искать решение задачи оптимизации вида:

\begin{equation}
  \label{GECO PF}
  \begin{aligned}
  & \min_A R(A) \\
  & s.t. \ rank(A) \le r
  \end{aligned}
\end{equation}
где $R: \mathbb{R}^{m \times n} \rightarrow \mathbb{R}$ - выпуклая и гладкая функция.


Алгоритм GECO основан на простом, но действенном наблюдении: вместо представления матрицы $A$ с использованием чисел $m \times n$ мы представляем ее с помощью бесконечномерного вектора $\lambda$, индексированного всеми парами $(u, v)$ взятыми из единичных сфер $\mathbb{R}^m$ и $\mathbb{R}^n$ соответственно. В этом представлении низкий ранг соответствует разреженности вектора $\lambda$.

Таким образом, мы можем свести задачу, заданную в уравнении \ref{GECO PF}, к задаче минимизации вектор-функции $f (\lambda)$ по множеству разреженных векторов, $\|\lambda\|_0 \le r$. Основываясь на этой редукции, мы применяем алгоритм жадной аппроксимации для минимизации выпуклой векторной функции с учетом ограничения разреженности. На первый взгляд, прямое применение этой редукции кажется невозможным, поскольку $\lambda$ - бесконечномерный вектор, и на каждой итерации жадного алгоритма нужно искать по бесконечному набору координат $\lambda$. Однако эту задачу поиска можно представить как задачу нахождения первых ведущих правых и левых сингулярных векторов матрицы.

Пусть $A \in R^{m \times n}$ -- матрица размера $m$ на $n$, и, не ограничивая общности, предположим, что $ m \le n $. Теорема SVD утверждает, что матрицу $A$ можно записать в виде $ A = \sum_{i = 1}^m \lambda_iu_iv_i^T$, где $u_1, \ldots, u_m$ лежат во множестве $\mathcal{U} = \{u \in R^m: \| u \| = 1 \}$, $v_1, \ldots, v_m$ -- во множестве $ \mathcal{V} = \{v \in R^n: \| v \| = 1 \}$, а $ \lambda_1, \ldots, \lambda_m$ - скаляры. Чтобы упростить представление, мы предполагаем, что каждое действительное число представлено с использованием конечного числа битов, поэтому наборы $\mathcal{U}$ и $\mathcal{V}$ являются конечными наборами. Отсюда следует, что мы также можем записать $A$ в виде $A = \sum_{(u, v) \in \mathcal{U} \times \mathcal{V}} \lambda_{u, v} uv^T$, где $\lambda \in R^{ | \mathcal{U} \times \mathcal{V} |}$, и мы индексируем элементы $ \lambda $, используя пары $ (u, v) \in \mathcal{U} \times \mathcal{V} $. Обратим внимание, что представление $ A $ с помощью вектора $ \lambda $ не единственно, но из теоремы SVD всегда существует представление $ A $, для которого количество ненулевых элементов $ \lambda $ не больше чем $ m $, т.е. $ \| \lambda \|_0 \le m $, где $ \| \lambda \|_0 = | \{(u, v): \lambda_ {u, v} \ne 0 \} | $. Кроме того, если $ rank (A) \le r $, то существует представление $ A $ с помощью вектора $ \lambda $, для которого $ \| \lambda \| _0 \le r $

Для (разреженного) вектора $ \lambda \in R^{| \mathcal{U} \times \mathcal{V} |} $ мы определяем соответствующую матрицу как

\begin{equation}
  A(\lambda) = \sum_{(u,v) \in \mathcal{U} \times \mathcal{V}} \lambda_{u,v} u v^T
\end{equation}

Обратим внимание, что $ A (\lambda) $ -- линейное отображение. Для функции $ R: R^{m \times n} \rightarrow R $ мы определяем функцию

\begin{equation}
  f(\lambda) = R(A(\lambda)) = R \left( \sum_{(u,v) \in \mathcal{U} \times \mathcal{V}} \lambda_{u,v} u v^T \right)
\end{equation}

Легко проверить, что если $ R $ -- выпуклая функция над $ R^{m \times n} $, то $ f $ выпукла над $ R ^ {| \mathcal{U} \times \mathcal{V} |} $ (поскольку $ f $ -- композиция $ R $ над линейным отображением). Таким образом, мы можем свести задачу, заданную в уравнении \ref{GECO PF}, к задаче

\begin{equation}
	\label{GECO MIN}
	\begin{aligned}
  	& \min_{\lambda} f(\lambda) \\
  	& s.t. \ \lambda \in \mathbb{R}^{|U \times R|}, \| \lambda \|_0 \le r
	\end{aligned}
\end{equation}

Хотя задача оптимизации, заданная в уравнении \ref{GECO MIN}, относится к произвольному большому пространству, прямая жадная процедура выбора может быть реализована эффективно. В начале жадного алгоритма полагается $ \lambda = (0, \ldots, 0) $. На каждой итерации мы сначала находим векторы $ (u, v) $, которые максимизируют величину частной производной $ f (\lambda $) по $ \lambda_ {u, v} $. Предполагая, что $ R $ дифференцируема, и используя цепное правило, получаем:

\begin{equation}
  \frac{\partial f(\lambda)}{\partial \lambda _{u,v}} = \langle \nabla R(A(\lambda)), uv^T \rangle = u^T \nabla R(A(\lambda))v
\end{equation}

где $ \nabla R (A (\lambda)) $ - матрица размера $ m \times n $ частных производных $ R $ по элементам матрицы $ A (\lambda) $. Векторы $ u, v $, которые максимизируют величину вышеуказанного выражения, являются левым и правым сингулярными векторами, отвечающими наибольшему сингулярному значению $ \nabla R (A (\lambda)) $. Следовательно, даже несмотря на то, что количество элементов в $ \mathcal{U} \times \mathcal{V} $ очень велико, мы все же можем эффективно выполнить жадный выбор одной пары $ (u, v) \in \mathcal{U} \times \mathcal{V} $.

В некоторых ситуациях даже вычисление ведущих сингулярных векторов может оказаться слишком дорогостоящим. Поэтому допустима приближенная максимизация. Обозначим через ApproxSV ($ \nabla R (A (\lambda)), \tau $) функцию, которая возвращает векторы, для которых

\begin{equation}
  u^T \nabla R(A(\lambda))v \ge (1 - \tau) \max_{p,q} p^T\nabla R(A(\lambda))q
\end{equation}

Пусть $ U $ и $ V $ - матрицы, столбцы которых содержат векторы $ u $ и $ v $, которые мы получили до этого момента. Второй шаг каждой итерации алгоритма устанавливает $ \lambda $ как решение следующей задачи оптимизации:

\begin{equation}
	\label{GECO MIN2}
  \min_{\lambda \in \mathbb R^{|\mathcal{U} \times \mathcal{V}|}} f(\lambda)  \ \text{s.t.} \ \text{supp}(\lambda) \subseteq \text{span}(U) \times \text{span}(V)
\end{equation}

где $ \text{supp} (\lambda) = \{(u, v): \lambda_{u, v} \ne 0 \} $ и $ \text {span} (U), \ \text {span} (V) $ - линейные оболочки столбцов  $ U $ и $ V $ соответственно. 

Далее опишем, как решить уравнение \ref{GECO MIN2}. Пусть $ s $ будет количеством столбцов в $ U $ и $ V $. Обратим внимание, что любой вектор $ u \in \text {span} (U) $ можно записать как $ Ub_u $, где $ b_u \in \mathbb R^s $, и аналогично любой $ v \in \text {span} (V) $ можно записать как $ Vb_v $. Следовательно, если носитель функции $ \lambda $ находится в $ \text {span} (U) \times \text {span} (V) $, получаем, что $ A (\lambda) $ может быть записано как

\begin{equation}
  A(\lambda) = \sum_{(u,v) \in \text{supp}(\lambda)} \lambda_{u,v}(Ub_u)(Vb_v)^T = U \left( \sum_{(u,v) \in \text{supp}(\lambda)} \lambda _{u,v}b_ub_v^T\right )V^T
\end{equation}

Таким образом, любая $ \lambda $, носитель функции которой находится в $ \text {span} (U) \times \text {span} (V) $, дает матрицу $ B (\lambda) = \sum_ {u, v} \lambda _ {u, v} b_ub_v^T $. Теорема SVD говорит, что обратное утверждение также верно, а именно, для любого $ B \in \mathbb R^{s \times s} $ существует $ \lambda $, носитель функции которого находится в $ \text {span} (U ) \times \text {span} (V) $, который порождает $ B $ (а также $ UBV^T $). Обозначим $ \tilde R (B) = R (UBV^T) $, из этого следует, что уравнение \ref{GECO MIN2} эквивалентно следующей задаче оптимизации без ограничений $ \min_{B \in \mathbb R^{s \times s}} \tilde R (B) $. Легко проверить, что $ \tilde R $ - выпуклая функция, и поэтому ее можно эффективно минимизировать. Как только мы получим матрицу B, которая минимизирует $ \tilde R (B) $, мы можем использовать ее SVD для нахождения соответствующей $ \lambda $.

На практике нам не нужно подстраивать $ \lambda $, а лишь подстраивать такие матрицы $ U, V $, что $ A (\lambda) = U V^T $.

\begin{algorithm}

1: \textbf{Ввод:} Выпукло-гладкая функция $R : R^{m \times n} \rightarrow R$; ограничение на ранг $r$; константа $\tau \in [0, 1/2]$

2: \textbf{Initialize:} $U = [], V = []$

3: \textbf{for} i=1,...,r \textbf{do}

4: $(u,v)$=ApproxSV($\nabla R(UV^T),\tau$)

5: \quad Положить $U = [U,u]$ and $V = [V,v]$

6: \quad Положить $B = \text{argmin}_{B : \in \mathbb R^{i \times i}} R(UBV^T)$

7: \quad Вычислить SVD: $B = PDQ^T$

8: \quad Обновить: $U=UPD,V=VQ$

9: \textbf{end for}
\caption{GECO}
\end{algorithm}


Время выполнения алгоритма следующее. Шаг 4 можно выполнить за время $ O (N \log (n) / \tau) $, где $ N $ - количество ненулевых элементов $ \nabla R (UV^T) $, используя степенной метод. Поскольку анализ алгоритма позволяет $ \tau $ быть константой (например, 1/2), это означает, что время выполнения равно $ O (N \log (n)) $. Время выполнения шага 6 зависит от структуры функции $ R $. Наконец, время выполнения шага 7 составляет не более $ r^3 $, шаг 8 можно выполнить за время $ O (r^2 (m + n)) $.

\chapter{Сравнение методов}

Методы оптимизации для низкоранговых матриц применимы для многих задач, например:
\begin{itemize}
  \item Приближение матрицей низкого ранга
  \item Заполнение матрицы
  \item Идентификация систем
  \item Обработка сигналов
  \item Редукция моделей
\end{itemize}

Протестируем первые два метода на второй задаче, то есть на задаче заполнения матрицы.

Эта задача может быть сформулирована в виде:

\begin{equation}
	\begin{aligned}
		&\min_X \sum_{(i,j) \in E} (Y_{i,j} - X_{i, j})^2 \\ &s.t. \ rank(X) \le r
	\end{aligned}
\end{equation}
где $E$ -- множество пар индексов известных значений данной матрицы $Y \in \mathbb{R}^{m \times n}$.

Для решения задачи минимизации в алгоритме AltMin будем использовать функцию minimize модуля optimize библиотеки scipy.

Для реализации выпуклого программирования в алгоритме IRM воспользуемся библиотекой cvxpy. В качестве решателя возьмём решатель MOSEK, который проявил себя лучше других в ходе исследования.

В качестве заполняемой матрицы возьмем случайную матрицу $Y$ размера $m$ на $n$, состоящую из целых чисел от 0 до 10. Для того, чтобы решить, какие значения мы знаем, возьмем булевскую матрицу того же размера, в которой на каждой из позиций стоит 1 с вероятностью $p$ и 0 с вероятностью $1-p$.

Проведем эксперименты для матрицы размера 15 на 20, которую будем приводить к рангу 10, при $p = 0.8$ и $p = 0.9$ и для матрицы размера 30 на 45, которую будем приводить к рангу 20, при $p = 0.8$.

Для первого и второго эксперимента положим параметры $\omega_0$ и $t$ алгоритма IRM равными 0.1 и 1.5 соответственно. Для третьего эксперимента положим $\omega_0 = 1$, $t = 1.1$.

В ходе исследования будем отслеживать средний квадрат ошибки и время выполнения итераций.

В случае алгоритма IRM также будем отслеживать величину $e_k$, которая ограничивает сверху $n-r$-е наименьшее собственное значение.

Метод AltMin устроен так, что с первой итерации его решение имеет заданный ранг, поэтому приводить для него $n-r$-е собственное значение не имеет смысла.

\begin{table}[ht]
	\centering
	\begin{tabular}{|c|c|c|}
		\hline Iteration & MSE & Time \\ \hline
		1 & 0.4877 & 2.21 \\ \hline
		2 & 0.1795 & 2.23 \\ \hline
		3 & 0.1161 & 2.03 \\ \hline
		4 & 0.0919 & 2.06 \\ \hline
		5 & 0.0786 & 1.96 \\ \hline
		6 & 0.0697 & 1.81 \\ \hline
		7 & 0.0632 & 1.66 \\ \hline
		8 & 0.0580 & 1.63 \\ \hline
		9 & 0.0537 & 2.01 \\ \hline
		10 & 0.0500 & 2.09 \\ \hline
		11 & 0.0468 & 2.06 \\ \hline
		12 & 0.0440 & 1.92 \\ \hline
		13 & 0.0417 & 1.98 \\ \hline
		14 & 0.0397 & 1.91 \\ \hline
		15 & 0.0381 & 1.64 \\ \hline
	\end{tabular}
	\begin{tabular}{|c|c|c|}
		\hline Iteration & MSE & Time \\ \hline
		1 & 1.0135 & 1.67 \\ \hline
		2 & 0.4895 & 1.53 \\ \hline
		3 & 0.4135 & 1.77 \\ \hline
		4 & 0.3719 & 1.44 \\ \hline
		5 & 0.3385 & 1.31 \\ \hline
		6 & 0.3068 & 1.43 \\ \hline
		7 & 0.2777 & 1.38 \\ \hline
		8 & 0.2538 & 1.39 \\ \hline
		9 & 0.2361 & 1.78 \\ \hline
		10 & 0.2241 & 1.68 \\ \hline
		11 & 0.2163 & 1.66 \\ \hline
		12 & 0.2113 & 1.67 \\ \hline
		13 & 0.2081 & 1.62 \\ \hline
		14 & 0.2059 & 1.64 \\ \hline
		15 & 0.2043 & 1.62 \\ \hline
	\end{tabular}
	\caption{AltMin (m = 15, n = 20, rank = 10 p = 0.8, 0,9)}
	\label{exp12AM}
\end{table}

\begin{table}[ht]
	\centering
	\begin{tabular}{|c|c|c|c|}
	\hline Iteration & $e_k$ & MSE & Time \\ \hline
	0 & - & 4.32e-28 & 0.011  \\ \hline
	1 & 0.9179 & 0.2179 & 1.54 \\ \hline
	2 & 0.2830 & 0.2333 & 1.85 \\ \hline
	3 & 0.1122 & 0.2324 & 1.52 \\ \hline
	4 & 0.0482 & 0.2303 & 1.44 \\ \hline
	5 & 0.0210 & 0.2290 & 1.63 \\ \hline
	6 & 0.0097 & 0.2278 & 1.27 \\ \hline
	7 & 0.0043 & 0.2271 & 1.67 \\ \hline
	8 & 0.0019 & 0.2268 & 1.25 \\ \hline
	9 & 0.000867 & 0.2265 & 1.47 \\ \hline
	10 & 0.000401 & 0.2263 & 1.26 \\ \hline
	11 & 0.000198 & 0.2262 & 1.62 \\ \hline
	12 & 7.6425e-5 & 0.2262 & 1.43 \\ \hline
	13 & 3.8299e-5 & 0.2261 & 1.38 \\ \hline
	14 & 1.9896e-5 & 0.226070 & 1.72 \\ \hline
	15 & 1.8754e-5 & 0.226018 & 1.86 \\ \hline
	\end{tabular}
	\begin{tabular}{|c|c|c|c|}
		\hline Iteration & $e_k$ & MSE & Time \\ \hline
		0 & - & 4.36e-28 & 0.012 \\ \hline
		1 & 0.8986 & 0.2565 & 1.55 \\ \hline
		2 & 0.2864 & 0.2881 & 1.64 \\ \hline
		3 & 0.1230 & 0.2982 & 1.50 \\ \hline
		4 & 0.0555 & 0.3039 & 1.57 \\ \hline
		5 & 0.0245 & 0.3081 & 1.91 \\ \hline
		6 & 0.0112 & 0.3106 & 1.35 \\ \hline
		7 & 0.0048 & 0.3127 & 1.27 \\ \hline
		8 & 0.0023 & 0.3137 & 1.30 \\ \hline
		9 & 0.001043 & 0.3146 & 1.66 \\ \hline
		10 & 0.000542 & 0.3150 & 1.18 \\ \hline
		11 & 0.000202 & 0.3156 & 1.60 \\ \hline
		12 & 0.000191 & 0.3153 & 1.19 \\ \hline
		13 & 9.8068e-5 & 0.3155 & 1.62 \\ \hline
		14 & 2.4713e-5 & 0.3159 & 1.65 \\ \hline
		15 & 1.0938e-5 & 0.3160 & 1.66 \\ \hline
	\end{tabular}
	\caption{IRM (m = 15, n = 20, rank = 10 p = 0.8, 0,9)}
	\label{exp12IRM}
\end{table}

\begin{table}[ht]
	\centering
	\begin{tabular}{|c|c|c|}
	\hline Iteration & MSE & Time \\ \hline
	1 & 0.6345 & 68.09 \\ \hline
	2 & 0.2117 & 56.78 \\ \hline
	3 & 0.10748 & 48.32 \\ \hline
	4 & 0.0682 & 44.02 \\ \hline
	5 & 0.0482 & 39.21 \\ \hline
	6 & 0.0359 & 40.14 \\ \hline
	\end{tabular}
	\begin{tabular}{|c|c|c|c|}
		\hline Iteration & $e_k$ & MSE & Time \\ \hline
		0 & - & 7.92e-23 & 0.063 \\ \hline
		1 & 0.8517 & 0.3419 & 86.92 \\ \hline
		2 & 0.1547 & 0.3760 & 100.83 \\ \hline
		3 & 0.0371 & 0.3853 & 107.84 \\ \hline
		4 & 0.0101 & 0.3893 & 79.91 \\ \hline
		5 & 0.0033 & 0.3904 & 77.42 \\ \hline
		6 & 0.0022 & 0.3892 & 84.20 \\ \hline
	\end{tabular}
	\caption{AltMin, IRM (m = 30, n = 45, rank = 20, p = 0.8)}
	\label{exp3}
\end{table}

\begin{figure}
	\centering
    \includegraphics[scale=.30]{AM_15x20_08_MSE.jpg}
    \includegraphics[scale=.30]{AM_15x20_08_time.jpg}
    \includegraphics[scale=.30]{IRM_15x20_08_e.jpg}
   	\includegraphics[scale=.30]{IRM_15x20_08_time.jpg}
    \caption{\label{15x20} AltMin, IRM (m = 15, n = 20, rank = 10, p = 0.8)}
  \end{figure}
  
 \begin{figure}
 	\centering
	\includegraphics[scale=.39]{AM_30x45_08_MSE.jpg}
    \includegraphics[scale=.39]{AM_30x45_08_time.jpg}
    \includegraphics[scale=.39]{IRM_30x45_08_e.jpg}
   	\includegraphics[scale=.39]{IRM_30x45_08_time.jpg}
   \caption{\label{30x45} AltMin, IRM (m = 30, n = 45, rank = 20, p = 0.8)}
 \end{figure}

На рисунках \ref{15x20} и \ref{30x45} видно, что, как и предполагалось, метод IRM постепенно уменьшает $n-r$-е наименьшее собственное значение. Достаточно малым для достижения заданного ранга оно становится на 7 итерации для матрицы 15 на 20 и на 5 итерации для матрицы 30 на 45. Также из графиков можно сделать вывод, что алгоритм AltMin постепенно уменьшает средний квадрат ошибки, причем в третьем эксперименте время выполнения каждой последующей операции уменьшается, что уменьшает время работы в целом.

В первом эксперименте (левые части таблиц \ref{exp12AM} и \ref{exp12IRM}) при более разреженной матрице алгоритм AltMin имеет значительно лучшую точность, по сравнению с IRM. IRM же был немного быстрее, однако стоит учитывать, что этот метод должен продолжать работу, пока не будет достигнут заданный ранг.

При $p = 0.9$ (правые части таблиц \ref{exp12AM} и \ref{exp12IRM}) точность методов и время их работы схожи. 

Во третьем эксперименте (таблица \ref{exp3}) используется матрица большего размера. В нем лучше себя показал метод AltMin, дав более высокую точность и скорость выполнения.

В ходе сравнения заметно лучше показал себя метод попеременной минимизации. В целом он дал более точное решение и показал более быструю сходимость. К тому же метод AltMin с первой итерации дает решение заданного ранга, то есть его $n - r$ наименьших собственных значений гарантированно равны 0.

Сравним полученное время работы методов с теоретической оценкой. При рассмотрении метода IRM приводилась оценка времени выполнения очередной итерации $O(n^6)$. В случае прямоугольной матрицы мы дополнительно приводим квадратную матрицу размера $m+n$ к рангу $m$. Поэтому оценкой в нашем случае является величина $O((m + n)^6)$. Для метода AltMin оценка определяется временем работы функции minimize, в которой используется алгоритм BFGS (Алгоритм Бройдена — Флетчера — Гольдфарба — Шанно). Этот алгоритм имеет оценку $O(n^2)$. Так как на вход функции подается максимум $nr$ переменных, оценкой времени выполнения итерации нашего метода будет величина $O((nr)^2)$

Для сравнения возьмем 1 и 3 эксперименты (матрицы одинаково разреженны и разного размера) Вместо второго эксперимента проведем новый с тем же $p = 0.8$ и матрицей размера 20 на 30, которую будем приводить к рангу 15.

Для метода AltMin среднее время выполнения итерации в первом эксперименте составило 1.95, во втором -- 10.504, в третьем -- 49.43.

Для метода IRM -- 1.53, 7.808 и 89.52 соответственно (не считая 0-ю итерацию).

Постоим графики оценок и посмотрим насколько они совпадают с полученными значениями. В качестве коэффициента перед $(m + n)^6$ возьмем отношение среднего времени выполнения итерации алгоритма IRM к значению функции $f(x) = x^6$ в точке $m + n$ для 2 эксперимента. В качестве коэффициента перед $(nr)^2$ возьмем отношение среднего времени выполнения итерации алгоритма AltMin к значению функции $f(x) = x^2$ в точке $nr$ для 2 эксперимента.

\begin{figure}
	\centering
    \includegraphics[scale=.39]{Graph_AltMin.jpg}
    \includegraphics[scale=.39]{Graph_IRM.jpg}
    \caption{\label{est} AltMin, IRM}
  \end{figure}

Из графиков на рисунке \ref{est} видно, что полученные результаты совпадают с теоретической оценкой.

\chapter{Применение}

\section{Заполнение матрицы}
Заполнение матрицы - это задача предсказания элементов некоторой неизвестной искомой матрицы $Y \in R^{m \times n}$ на основе случайного подмножества наблюдаемых элементов $E \subset [m] \times [n]$. Например, в известной задаче Netflix $m$ представляет количество пользователей, $n$ представляет количество фильмов, а $Y_{i, j}$ - рейтинг, который пользователь $i$ дает фильму $j$. Один из подходов к заполнению матрицы $Y$ - найти матрицу $A$ низкого ранга, которая приблизительно совпадает с $Y$ на элементах $E$ (в смысле среднего квадрата ошибки).

Тогда можно сформулировать задачу:
\begin{equation}
	\begin{aligned}
		&\min_A \frac{1}{|E|} \sum_{(i,j) \in E} (A_{i, j} - Y_{i,j})^2 \\
		& s.t. \ rank(A) \le r
	\end{aligned}
\end{equation}

\section{Аппроксимация матрицей низкого ранга}
Очень распространенная задача при анализе данных -- найти матрицу $A$ низкого ранга, которая аппроксимирует данную матрицу $Y$, а именно решение 
\begin{equation}
	\begin{aligned}
		& \min_A d(A, Y) \\ 
		& s.t. \ rank (A) \le r
  	\end{aligned}
\end{equation}
где $d$ -- некоторая мера несоответствия. Для простоты предположим, что $Y \in R^{n \times n}$. Когда $d (A, V)$ -- нормированная норма Фробениуса $d (A, V) = \frac{1}{n^2} \sum_{i, j} (A_{i,j} - Y_{i,j})^2$, эта задача решается эффективно через SVD. Однако хорошо известно, что из-за использования нормы Фробениуса эта процедура чувствительна к выбросам.

\section{Задача идентификации систем}
Рассматривая линейную инвариантную во времени систему с входом $u \in \mathbb{R}^m$ и выходом $y \in \mathbb{R}^p$, в этой задаче требуется идентифицировать линейную систему через $T$ выборок $w_d := (u_d, y_d) \in (\mathbb{R}^m \times \mathbb{R}^p)^T$. Такая задача идентификации эквивалентна поиску полной матрицы рангов строк $R = [R_0, R_1, ..., R_l]$ такой, что
\begin{equation}
	\label{SIP}
	R_0w(t) + R_1w(t+1) + \ldots + R_lw(t+l)=0, \ \forall t = 1,2,\ldots, T-l
\end{equation}
где l - отставание идентифицированной модели. Выражение в \ref{SIP} можно переписать в виде $R \mathscr{H}_{l+1,T-l}(w) = 0$, где $\mathscr{H}_{l+1,T-l}(w) = \begin{bmatrix} w(1) && w(2) && \dots && w(T-l) \\ w(2) && w(3) && \dots && w(T-l + 1) \\ \vdots && \vdots && \dots && \vdots \\w(2) && w(3) && \dots && w(T-l + 1) \end{bmatrix}$

Поскольку размерность строки $R$ равна количеству выходов $p$ модели, ранг $\mathscr{H}_{l + 1, T - l} (w)$ ограничен $p$. В результате проблема идентификации эквивалентна следующей задаче аппроксимации при условии низкого ранга:
\begin{equation}
  \begin{aligned}
  	&\min_w \|w - w_d\|_F^2 \\ &s.t. \ rank(\mathscr{H}_{l+1, T-l}(w)) \le r
  \end{aligned}
\end{equation}
где $r = (m + p)(l + 1) - p$.

%\section{Structured H2 Controller Design}
%\section{Cardinality Minimization Problem}


\begin{thebibliography}{00}
	\bibitem{RCOP} Chuangchuang Sun, Ran Dai
	\emph{Rank-Constrained Optimization: Algorithms and Applications}, 2018

	\bibitem{GECO} Shai Shalev-Shwartz, Alon Gonen, Ohad Shamir
	 \emph{Large-Scale Convex Minimization with a Low-Rank Constraint
}, June 2011
	\bibitem{AMalgo} Prateek Jain, Praneeth Netrapalli, Sujay Sanghavi
        \emph{Low Rank Matrix Completion using Alternating Minimization}, November 2016

	\bibitem{trunc:norm} Prateek Jain, Praneeth Netrapalli, Sujay Sanghavi
        \emph{ Low-rank Matrix Completion using Alternating Minimization}, December 4, 2012

	\bibitem{math:ts} Moritz Hardt
        \emph{Understanding Alternating Minimization for Matrix Completion}, May 15, 2014

\end{thebibliography}
\addcontentsline{toc}{chapter}{Литература}
\end{document}