SVD-разложение и его практические приложения с исходным кодом на языке Matlab.
(SVD - decomposition and its practical applications with source code in the language Matlab.
\documentclass[a4paper,12pt]{article}
\usepackage[noTeX]{mmap} % T1 fonts + mathematics (Unicode)
\usepackage[utf8]{inputenc} % translates various standard and other input encodings into a ‘LATEX internal language’
\usepackage[english,russian]{babel} % english and russian language
\usepackage{amssymb,amsmath,amsthm,amsfonts} % math symbols
\usepackage[framed,numbered,autolinebreaks]{mcode} % lstlisting for matlab code
\usepackage{color} % text color
\usepackage{wrapfig, graphicx} % figures
\usepackage{listings} % for source code
\usepackage{indentfirst} % indent at the beginning of the paragraph
\usepackage{hyperref} % links
\usepackage{caption,subcaption} % subfigures
\usepackage{tikz} % package for creating graphics programmatically
\usetikzlibrary{matrix,positioning,decorations.pathreplacing,calc} % used libraries for tikz package
\usepackage[margin=2cm]{geometry} % document margins
\usepackage{multirow}
\usepackage{kbordermatrix}
\usepackage{stackengine}
\newcommand{\HRule}{\rule{\linewidth}{0.5mm}}
\newcommand{\quotes}[1]{``#1''}
\numberwithin{equation}{subsection}
\newcommand{\noindex}{\hspace*{-0.8em}}
\DeclareMathOperator*{\argmax}{\arg\!\max}
\DeclareMathOperator*{\sign}{sign}
\DeclareMathOperator*{\rank}{rank}
\DeclareMathOperator{\Span}{span}
\newcommand{\matlabel}[2]{% \matlabel{<label>}{<stuff>}
\begin{array}{@{}c@{}} \mbox{\small$#1$} \\ #2 \end{array}
}
\DeclareMathOperator{\Mcol}{col}
\DeclareMathOperator{\Mrow}{row}
\DeclareMathOperator{\Mnull}{null}
\begin{document}
\input{title.tex}
\newpage
\tableofcontents
\newpage % =============================================================================================================================================
\section{Введение} \label{intro}
\subsection{Краткое описание SVD-разложения} \label{intro:description}
Одной из самых плодотворных идей в теории матриц является матричное разложение или каноническая форма. В последнее время матричные разложения стали оплотом численных методов линейной алгебры, которые служат основой решения множества проблем.
Из многочисленных разложений матриц, сингулярному разложению, которое представляет из себя факторизацию матрицы $\mathbf{A}$ в произведение $\mathbf{U\Sigma V^H}$, где $\mathbf{U,V}$ -- унитарные матрицы, а $\mathbf{\Sigma}$ -- диагональная матрица, отводится особое место. На это есть несколько причин:
\begin{itemize}
\item Во-первых, тот факт, что в разложении участвуют унитарные матрицы делает его идеальным механизмом для геометризации преобразования $\mathbf{A}$ в $n$-мерном пространстве.
\item Во-вторых, сингулярное разложение является устойчивым, т.е. малым возмущениям матрицы $\mathbf{A}$ соответствуют малые возмущения матрицы $\mathbf{\Sigma}$ и наоборот.
\item В-третьих, диагональная матрица $\mathbf{\Sigma}$ позволяет легко понять является ли матрица A почти вырожденной и, если она таковой является, сингулярное разложение дает возможность понизить ранг матрицы $\mathbf{A}$ с наименьшей погрешностью.
\item Последним по порядку, но не по значению является тот факт, что благодаря усилиям математика Gene Golub были получены эффективные, устойчивые алгоритмы вычисления сингулярного разложения на ЭВМ, которые используются уже на протяжении более 40 лет и включены во все математические пакеты сегодняшнего времени.
\end{itemize}
\subsection{Историческая справка} \label{intro:history}
Интригующим является тот факт, что большинство классических матричных разложений были получены задолго до использования матриц и матричной нотации: они задавались с помощью определителей, линейных систем уравнений и, особенно, билинейных и квадратичных форм. Отцом развития данного направления исследований стал Gauss. Один из его алгоритмов, разработанный им в 1823 году, раскладывает матрицу квадратичной формы $\mathbf{x^TAx}$ в произведение $\mathbf{RD^{-1}R}$, где $\mathbf{D}$ -- диагональная матрица, а $\mathbf{R}$ -- верхнетреугольная матрица с элементами матрицы $\mathbf{D}$ на главной диагонали. Gauss также разработал алгоритм для эффективного нахождения обратной матрицы, который преобразовывал систему вида $\mathbf{y = Ax}$ к системе вида $\mathbf{x = By}$. Его умение манипулировать квадратичными формами и системами уравнений позволили ему разработать метод наименьших квадратов, который лег в основу методов машинного обучения около 30 лет назад.
Cauchy в 1829 году установил свойства собственных значений и собственных векторов симметричных матриц, рассматривая соответствующие им однородные системы уравнений. В 1846 году Jacobi опубликовал алгоритм диагонализации симметричных матриц, а в посмертной статье 1857 года он получил $LU$-разложение билинейных форм. Weierstrass в 1868 году получил канонические формы для пар билинейных функций -- то, что в наше время называется обобщенной проблемой собственных значений: $\mathbf{Ax = \mu Bx}$, где $\mathbf{B}$ -- некоторая матрица. Подводя итог, иожно сказать, что разработка сингулярного разложения -- результат, полученный в результате многолетних исследований билинейных форм.
Сингулярное разложение было первоначально разработано в дифференциальной геометрии при изучении свойств билинейных форм учеными Eugenio Beltrami и Camille Jordan независимо в 1873 и 1874 годах соответственно. James Joseph Sylvester пришел к понятию сингулярного разложения для квадратных матриц в 1889 году и, вероятно, также независимо от Eugenio Beltrami и Camille Jordan. Sylvester называл сингулярные значения каноническими множителями матрицы. Четвертый математик, который открыл сингулярное разложение независимо от других -- Autonne в 1915 году, который пришел к сингулярному разложению с помощью полярного разложения. Первое доказательство сингулярного разложения для прямоугольных и комплексных матриц было осуществлено математиками Carl Eckart и Gale Young в 1936 году.
В 1907 году Erhard Schmidt определил аналог сингулярного разложения для интегральных операторов, не зная о том, что параллельно ведется работа над сингулярным разложением для конечномерных матриц. Его теория была дальше развита математиком Émile Picard в 1910 году. Именно Picard первым назвал $\sigma_k$ сингулярными значениями.
Практические методы вычисления SVD-разложения восходят к работам Kogbetliantz в 1954, 1955 и Hestenes в 1958, алгоритм которых сильно напоминает метод Якоби для вычисления собственных значений с использованием поворотов Гивенса. Однако этот алгоритм был заменен методом, разработанным математиками Gene Golub и William Kahan в 1965, который использует преобразования Хаусхолдера. В 1970 году Golub и Christian Reinsch опубликовали вариант алгоритма Golub/Kahan, который до сих пор является одним из самых используемых.
\section{Теория} \label{theory}
% Доказательство SVD-разложения математика Beltrami
\subsection{Beltrami} \label{theory:beltrami}
Beltrami и Jordan считаются основателями теории сингулярного разложения. Beltrami -- за то, что он первым опубликовал работу о сингулярном разложении, а Jordan -- за элегантность и полноту своей работы. Работа Beltrami появилась в журнале \quotes{Journal of Mathematics for the Use of the Students of the Italian Universities} в 1873 году, основная цель которой заключалась в том, чтобы ознакомить студентов с билинейными формами.
\subsubsection{Вывод SVD-разложения}
Beltrami начинает свое доказательство с рассмотрения билинейной формы
\[
f(x,y) = x^TAy,
\]
где $A$ -- действительная квадратная матрица $n$-го порядка. Если сделать замену вида:
\[
x = U\xi \text{ и } y = V\eta,
\]
то билинейную форму можно переписать в виде:
\[
f(x,y) = \xi^TS\eta,
\]
где
\begin{equation} \label{beltrami-UtAV}
S = U^TAV.
\end{equation}
Если потребовать, чтобы матрицы $U,V$ были ортогональными, то существует всего $n^2-n$ степеней свободы при их задании, которые можно использовать, чтобы обнулить внедиагональные элементы матрицы $S$.
Докажем это. Будем выяснять количество параметров, задающих ортогональный базис или, что то же самое, элементы ортогональной матрицы. Рассмотрим 2-мерное пространство. Единичный вектор в двумерном пространстве имеет всего 1 степень свободы -- угол вращения вокруг начала координат, т.к. длина у него постоянна. Учитывая тот факт, что второй ортогональный вектор, перпендикулярный первому задается однозначно с точностью до знака, но знак не является степенью свободы, мы приходим к выводу, что репер в 2-мерном пространстве задается одной степенью свободы. Рассмотрим 3-мерное пространство. Единичный вектор в нем задается с помощью 2-ух степеней свободы (широта и долгота), а остальные 2 вектора лежат в плоскости, перпендикулярной этому вектору, т.е. в 2-мерном пространстве, которое уже было нами рассмотрено. Т.о. репер в 3-мерном пространстве задается с помощью $2+1=3$ степеней свободы. Продолжая размышления по индукции, мы приходим к выводу, что репер в $n$-мерном пространстве задается с помощью $1+2+...+(n-1) = \dfrac{n(n-1)}{2}$ степеней свободы. Полученный результат показывает, что одна ортогональная матрица задается с помощью $\dfrac{n(n-1)}{2}$ степеней свободы, а так как в нашем случае 2 такие матрицы, то мы имеем $n^2 - n$ степеней свободы при их задании.
Предположим, что $S$ -- диагональная матирца, т.е. $S = \Sigma = diag(\sigma_1,...,\sigma_n)$. Тогда из (\ref{beltrami-UtAV}) и ортогональности $V,U$ следует, что
\begin{eqnarray}
U^TA = \Sigma V^T \label{UtA=SVt} \\
AV = U\Sigma \label{AV=US}
\end{eqnarray}
Выражая из (\ref{AV=US}) отдельно матрицы $U$ и $V$ и, подставляя полученные выражения в (\ref{UtA=SVt}), получаем уравнения
\begin{eqnarray}
U^T(AA^T) = \Sigma^2U^T \label{UtAAt=S2Ut} \\
(A^TA)V = V\Sigma^2 \label{AtAV=VS2}
\end{eqnarray}
Покажем, что $\sigma_i$ -- корни уравнений
\begin{eqnarray}
\det(AA^T - \sigma^2I) = 0, \label{AAt=s2I} \\
\det(A^TA - \sigma^2I) = 0. \label{AtA=s2I}
\end{eqnarray}
Для этого покажем, что набор собственных значений инвариантен относительно преобразования подобия $B = U\Lambda U^{-1}$:
\begin{eqnarray}
\det(B-\mu I) = \det(U\Lambda U^{-1} - \mu I) = \det(U(\Lambda - \mu I)U^{-1}) = \nonumber \\
= \det(\Lambda - \mu I) \cdot \det(U) \cdot \det(U)^{-1} = \det(\Lambda - \mu I). \nonumber
\end{eqnarray}
Покажем, что (\ref{AAt=s2I}), (\ref{AtA=s2I}) являются эквивалентными алгебраическими уравнениями. Для этого необходимо доказать, что множество ${\sigma_i^2}$ является множеством собственных значений матриц $AA^T$ и $A^TA$.
Предположим, что $\lambda$ -- собственное значение матрицы $A^TA$, т.е.\[ A^TAx = \lambda x. \] Тогда, домножив обе части на $A$, мы получим следующее выражение: \[ AA^TAx = \lambda Ax. \] Если сделать замену переменных $y := Ax$, то мы докажем, что если $\lambda$ -- собственное значение матрицы $A^TA$, то оно также будет и собственным значением матрицы $AA^T$: \[ AA^Ty = \lambda y. \]
Таким образом, мы показали, что (\ref{AAt=s2I}), (\ref{AtA=s2I}) являются эквивалентными алгебраическими уравнениями, т.е. корни $\sigma_i$ будут одинаковыми в обоих случаях.
Необходимо заметить, что вывод, предложенный Beltrami предполагает, что матрица $\Sigma$, а следовательно и матрица $A$ -- невырожденные.
Т.к. матрицы $A^TA$ и $AA^T$ симметричные и положительно определенные, то они имеют полный набор действительных, положительных собственных значений.
На основании вышеизложенных размышлений Beltrami написал алгоритм нахождения сингулярного разложения:
\begin{enumerate}
\item Найти корни уравнения (\ref{AAt=s2I}).
\item Определить $U$ из (\ref{UtAAt=S2Ut}). Здесь Beltrami добавляет, что столбцы $U$ определены с точностью до фактора $\pm 1$. Этот факт и то, что матрица $U$ ортогональна, справедливы только в том случае, когда все $\sigma_i$ различные.
\item Определить $V$ из (\ref{UtA=SVt}). Этот шаг требует, чтобы $\Sigma$ была невырожденной.
\end{enumerate}
Таким образом Beltrami вывел SVD-разложение для действительных, квадратных, невырожденных матриц с различными сингулярными значениями.
% Доказательство SVD-разложения математика Sylvester
\subsection{Sylvester}\label{theory:sylvester}
Sylvester написал 2 статьи по сингулярному разложению. В первой статье он привел итеративный алгоритм редуцирования квадратичной формы к диагональному виду, а в примечании к статье указал, что аналогичный итерационный метод может диагонализовать билинейную форму. Во-второй статье он написал правило, по которому задача диагонализации билинейной формы сводилась к диагонализации квадратичной формы. Это правило оказалось ничем иным, как алгоритмом Beltrami.
\subsubsection{Правило нахождения сингулярного разложения}
Sylvester начинает с билинейной формы
\[
B(x,y) = x^TAy
\]
и рассматривает квадратичную форму
\[
M(x) = \sum_i \bigg( \dfrac{dB}{dy_i} \bigg)^2,
\]
которая является квадратичной формой вида
\[
M(x) = x^TAA^Tx.
\]
Пусть $M = \sum \lambda_i \xi_i^2$ -- каноническая форма матрицы $M$. Если $B$ имеет каноническую форму $B = \sum \sigma_i \xi_i \eta_i$, тогда $\sum (\sigma_i \xi)^2$ подобно $M = \sum \lambda_i \xi_i^2$, из чего следует, что $\lambda_i = \sigma_i^2$.
Чтобы найти змену, Sylvester вводит матрицы $M = AA^T$ и $N = A^TA$ и утверждает, что замена $x$ -- такая замена, которая диагонализирует матрицу $M$, а замена для $y$ -- такая замена, которая диагонализирует матрицу $N$. Это верно только в том случае, когда все сингулярные значения матрицы $A$ различны.
% Доказательство SVD-разложения математика Autonne
\subsection{Autonne} \label{theory:autonne}
\subsubsection{Полярное разложение}
Теорема: Любую матрицу A можно представить в виде произвдения ортогональной матрицы Q и симметричной матрицы S:
\begin{equation} \label{qs-decomposition}
A = QS,
\end{equation}
причем все собственные векторы симметричной матрицы S -- неотрицательные числа.
Геометрический смысл: любое линейное преобразование может быть представлено в виде двух последовательных преобразований:
\begin{enumerate}
\item масштабирование вдоль ортогональных направлений на неотрицательные коэффициенты
\item последовательность поворотов с возможными отражениями между ними
\end{enumerate}
Для того, чтобы доказать эту теорему, приведем пример алгоритма для построения данного разложения и покажем, что матрицы Q и S уникальны для каждой матрицы A.
\begin{equation} \label{qs-find-s}
\begin{aligned}
&A^TA = SQ^TQS = S^2 \\
&S = \sqrt{A^TA}
\end{aligned}
\end{equation}
Зная S и A можно найти Q:
\begin{equation} \label{qs-find-q}
Q = AS^{-1}
\end{equation}
Матрицы Q и S уникальны для каждой матрицы A. Покажем, что матрица Q является ортогональной:
\begin{equation} \nonumber
\begin{aligned}
Q^TQ &= S^{-1}A^TAS^{-1} = S^{-1}S^2S^{-1} = I \\
QQ^T &= AS^{-1}S{-1}A^T = A(S^2)^{-1}A^T = AA^{-1}(A^T)^{-1}A^T = I
\end{aligned}
\end{equation}
\subsubsection{Вывод SVD-разложения}
Воспользуемся полярным разложением и представим матрицу $A$ в виде:
\[
A = QS
\]
Воспользуемся свойством спектрального разложения симметричных матриц:
\[
S = X \Lambda X^T
\]
Симметричные матрицы имеют полный набор собственных значений, а собственные векторы могут быть выбраны ортогональными. Таким образом, приходим к следующей формуле:
\[
A = QX \Lambda X^T
\]
В данной формуле Q,X -- ортогональные матрицы, поэтому их произведением является также ортогональная матрица $Y = QX$.
Подставляя все в одну формулу, получаем SVD-разложение:
\begin{equation} \label{svd-ylx}
A = Y \Lambda X^T
\end{equation}
Или более традиционная запись:
\begin{equation} \label{svd}
A = U \Sigma V^T
\end{equation}
Коэффициенты $\sigma_1,...,\sigma_n$, стоящие на главной диагонали матрицы $\Sigma$, называются сингулярными значениями и эквивалентны собственным значениям матрицы $A^TA$.
\subsection{Сингулярные векторы} \label{theory:sing_vectors}
SVD-разложение -- аналог спектрального разложения, который работает с любыми матрицами.
Столбцы матрицы $V$ называются правыми сингулярными векторами и всегда ортогональны друг другу. Столбцы матрицы U называются левыми сингулярными векторами и также ортогональны друг другу.
Чтобы лучше понять, что такое SVD-разложение, представим матрицу $А$ размеров $n \times d$ как n точек $\{a^{(i)}\}_{i=1}^n$ в d-мерном пространстве и поставим задачу нахождения наилучшего k-размерного подпространства относительно этих $n$ точек.
\begin{equation} \nonumber
A_{n \times d} =
\begin{pmatrix}
a_{1}^{(1)} & a_{2}^{(1)} & \cdots & a_{d}^{(1)} \\
a_{1}^{(2)} & a_{2}^{(2)} & \cdots & a_{d}^{(2)} \\
\vdots & \vdots & \ddots & \vdots \\
a_{1}^{(n)} & a_{2}^{(n)} & \cdots & a_{d}^{(n)}
\end{pmatrix}
\end{equation}
В данном контексте \quotes{наилучшим} называется пространство с минимальным расстоянием до этих точек. Альтернативной задачей является задача наилучшего приближения, которая обычно решается с помощью МНК. В задаче наилучшего приближения минимизируется вертикальная ошибка, в то время как при нахождении наилучшего подпространства минимизируется перпендикулярное расстояние до точек.
\begin{figure}[!h]
\caption{Минимизируемые ошибки}
\centering
\includegraphics[width=0.95\textwidth]{pcalsa.png}
\end{figure}
Рассмотрим проекцию точки $x^{(i)}$ на подпространство, проходящее через начало координат:
\[
\sum_{k=1}^{d} (x^{(i)}_k)^2 = (l^{(i)})^2 + (h^{(i)})^2,
\]
где $l^{(i)}$ -- проекция точки $x^{(i)}$ на подпространство, а $h^{(i)}$ -- перпендикуляр от точки до подпространства.
\begin{align*}
\begin{tikzpicture}[scale=4]
\draw[<->] (0,1.5) -- (0,0) -- (1.5,0);
\draw[red,thick] (-0.125,-0.0625) -- (1.5,0.75);
\draw[blue,thick,->] (0,0) -- (1,0.5) node[anchor=base,black] at (0.5,0.28) {$l^{(i)}$};
\draw[fill] (0.65,1.2) circle [radius=0.02];
\draw[blue,thick,->] (1,0.5) -- (0.662,1.183) node[anchor=base,black] at (0.93,0.85) {$h^{(i)}$};
\draw[blue,thick,->] (0,0) -- (0.638,1.183) node[anchor=base,black] at (0.3,0.7) {$x^{(i)}$};
\end{tikzpicture}
\end{align*}
\[
(h^{(i)})^2 = \sum_{k=1}^{d} (x^{(i)}_k)^2 - (l^{(i)})^2
\]
Чтобы минимизировать сумму расстояний от точек до прямой можно минимизировать $\sum_{i=1}^n \big[\sum_{k=1}^{d} (x^{(i)}_k)^2 - (l^{(i)})^2 \big]$, однако т.к. $\sum_{i=1}^n \sum_{k=1}^{d} (x^{(i)}_k)^2$ -- константа, то задача минимизации суммы расстояний от точек до прямой эквивалентна задаче максимизации суммы проекций на прямую $\sum_{i=1}^n (l^{(i)})^2$.
Определим теперь сингулярные векторы матрицы $A_{n \times d}$. Представим прямую, проходящую через начало координат. Предположим, что $\mathbf{v}$ -- единичный вектор вдоль этой прямой. Длину проекции $\mathbf{a^{(i)}}$ ($i$-ой строки матрицы $A$) на $\mathbf{v}$ можно представить в виде $|\mathbf{a^{(i)} \cdot v}|$. Отсюда видно, что сумму квадратов проекции можно представить в виде $||A\mathbf{v}||^2$. Следовательно, наилучшая прямая будет та, которая максимизирует $||A\mathbf{v}||^2$ (по т. Пифагора) и, соответственно, минимизирует сумму квадратов расстояний от точек до прямой.
На основе предыдущих размышлений определим первый собственный вектор $\mathbf{v_1}$ матрицы $A$, который является вектором-столбцом, как наилучшую прямую, заданную этим вектором, для n точек в d-мерном пространстве.
\[
\mathbf{v_1} = \argmax_{||\mathbf{v}|| = 1}||A\mathbf{v}||
\]
Значение $\sigma_1(A) = ||A\mathbf{v_1}||$ называется первым сингулярным значением матрицы A. Отметим также, что $\sigma_1^2$ -- сумма квадратов проекций точек на прямую, заданную вектором $\mathbf{v_1}$.
\quotes{Жадный} подход к нахождению 2-мерного подпространства для матрицы $A$ берет $\mathbf{v_1}$ в качестве первого базисного вектора и находит наилучшее 2-мерное подпространство, содержащее $\mathbf{v_1}$. Для любого двумерного подпространства сумма квадратов проекций точек эквивалентна сумме квадратов проекций на $\mathbf{v_1}$ плюс сумма квадратов проекций на $\mathbf{v_2}$. Нет заранее никакой гарантии, что \quotes{жадный} подход даст наилучшее решение, однако в данном случае это так.
Второй сингулярный вектор находится следующим образом:
\[
\mathbf{v_2} = \argmax_{\mathbf{v} \bot \mathbf{v_1}, ||\mathbf{v}|| = 1} ||A\mathbf{v}||
\]
Значение $\sigma_2(A) = ||A\mathbf{v_2}||$ называется вторым сингулярным значением матрицы $A$.
Аналогично находятся и все оставшиеся векторы $d$-мерного пространства:
\[
\mathbf{v_{k+1}} = \argmax_{\mathbf{v} \bot \mathbf{v_1},...,\mathbf{v_k}, ||\mathbf{v}|| = 1} ||A\mathbf{v}||
\]
Докажем теорему о том, что жадный алгоритм в действительности находит наилучшее подпространство любой размерности.
Пусть $A$ -- матрица размерности $n\times d$, где $\mathbf{v_1,...,v_r}$ -- сингулярные векторы, которые заданы выше. Для $1 \leq k \leq r$, пусть $V_k$ -- подпространство, порожденное базисными векторами $\mathbf{v_1,...,v_r}$. Тогда для любого $k, V_k$ является наилучшим $k$-размерным подпространством $A$.
Будем проводить доказательство по индукции. Для $k=1$ утверждение, очевидно, выполняется. Предположим, что $V_{k-1}$ -- наилучшее $(k-1)$-мерное подпространство. Предположим, что $W_k$ -- наилучшее $k$-мерное подпространство. Выберем базис ${\mathbf{w_1,...,w_k}}$ таким образом, чтобы $\mathbf{w_k}$ было перпендикулярно $V_{k-1}$. Тогда
\[
\sum_{i=1}^k ||A \mathbf{w_i}||^2 \leq \sum_{i=1}^{k-1} ||A \mathbf{v_i}||^2 + ||A \mathbf{w_k}||^2,
\]
т.к. $V_{k-1}$ -- оптимальное $(k-1)$-мерное подпространство. Т.к. $\mathbf{w_k}$ перпендикулярно $\mathbf{v_1,...,v_{k-1}}$, то по определению $\mathbf{v_k}$: $||A\mathbf{w_k}||^2 \leq ||A\mathbf{v_k}||^2$. Таким образом
\[
\sum_{i=1}^k ||A \mathbf{w_i}||^2 \leq \sum_{i=1}^k ||A \mathbf{v_i}||^2,
\]
что показывает, что $V_k$ по крайней мере не хуже $W$, а следовательно -- оптимально.
Заметим, что $n$-вектор $A\mathbf{v_i}$ является списком длин (со знаком) проекций строк матрицы $A$ на вектор $\mathbf{v_i}$. Будем смотреть на $||A\mathbf{v_i}|| = \sigma_i(A)$ как на некую \quotes{компоненту} матрицы $A$ вдоль вектора $\mathbf{v_i}$. Для того, чтобы такая интерпретация имела смысл, необходимо чтобы сумма квадратов \quotes{компонент} матрицы $A$ вдоль каждой $\mathbf{v_i}$ давала бы квадрат всего \quotes{содержимого} матрицы $A$. Это на самом деле выполняется и является матричной аналогией разложения вектора по компонентам ортогонального базиса.
Рассмотрим строку $\mathbf{a_j}$ матрицы $A$. Т.к. $\mathbf{v_1,...,v_r}$ образует линейную оболочку всех строк матрицы $A$, то $\forall \mathbf{v} \perp \mathbf{v_1,...,v_r}:$ $\mathbf{a_j \cdot v} = 0$. Тогда для каждого $\mathbf{a_j}$ выполняется равенство $\sum_{i=1}^r (\mathbf{a_j \cdot v_i})^2 = ||\mathbf{a_j}||^2$. Суммируя по всем строкам получаем:
\[
\sum_{j=1}^n ||\mathbf{a_j}||^2 = \sum_{j=1}^n \sum_{i=1}^r (\mathbf{a_j \cdot v_i})^2 = \sum_{i=1}^r \sum_{j=1}^n (\mathbf{a_j \cdot v_i})^2 = \sum_{i=1}^r ||A \mathbf{v_i}||^2 = \sum_{i=1}^r \sigma_i^2(A).
\]
При этом следует учесть, что $\sum_{j=1}^n ||\mathbf{a_j}||^2 = \sum_{j=1}^n \sum_{k=1}^d a_{jk}^2$ -- сумма квадратов всех элементов матрицы $A$. Таким образом, сумма квадратов сингулярных значений матрицы $A$ на самом деле эквивалентна квадрату всего \quotes{содержимого} матрицы $A$. Существует специальная норма, связанная с этой величиной -- норма Фробениуса, обозначаемая $||A||_F$ и определяющаяся
\[
||A||_F = \sqrt{\sum_{j,k}a_{jk}^2}
\]
Следствием вышеизложенных выкладок служит утверждение о том, что сумма квадратов сингулярных значений матрицы $A$ равна квадрату нормы Фробениуса.
\[
\sum_i \sigma_i^2(A) = ||A||_F^2
\]
Матрица $A$ может быть полностью описана тем, как она трансформирует векторы $\mathbf{v_i}$. Любой вектор $\mathbf{v}$ может быть представлен в виде линейной комбинации векторов $\mathbf{v_1,...,v_r}$ и вектора, перпендикулярного всем $\mathbf{v_i}$. $A\mathbf{v}$ -- такая же линейная комбинация $A\mathbf{v_1},...,A\mathbf{v_r}$, как и $\mathbf{v}$ -- линейная комбинация векторов $\mathbf{v_1},...,\mathbf{v_r}$. Таким образом, $A\mathbf{v_1},...,A\mathbf{v_r}$ образуют фундаментальное множество векторов, ассоциированных с матрицей $A$. После нормализации получаем:
\[
\mathbf{u_i} = \dfrac{1}{\sigma_i(A)} A\mathbf{v_i}.
\]
Векторы $\mathbf{u_1,...,u_r}$ называются левыми сингулярными векторами матрицы $A$, векторы $\mathbf{v_1,...,v_r}$ -- правыми сингулярными векторами.
Покажем, что левые сингулярные векторы также являются ортогональными.\\
Для доказательства данной теоремы забежим несколько вперед и воспользуемся теоремой из (\ref{theory:sing_decomp}), которая утверждает, что
\[ A = \sum_{i=1}^r \sigma_i\mathbf{u_iv_i^T}. \]
Составим набор матриц $\{B_{k}\}$ таких, что
\begin{equation} \label{eq:B_k}
B_{k} = A - \sum_{i=1}^k \sigma_i\mathbf{u_iv_i^T} = \sum_{i=k+1}^r \sigma_i\mathbf{u_iv_i^T}.
\end{equation}
Из вида уравнения (\ref{eq:B_k}) видно, что правыми сингулярными векторами матрицы $B_k$ являются векторы $\{\mathbf{v_{k+1}},...,\mathbf{v_r}\}$, которые в свою очередь являются также подмножеством множества правых сингулярных векторов матрицы $A$.
Таким образом,
\[ B_{r-1} = \sigma_r \mathbf{u_rv_r^T} \]
имеет единственный левый сингулярный вектор $\mathbf{u_r}$. В данном случае система левых сингулярных векторов $\{\mathbf{u_r}\}$ является ортонормированной системой. По индукции предположим, что матрица $B_{i}$ имеет ортонормированную систему левых сингулярных векторов $\{\mathbf{u_{i+1}},...,\mathbf{u_r}\}$. Необходимо показать, что матрица $B_{i-1}$ также имеет ортонормированную систему левых сингулярных векторов. Для этого необходимо показать, что вектор $\mathbf{u_i}$ ортогонален множеству векторов $\{\mathbf{u_{i+1}},...,\mathbf{u_r}\}$.\\
Предположим, что это не так, что $\exists\ j > i$ такие, для которых $\mathbf{u_i^T}\mathbf{u_j} \neq 0$. Без ограничения общности предположим, что $\mathbf{u_i^T}\mathbf{u_j} > 0$. В таком случае для бесконечно малого $\varepsilon > 0$ \[ A\bigg( \dfrac{\mathbf{v_i}+\varepsilon\mathbf{v_j}}{||\mathbf{v_i}+\varepsilon\mathbf{v_j||_2}} \bigg) = \dfrac{\sigma_i\mathbf{u_i} + \varepsilon\sigma_j\mathbf{u_j}}{||\mathbf{v_i}+\varepsilon\mathbf{v_j}||_2} = \dfrac{\sigma_i\mathbf{u_i} + \varepsilon\sigma_j\mathbf{u_j}}{\sqrt{||\mathbf{v_i}||_2^2+||\varepsilon\mathbf{v_j}||_2^2}} = \dfrac{\sigma_i\mathbf{u_i} + \varepsilon\sigma_j\mathbf{u_j}}{\sqrt{1+\varepsilon^2}}. \]
Из построения правых сингулярных векторов известно, что \[ \mathbf{v_i} = \argmax_{\mathbf{v} \bot \mathbf{v_1},...,\mathbf{v_{i-1}}, ||\mathbf{v}|| = 1} ||A\mathbf{v}||_2\ \Rightarrow\ \forall \mathbf{x} \in \Span\{\mathbf{v_i},...,\mathbf{v_r}\},\ \mathbf{x} \neq \mathbf{v_i},\ ||A\mathbf{x}||_2 < ||A\mathbf{v_i}||_2 = \sigma_i. \]
Вспомним, что $\sigma_i \mathbf{u_i} = A\mathbf{v_i}$, тогда из вышеприведенных рассуждений получим неравенство \[ \forall \mathbf{x} \in \Span\{\mathbf{v_i},...,\mathbf{v_r}\},\ \mathbf{x} \neq \mathbf{v_i}:\ ||A\mathbf{x}||_2 < ||\sigma_i \mathbf{u_i}||_2 \Leftrightarrow \mathbf{u_i^T}A\mathbf{x} < \sigma_i \mathbf{u_i^T} \mathbf{u_i} \Leftrightarrow \mathbf{u_i^T}A\mathbf{x} < \sigma_i \]
\begin{equation} \nonumber
\begin{aligned}
\mathbf{u_i^T}A\bigg( \dfrac{\mathbf{v_i}+\varepsilon\mathbf{v_j}}{||\mathbf{v_i}+\varepsilon\mathbf{v_j||_2}} \bigg) = \mathbf{u_i^T} \dfrac{\sigma_i\mathbf{u_i} + \varepsilon\sigma_j\mathbf{u_j}}{\sqrt{1+\varepsilon^2}} &= (\sigma_i\ + \varepsilon\sigma_j\mathbf{u_i^T}\mathbf{u_j})\big( 1 - \frac{\varepsilon^2}{2} + O(\varepsilon^4) \big)\\
&= \sigma_i\ + \varepsilon\sigma_j\mathbf{u_i^T}\mathbf{u_j} + O(\varepsilon^2) > \sigma_i\
\end{aligned}
\end{equation}
Полученное противоречие доказывает, что система левых сингулярных векторов $\{\mathbf{u_i},...,\mathbf{u_r}\}$ ортонормирована. Продолжая по индукции, получим, что все левые сингулярные векторы матрицы $A$ образуют ортонормированную систему векторов.
\subsection{Сингулярное разложение} \label{theory:sing_decomp}
Пусть $A$ -- матрица размерности $n \times d$, $\mathbf{v_1,...,v_r}$ -- ее правые сингулярные векторы, которым соответствуют собственные значения $\sigma_1,...,\sigma_r$. Тогда $\mathbf{u_i} = \dfrac{1}{\sigma_i}A\mathbf{v_i}, i = 1,...,r$ -- левые сингулярные веторы. Покажем, что
\[
A = \sum_{i=1}^r \sigma_i\mathbf{u_iv_i^T}.
\]
Для доказательства этой теоремы, докажем сперва следующую лемму: матрицы $A$ и $B$ идентичны тогда и только тогда, когда $\forall \mathbf{v}: A\mathbf{v} = B\mathbf{v}$.
Необходимость: ясно, что если $A=B$, то $\forall \mathbf{v}: A\mathbf{v} = B\mathbf{v}$.
Достаточность: предположим, что $\forall \mathbf{v}: A\mathbf{v} = B\mathbf{v}$. Пусть $\mathbf{e_i}$ -- вектор, все элементы которого нули, за исключением $i$-ой компоненты, которая равна 1. В таком случае $A\mathbf{e_i}$ -- $i$-ый столбец матрицы $A$ и, следовательно, $A=B$ если $\forall i: A\mathbf{e_i} = B\mathbf{e_i}$.
Воспользуемся только что доказанной леммой для того, чтобы доказать равенство
\[
A = \sum_{i=1}^r \sigma_i\mathbf{u_iv_i^T}.
\]
Для любого правого сингулярного вектора $\mathbf{v_j}$ выполняется равенство:
\[
A\mathbf{v_j} = \sum_{i=1}^r \sigma_i \mathbf{u_i} \mathbf{v_i^T} \mathbf{v_j},
\]
т.к. $\forall i \neq j: \mathbf{v_i^T} \mathbf{v_j} = 0$, т.к. сингулярные векторы ортогональны. Учитывая тот факт, что любой вектор $\mathbf{v}$ может быть представлен в виде линейной комбинации сингулярных векторов и вектора, перпендикулярного им, то $A\mathbf{v} = \sum_{i=1}^r \sigma_i \mathbf{u_i} \mathbf{v_i^T} \mathbf{v}$. Используя только что доказанную лемму получаем, что
\[
A = \sum_{i=1}^r \sigma_i\mathbf{u_iv_i^T}.
\]
\subsection{Геометрическая интерпретация} \label{theory:interpret_transform}
Геометрическая интерпретация сингулярного разложения заключается в следующем факте из геометрии: \textit{Образом любого линейного преобразования, заданного с помощью матрицы $m \times n$, примененного к единичной сфере является гиперэллипсоид.}
Гиперэллипсоид в $\mathbb{R}^m$ можно определить как поверхность, полученную в результате растяжения единичной сферы в $\mathbb{R}^m$ в $\sigma_1,...,\sigma_m$, возможно, в ноль раз вдоль каких-то ортогональных направлений $u_1,...,u_m \in \mathbb{R}^m$. Для определенности предположим, что $u_i$ -- единичные векторы. Тогда векторы $\{\sigma_i u_i\}$ -- главные полуоси гиперэллипсоида, длины которых $\sigma_1,...,\sigma_m$. Если $\rank{(A)} = r$, то ровно $r$ длин $\sigma_i$ ненулевые, и, в частности, если $m \ge n$, то не более, чем $n$ из них ненулевые.
Под единичной сферой подразумевается обычная евклидова сфера в $n$-мерном пространстве, т.е. единичная сфера для 2-нормы; обозначим ее как $S$. Тогда $AS$, образ $S$ под действием преобразования $A$ -- гиперэллипсоид.
\begin{figure}[!h]
\centering
\includegraphics[width=0.85\textwidth]{img/svd_geom_interp.png}
\caption{Преобразование из единичной окружности в эллипс}
\label{fig:interpret_transform}
\end{figure}
Первым делом определим $n$ сингулярных значений преобразования $A$. Ими являются длины $n$ главных полуосей эллипса $AS$, обозначенных $\sigma_1,\sigma_2$. Обычно принято нумеровать сингулярные значения в убывающем порядке, $\sigma_1 \ge \sigma_2 \ge ... \ge \sigma_n > 0$.
Далее определим $n$ левых сингулярных векторов преобразования $A$. Ими являются единичные векторы $\{u_1,u_2\} \in AS$, направленные вдоль главных полуосей эллипса $AS$, пронумерованные таким образом, чтобы соответствовать номерам сингулярных значений. Поэтому $\sigma_iu_i$ -- $i$-ая по величине полуось.
Последними определим $n$ правых сингулярных векторов преобразования $A$. Ими являются единичные векторы $\{v_1,v_2\} \in S$, которые являются прообразами главных осей эллипса $AS$, пронумерованные таким образом, чтобы $Av_j = \sigma_ju_j$.
\subsection{Существование и единственность} \label{theory:uniqueness}
В (\ref{theory:beltrami}), (\ref{theory:sylvester}), (\ref{theory:autonne}) и (\ref{theory:sing_decomp}) было показано, как можно найти сингулярное разложение. Следующим необходимо доказать, что любая матрица может быть разложена таким образом.
Доказательство будем производить по индукции по размерности матрицы $A_{m \times n}$, где $A$ -- комплексная матрица размерности $m \times n$. Все нижеизложенные выкладки справедливы для $m \ge n$ (если $m < n$, то необходимо сперва транспонировать матрицу).
Для $n=1$ (и любой $m$) матрица $A$ является ветором-столбцом. Возьмем $V=1,\ \hat{\Sigma} = ||A||_2,\ \hat{U} = \dfrac{A}{||A||_2}$. Тогда очевидно, что мы нашли \quotes{сжатое} сингулярное разложение, т.е. \[ A = \hat{U} \hat{\Sigma} V^H. \] \quotes{Полное} сингулярное разложение можно получить, расширив матрицу $\hat{U}$ до $U$ с помощью ортогонализации Грамма-Шмидта и добавив несколько нулей к $\hat{\Sigma}$.
Предположим теперь, что SVD-разложение существует для матрицы размерности $(m-1) \times (n-1)$ и покажем, что оно также существует для матрицы размерности $m \times n$. Возьмем $\mathbf{v}_1 \in \mathbb{C}^n$ таким образом, чтобы \begin{equation} \nonumber
\left\{
\begin{aligned}
& ||\mathbf{v}_1||_2 = 1, \\
& ||A||_2 = \sup_{\mathbf{v}_1 \in \mathbb{C}^n} ||A\mathbf{v}_1||_2 > 0.
\end{aligned}
\right.
\end{equation}
Теперь возьмем \begin{equation} \label{theory:uniqueness:u_1}
\mathbf{u_1} = \dfrac{A\mathbf{v}_1}{||A\mathbf{v}_1||_2}.
\end{equation}
Далее воспользуемся ортогонализацией Грамма-Шмидта для того, чтобы расширить $\mathbf{u}_1$ и $\mathbf{v}_1$ до унитарных матриц, добавляя столбцы $\widetilde{U}_1$ и $\widetilde{V}_1$ соответственно. Т.о. матрицы $U_1,\ V_1$ выглядят следующим образом: \[ U_1 = \left[ \begin{array}{cc} \mathbf{u}_1 & \widetilde{U}_1 \end{array} \right], \quad V_1 = \left[ \begin{array}{cc} \mathbf{v}_1 & \widetilde{V}_1 \end{array} \right]. \] В таком случае получаем следующее выражение: \[ U_1^HAV_1 = \left[ \begin{array}{c} \mathbf{u}_1^H \\ \widetilde{U}_1 \end{array} \right] A \left[ \begin{array}{cc} \mathbf{v}_1 & \widetilde{V}_1 \end{array} \right] = \left[ \begin{array}{cc} \mathbf{u}_1^HA\mathbf{v}_1 & \mathbf{u}_1^HA\widetilde{V}_1 \\ \widetilde{U}_1^HA\mathbf{v}_1 & \widetilde{U}_1^HA\widetilde{V}_1 \end{array} \right]. \]
Используя (\ref{theory:uniqueness:u_1}) и \quotes{специальный} выбор $\mathbf{v}_1$ имеем: \[ \mathbf{u}_1^HA\mathbf{v}_1 = \dfrac{(A\mathbf{v}_1)^H}{||A\mathbf{v}_1||_2}A\mathbf{v}_1 = \dfrac{||A\mathbf{v}_1||_2^2}{||A\mathbf{v}_1||_2} = ||A\mathbf{v}_1||_2 = ||A||_2. \] Для этой величины вводится обозначение $\sigma_1 := ||A||_2$.
Опять, используя (\ref{theory:uniqueness:u_1}), получаем \[ \widetilde{U}_1^HA\mathbf{v}_1 = \widetilde{U}_1^H\mathbf{u}_1||A\mathbf{v}_1||_2 = \mathbf{0}, \] т.к. $\widetilde{U}_1^H$ построено с помощью ортогонализации Грамма-Шмидта таким образом, чтобы $\widetilde{U}_1^H\mathbf{u}_1 = \mathbf{0}$.
Теперь покажем, что $\omega^H := \mathbf{u}_1^HA\widetilde{V}_1 = \mathbf{0}^T$. Воспользуемся тем фактом, что $||A||_2 = \sigma_1$. Используя только что доказанные факты, запишем матрицу $\Sigma$ следующим образом: \[ \left\{ \begin{aligned}
& \Sigma = U_1^HAV_1 = \left[ \begin{array}{cc}
\sigma_1 & \omega^H \\
\mathbf{0} & \widetilde{A}
\end{array}
\right] \\
& A = U_1 \Sigma V_1^T,\ ||A||_2 = \sigma_1\ \Rightarrow\ \Sigma = U_1^HAV_1,\ ||\Sigma||_2 = \sigma_1.
\end{aligned} \right.
\]
Воспользуемся определением нормы матрицы $||A||_2 = \sup_{x} ||Ax||_2$. Возьмем специальный вектор $x = \left[ \begin{array}{c} \sigma_1 \\ \omega \end{array} \right]$, тогда
\[ \left\lVert \left[ \begin{array}{cc}
\sigma_1 & \omega^H \\
\mathbf{0} & \widetilde{A}
\end{array}
\right] \left[ \begin{array}{c}
\sigma_1 \\ \omega
\end{array}
\right] \right\rVert_2 \ge \sigma_1^2 + \omega^H\omega = \sqrt{\sigma_1^2 + \omega^H\omega} \left\lVert \left[ \begin{array}{c} \sigma_1 \\ \omega \end{array}
\right] \right\rVert_2, \] из чего следует, что $\sigma_1 \ge \sqrt{\sigma_1^2 + \omega^H\omega}$, откуда получается, что $\omega^H = \mathbf{0}^T$.
Т.о. мы показали, что \[ U_1^HAV_1 = \left[ \begin{array}{cc}
\sigma_1 & \mathbf{0}^T \\
\mathbf{0} & \widetilde{A}
\end{array}
\right], \] т.е. свели задачу с размерности $m \times n$ к $(m-1) \times (n-1)$, которая по условию индукции уже решена: $\widetilde{A} = U_2\Sigma_2V_2^H$. Тогда решение исходной задачи можно записать следующим образом: \[ U_1^HAV_1 = \left[ \begin{array}{cc}
\sigma_1 & \mathbf{0}^T \\
\mathbf{0} & U_2\Sigma_2V_2^H
\end{array}
\right] = \left[ \begin{array}{cc}
1 & \mathbf{0}^T \\
\mathbf{0} & U_2
\end{array}
\right]\left[ \begin{array}{cc}
\sigma_1 & \mathbf{0}^T \\
\mathbf{0} & \Sigma_2
\end{array}
\right]\left[ \begin{array}{cc}
1 & \mathbf{0}^T \\
\mathbf{0} & V_2
\end{array}
\right]^H, \] или, если выразить матрицу $A$: \[ A = U_1\left[ \begin{array}{cc}
1 & \mathbf{0}^T \\
\mathbf{0} & U_2
\end{array}
\right]\left[ \begin{array}{cc}
\sigma_1 & \mathbf{0}^T \\
\mathbf{0} & \Sigma_2
\end{array}
\right]\left[ \begin{array}{cc}
1 & \mathbf{0}^T \\
\mathbf{0} & V_2
\end{array}
\right]^H V_1^H, \] что также является SVD-разложением, т.к. произведение унитраных матриц -- унитарная матрица.
Сингулярное разложение в общем случае не является единственным. Однако сингулярные значения $\{\sigma_j\}$ уникальны для каждой матрицы $A$. Если все сингулярные значения различные, то левые и правые сингулярные векторы определены однозначно с точностью до комплексных знаков. Если все длины главных полуосей гиперэллипса различны, то эти оси однозначно определяют сингулярные векторы и значения с геометрической точки зрения.
\newpage % =======================================================================================================================================
\section{Математический метод}
\subsection{Описание алгоритма}
\subsubsection{Для матриц $n\times n$, $r = n$}
Рассмотрим невырожденную квадратную матрицу $A_{n \times n}$, $r = \rank(A) = n$. Самым простым методом нахождения сингулярного разложения является сведение задачи нахождения сингулярного разложения к задаче нахождения спектрального разложения.
В разделах (\ref{theory:beltrami}), (\ref{theory:sylvester}), (\ref{theory:autonne}) было доказано, что любая квадратная невырожденная матрица $A$ может быть представлена в виде
\begin{equation} \label{A=USV^T}
A = U\Sigma V^T,
\end{equation}
где $U,V$ -- ортогональные матрицы, а $\Sigma$ -- диагональная матрица.
Составим матрицу $B = A^TA$. Она симметричная и положительно определенная, следовательно имеет полный набор действительных положительных собственных значений. Учитывая, что $A = U\Sigma V^T$, перепишем $B$ в виде \[ B = (U\Sigma V^T)^TU\Sigma V^T = V \Sigma^T \underbrace{U^TU}_{I} \Sigma V^T = V \Sigma^2 V^T. \] Без особого труда можно заметить, что $B = V \Sigma^2 V^T$ -- спектральное разложение матрицы $B$, где $V$ -- матрица, состоящая из столбцов собственных векторов, а $\Lambda := \Sigma^2$ -- диагональная матрица, состоящая из собственных значений.
Таким образом, найдя спектральное разложение матрицы $B$ (например, двусторонним методом вращений Якоби, рассмотренном на лекции), мы получим матрицы $V,\Sigma^2$. Матрица $\Sigma = \sqrt{\Sigma^2}$ является диагональной матрицей с элементами равными корню квадратному из элементов матрицы $\Sigma^2$.
Оставшаяся матрица $U$ находится явным способом из (\ref{A=USV^T}): \[ U = AV\Sigma^{-1}, \] где $\Sigma^{-1}$ -- диагональная матрица с элементами равными обратным значениям элементов матрицы $\Sigma$.
\subsubsection{Исходный код}
\lstinputlisting{source_code/svd_math_1.m} \label{source:svd_math_1}
\subsubsection{Для матриц $m\times n$, $r = \min{(m,n)}$}
В разделах (\ref{theory:sing_vectors}), (\ref{theory:sing_decomp}) было неявно доказано, что любая матрица $A_{m \times n}$ такая, что $r = \rank(A) = \min{(m,n)}$ может быть разложена в произведение $A = U\Sigma V^T$.
Основное отличие этого подхода от предыдущего заключается в том, что используемый нами двусторонний алгоритм вращений Якоби для нахождения спектрального разложения не работает для вырожденных матриц.
Предположим, что матрица $A$ имеет размерность $m \times n$ и $\rank(A) = m$. Если мы воспользуемся вышеизложенным алгоритмом, то получим следующий результат: \[ (A^T)_{n\times m} A_{m\times n} = B_{n\times n}, \] где $B_{n\times n}$ -- симметрическая положительная полуопределенная матрица, причем $\rank(B_{n\times n}) = m, \text{ } m < n$, т.е. матрица $B$ является вырожденной.
Чтобы обойти эту проблему достаточно рассмотреть два случая:
\begin{itemize}
\item Если $\rank(A_{m \times n}) = n$, то $B = A^TA$ и весь алгоритм аналогичен вышеизложенному.
\item Если $\rank(A_{m \times n}) = m$, то $B = AA^T$ и алгоритм слегка изменяется.
Учитывая (\ref{A=USV^T}), найдем представление $B$: \[ B = U\Sigma V^T(U\Sigma V^T)^T = U \Sigma \underbrace{VV^T}_{I} \Sigma^T U^T = U \Sigma^2 U^T. \] Дальнейшие вычисления аналогичны предыдущему алгоритму, за исключением того, что, разложив $B$ с помощью спектрального разложения, мы получим $U,\Sigma$ и после этого из (\ref{A=USV^T}) найдем $V$ с помощью явной подстановки $V = (\Sigma^{-1}U^TA)^T$.
\end{itemize}
Рассмотрим правило немного подробнее. \\
Если $\rank(A_{m \times n}) = n$, то \[ B = (A^T)_{n \times m}A_{m \times n} = B_{n \times n}, \] где $B$ -- невырожденная симметричная матрица. Сингулярное разложение в таком случае будет выглядеть так: \[ A_{m \times n} = U_{m \times n} \Sigma_{n \times n} (V^T)_{n \times n}. \]
Если $\rank(A_{m \times n}) = m$, то \[ B = A_{m \times n}(A^T)_{n \times m} = B_{m \times m}, \] где $B$ -- невырожденная симметричная матрица. Сингулярное разложение в таком случае будет выглядеть так: \[ A_{m \times n} = U_{m \times m} \Sigma_{m \times m} (V^T)_{m \times n}. \]
\subsubsection{Исходный код}
\lstinputlisting{source_code/svd_math_2.m} \label{source:svd_math_2}
\subsubsection{Для матриц $m\times n$, $r < \min{(m,n)}$}
В разделах (\ref{theory:sing_vectors}), (\ref{theory:sing_decomp}) было неявно доказано, что любая матрица $A_{m \times n}$ такая, что $r = \rank(A) < \min{(m,n)}$, может быть разложена в произведение $A = U\Sigma V^T$.
Используемый нами двусторонний алгоритм вращений Якоби для нахождения спектрального разложения не работает для вырожденных матриц. И в данном случае никак не получится создать невырожденную симметричную матрицу $B$. Поэтому единственным решением возникшей проблемы является использование другого алгоритма вычисления спектрального разложения, который будет правильно обрабатывать вырожденность матрицы.
Из вышесказанного заключаем, что алгоритм вычисления SVD разложения для всех возможных действительных матриц будет аналогичен (\ref{source:svd_math_1}), за тем исключением, что функция \mcode{eig(B)} вычисляет спектральное разложение симметричной положительной полуопределенной матрицы $B$ с учетом вырожденности.
После изучения документации языка Matlab становится ясно, какие алгоритмы спектрального разложения используются на практике. По умолчанию, если матрица $B$ симметричная, то используется разложение Холецкого для решения задачи собственных значений и векторов, в противном случае используется обобщенное разложение Шура.
В конкретной задаче известно, что матрица $B$ -- симметричная и положительно полуопределенная, следовательно для нахождения собственных векторов и значений можно использовать разложение Холецкого. Более подробно об этом можно почитать в \cite{citation4}, \cite{citation5}.
\subsubsection{Исходный код} \label{source:svd_math}
\lstinputlisting{source_code/svd_math.m}
\subsection{Виды SVD разложений}
Разница SVD разложений заключается в размерностях получающихся матриц $U,\Sigma,V$.
\subsubsection{\quotes{Тонкий} SVD}
Данный тип разложения использовался в (\ref{source:svd_math_1}) и (\ref{source:svd_math_2}). \\
Предположим, что имеется матрица $A_{m \times n}$, причем $m \ge n$. Тогда сингулярное разложение можно представить в виде произведения матриц следующих размерностей: \[ A_{m \times n} = U_{m \times n} \Sigma_{n \times n} (V^T)_{n \times n}. \] Представим данное разложение графически:
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (2,3) node[pos=.5] {$A_{m \times n}$};
\node at(2.5,1.5){=};
\draw (3,0) rectangle (5,3) node[pos=.5] {$U_{m \times n}$};
\draw (5.5,1) rectangle (7.5,3) node[pos=.5] {$\Sigma_{n \times n}$};
\draw (8,1) rectangle (10,3) node[pos=.5] {$(V^T)_{n \times n}$};
\end{tikzpicture}
\end{align*}
Теперь предположим, что имеется матрица $A_{m \times n}$, причем $m \le n$. Тогда графическое представление будет выглядеть таким образом:
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (3,2) node[pos=.5] {$A_{m \times n}$};
\node at(3.5,1){=};
\draw (4,0) rectangle (7,2) node[pos=.5] {$U_{m \times n}$};
\draw (7.5,-1) rectangle (10.5,2) node[pos=.5] {$\Sigma_{n \times n}$};
\draw (11,-1) rectangle (14,2) node[pos=.5] {$(V^T)_{n \times n}$};
\end{tikzpicture}
\end{align*}
Однако в (\ref{jacobi_svd}) при условии, что $m \le n$ результирующие матрицы имеют вид: \[ A_{m \times n} = U_{m \times m} \Sigma_{m \times m} (V^T)_{m \times n}. \]
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (3,2) node[pos=.5] {$A_{m \times n}$};
\node at(3.5,1){=};
\draw (4,0) rectangle (6,2) node[pos=.5] {$U_{m \times m}$};
\draw (6.5,0) rectangle (8.5,2) node[pos=.5] {$\Sigma_{m \times m}$};
\draw (9,-1) rectangle (11,2) node[pos=.5] {$(V^T)_{m \times n}$};
\end{tikzpicture}
\end{align*}
\subsubsection{\quotes{Полный} SVD}
Данный тип разложения использовался в (\ref{source:svd_math}).\\
Предположим, что имеется матрица $A_{m \times n}$. Тогда сингулярное разложение можно представить в виде произведения матриц следующих размерностей: \[ A_{m \times n} = U_{m \times m} \Sigma_{m \times n} (V^T)_{n \times n}. \] Тогда графически это разложение можно представить двумя способами в зависимости от того, что больше: $m$ или $n$.
Если $m < n$, то
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (3,2) node[pos=.5] {$A_{m \times n}$};
\node at(3.5,1){=};
\draw (4,0) rectangle (6,2) node[pos=.5] {$U_{m \times m}$};
\draw (6.5,0) rectangle (9.5,2) node[pos=.5] {$\Sigma_{m \times n}$};
\draw (10,-1) rectangle (13,2) node[pos=.5] {$(V^T)_{n \times n}$};
\end{tikzpicture}
\end{align*}
Если $m > n$, то
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (2,3) node[pos=.5] {$A_{m \times n}$};
\node at(2.5,1.5){=};
\draw (3,0) rectangle (6,3) node[pos=.5] {$U_{m \times m}$};
\draw (6.5,0) rectangle (8.5,3) node[pos=.5] {$\Sigma_{m \times n}$};
\draw (9,1) rectangle (11,3) node[pos=.5] {$(V^T)_{n \times n}$};
\end{tikzpicture}
\end{align*}
\subsubsection{\quotes{Сжатый} SVD}
Данный тип разложения является очень важным и часто используемым по двум причинам:
\begin{enumerate}
\item алгоритм имеет меньшую вычислительную сложность
\item алгоритм экономит память
\end{enumerate}
Для упрощения понимания рассмотрим матрицу $A_{m \times n}$, причем $m > n$. Предположим, что $\rank(A) = r$. Тогда сингулярное разложение в \quotes{полном} виде будет выглядеть следующим образом:
\begin{align*}
\begin{tikzpicture}[
baseline,
mymat/.style={
matrix of math nodes,
ampersand replacement=\&,
left delimiter=(,
right delimiter=),
nodes in empty cells,
nodes={outer sep=-\pgflinewidth,text depth=0.4ex,text height=1.9ex,text width=1.1em}
}
]
\begin{scope}[every right delimiter/.style={xshift=-3ex}]
\matrix[mymat] (matu)
{
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
};
\node
at ([shift={(3pt,-7pt)}]matu-3-2.west)
{$\cdots$};
\node
at ([shift={(3pt,-7pt)}]matu-3-5.west)
{$\cdots$};
\foreach \Columna/\Valor in {1/1,3/r,4/{r+1},6/m}
{
\draw
(matu-1-\Columna.north west)
rectangle
([xshift=4pt]matu-6-\Columna.south west);
\node[above]
at ([xshift=2pt]matu-1-\Columna.north west)
{$u_{\Valor}$};
}
\draw[decorate,decoration={brace,mirror,raise=3pt}]
(matu-6-1.south west) --
node[below=4pt] {$\Mcol(A)$}
([xshift=4pt]matu-6-3.south west);
\draw[decorate,decoration={brace,mirror,raise=3pt}]
(matu-6-4.south west) --
node[below=4pt] {$\Mnull(A^T)$}
([xshift=4pt]matu-6-6.south west);
\end{scope}
\begin{scope}[every right delimiter/.style={xshift=-0.5ex}]
\matrix[mymat,right=10pt of matu] (matsigma)
{
\sigma_{1} \& \& \& \& \& \\
\& \ddots \& \& \& 0 \& \\
\& \& \sigma_{r} \& \& \& \\
\& \& \& 0 \& \& \\
\& 0 \& \& \& \ddots \& \\
\& \& \& \& \& 0\\
0\& \hdots \& \hdots \& \hdots \& \hdots \& 0 \\
\vdots \& \& \& \& \& \vdots \\
0 \& \hdots \& \hdots \& \hdots \& \hdots \& 0 \\
};
\draw[gray] (matsigma-3-1.south west) -- (matsigma-3-6.south east);
\draw[gray] (matsigma-1-3.north east) -- (matsigma-6-3.south east);
\draw[gray] (matsigma-6-1.south west) -- (matsigma-6-6.south east);
\end{scope}
\begin{scope}[every right delimiter/.style={xshift=-1ex}]
\matrix[mymat,right=22pt of matsigma] (matv)
{
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
\& \& \& \& \& \\
};
\foreach \Fila/\Valor in {1/1,3/r,4/{r+1},6/n}
{
\draw
([yshift=-6pt]matv-\Fila-1.north west)
rectangle
([yshift=-10pt]matv-\Fila-6.north east);
\node[right=6pt] % двигаются v_i
at ([yshift=-8pt]matv-\Fila-6.north east)
{$v^{T}_{\Valor}$};
}
\draw[decorate,decoration={brace,raise=32pt}]
([yshift=-7pt]matv-1-6.north east) --
node[right=34pt] {$\Mrow(A)$}
([yshift=-11pt]matv-3-6.north east);
\draw[decorate,decoration={brace,raise=32pt}]
([yshift=-7pt]matv-4-6.north east) --
node[right=34pt] {$\Mnull(A)$}
([yshift=-11pt]matv-6-6.north east);
\end{scope}
\end{tikzpicture}
\end{align*}
Из приведенной иллюстрации легко можно заметить, что матрицу $\Sigma$ можно представить в виде блочной матрицы
\begin{equation}
\Sigma_{m \times n} = \left(
\begin{array}{r@{}c|c@{}l}
& \Sigma_{r \times r} & O_{r \times (n-r)} & \\
\cline{1-4}
& O_{(n-r) \times r} & O_{(n-r) \times (n-r)} & \\
\cline{1-4}
& O_{(m-n+r) \times r} & O_{(m-n+r) \times (n-r)} &
\end{array}
\right),
\end{equation}
поэтому матрицу $\Sigma_{m \times n}$ можно задать матрицей меньшей размерности $\Sigma_{r \times r}$. Векторы $u_{r+1},...,u_{m}$ задают ядро пространства строк, а векторы $v_{r+1},...,v_{n}$ задают ядро пространства столбцов, поэтому их также можно не считать. Т.о. матрицы $U,\Sigma,V$ свелись к матрицам $\widetilde{U},\widetilde{\Sigma},\widetilde{V}$:
\begin{equation} \nonumber
\begin{aligned}
U_{m \times m} & \rightarrow \widetilde{U}_{m \times r} \\
\Sigma_{m \times n} & \rightarrow \widetilde{\Sigma}_{r \times r} \\
V_{n \times n} & \rightarrow \widetilde{V}_{n \times r}
\end{aligned}
\end{equation}
Сингулярное разложение можно представить в виде произведения матриц следующих размерностей: \[ A_{m \times n} = U_{m \times m} \Sigma_{m \times n} (V^T)_{n \times n} = U_{m \times r} \Sigma_{r \times r} (V^T)_{r \times n}. \]
Графически SVD-разложение теперь можно представить следующими двумя изображениями:
Если $m < n$, то
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (3,2) node[pos=.5] {$A_{m \times n}$};
\node at(3.5,1){=};
\draw[dashed] (4,0) rectangle (6,2) node[pos=.5] {};
\draw[dashed] (6.5,0) rectangle (9.5,2) node[pos=.5] {};
\draw[dashed] (10,-1) rectangle (13,2) node[pos=.5] {};
\draw[red] (4,0) rectangle (5.5,2) node[pos=.5] {$U_{m \times r}$};
\draw[red] (6.5,0.5) rectangle (8,2) node[pos=.5] {$\Sigma_{r \times r}$};
\draw[red] (10,0.5) rectangle (13,2) node[pos=.5] {$(V^T)_{r \times n}$};
\end{tikzpicture}
\end{align*}
Если $m > n$, то
\begin{align*}
\begin{tikzpicture}
\draw (0,0) rectangle (2,3) node[pos=.5] {$A_{m \times n}$};
\node at(2.5,1.5){=};
\draw[dashed] (3,0) rectangle (6,3) node[pos=.5] {};
\draw[dashed] (6.5,0) rectangle (8.5,3) node[pos=.5] {};
\draw[dashed] (9,1) rectangle (11,3) node[pos=.5] {};
\draw[red] (3,0) rectangle (4.5,3) node[pos=.5] {$U_{m \times r}$};
\draw[red] (6.5,1.5) rectangle (8,3) node[pos=.5] {$\Sigma_{r \times r}$};
\draw[red] (9,1.5) rectangle (11,3) node[pos=.5] {$(V^T)_{r \times n}$};
\end{tikzpicture}
\end{align*}
\newpage % =======================================================================================================================================
\section{Односторонний алгоритм вращений Якоби} \label{jacobi_svd}
\subsection{Описание алгоритма}
Пусть дана матрица $A$ размерности $m \times n$, причем $m \ge n$ и $\rank{(A)} = n$. Если условие $m \ge n$ не выполняется, то находится сингулярное разложение матрицы $B = A^T$: \[ B = \bar{U}\bar{\Sigma}\bar{V}^T, \]
после чего оригинальное разложение матрицы $A$ получается транспонированием:
\[ A = B^T = (\bar{U}\bar{\Sigma}\bar{V}^T)^T = \underbrace{\bar{V}}_{U}\underbrace{\bar{\Sigma}}_{\Sigma}\underbrace{\bar{U}^T}_{V^T}. \]
Условие $\rank{(A)} = n$ является необходимым.
Предположим, что нам известна матрица $M$ размерности $2 \times 2$ такая, что
\begin{equation} \nonumber
M^TM = \begin{pmatrix}
\alpha & \gamma \\
\gamma & \beta
\end{pmatrix}.
\end{equation}
Симметричность выполняется всегда: $(M^TM)^T = M^T(M^T)^T = M^TM$. В таком случае с помощью вращений Якоби $\Theta (M^T M) \Theta^T$,
\begin{equation} \nonumber
\Theta = \begin{pmatrix}
\cos{\theta} & -\sin{\theta} \\
\sin{\theta} & \cos{\theta}
\end{pmatrix} := \begin{pmatrix}
c & -s \\
s & c
\end{pmatrix}
\end{equation}
симметричную матрицу $M^TM$ можно привести к диагональному виду:
\[ D = \Theta M^T M \Theta^T = (M \Theta^T)^T (M \Theta^T). \]
Перемножив матрицы $\Theta M^T M \Theta^T$, получим следующее выражение:
\begin{equation} \nonumber
\begin{pmatrix}
\alpha c^2 - 2\gamma cs + \beta s^2 & -\gamma s^2 + (\alpha - \beta) cs + \gamma c^2 \\
-\gamma s^2 + (\alpha - \beta) cs + \gamma c^2 & \beta c^2 +2\gamma cs + \alpha s^2
\end{pmatrix} = \begin{pmatrix}
d_1 & 0 \\
0 & d_2
\end{pmatrix}.
\end{equation}
Оно имеет всего одну степень свободы $\theta$, поэтому достаточно решить единственное уравнение:
\[ -\gamma s^2 + (\alpha - \beta) cs + \gamma c^2 = 0, \]
чтобы задать матрицу вращений Гивенса $\Theta$. Поделив обе части уравнения на $-\gamma c^2$, получим:
\[ t^2 + \dfrac{(\beta - \alpha)}{\gamma}t - 1 = 0; \text{ } t := \dfrac{s}{c} = \tan{\theta}. \]
\begin{equation} \label{eq:quadratic}
t^2 + 2\zeta t - 1 = 0; \text{ } \zeta = \dfrac{(\beta - \alpha)}{2\gamma}.
\end{equation}
Решая уравнение, получаем корни \[ t = -\zeta \pm \sqrt{1 + \zeta^2}. \]
Заметим, что \[ (-\zeta \pm \sqrt{1 + \zeta^2})(\zeta \pm \sqrt{1 + \zeta^2}) = -\zeta^2 + 1 + \zeta^2 = 1, \]
откуда следует, что \[ t = -\zeta \pm \sqrt{1 + \zeta^2} = \dfrac{1}{\zeta \pm \sqrt{1 + \zeta^2}}. \]
Возьмем в качестве корня уравнения \[ t = \dfrac{\sign{(\zeta)}}{|\zeta|+\sqrt{1+\zeta^2}}. \]
Из тангенса выразим необходимые нам синус и косинус:
\begin{equation} \nonumber
\left\{
\begin{aligned}
c & = \dfrac{1}{\sqrt{1+t^2}} \\
s & = ct
\end{aligned}
\right.
\end{equation}
Таким образом мы полностью задали $\Theta$.
Введем замену $U_0 := A; V_0 := I_{n \times n}$. Предположим, что $R_0$ -- какая-то матрица вращений Гивенса. Тогда, вращая матрицу $U_0$ последовательно с помощью $R_0$ и $R_0^{-1}$, мы получим исходную матрицу $U_0$:
\[ A \equiv U_0 = U_0R_0^{-1}R_0. \]
Воспользовавшись свойством ортогональных матриц ($R^{-1} = R^T$), перепишем равенство в виде:
\[ A = U_0R_0^TR_0. \]
Сделаем замену
\begin{equation} \nonumber
A = \underbrace{U_0R_0^T}_{U_1}\underbrace{R_0}_{V_1} = U_1V_1.
\end{equation}
Продолжим рассматривать уже $U_1$. Проведем аналогичные рассуждения: предположим, что $R_1$ -- какая-то матрица вращений Гивенса. Тогда, вращая матрицу $U_1$ последовательно с помощью $R_1$ и $R_1^{-1}$, мы получим исходную матрицу $U_1$:
\[ A \equiv U_1 = U_1R_1^{-1}R_1. \]
Опять сделаем замену
\begin{equation} \nonumber
A = \underbrace{U_1R_1^T}_{U_2}\underbrace{R_1V_1}_{V_2} = \underbrace{U_0R_0^TR_1^T}_{U_2}\underbrace{R_1R_0}_{V_2} = U_2V_2.
\end{equation}
Продолжая аналогичные рассуждения, напишем формулу для $n$-ой итерации:
\begin{equation} \nonumber
A = \underbrace{U_{n-1}R_{n-1}^T}_{U_n}\underbrace{R_{n-1}V_{n-1}}_{V_n} = \underbrace{U_0R_0^T...R_{n-1}^T}_{U_n}\underbrace{R_{n-1}...R_0}_{V_n} = U_nV_n.
\end{equation}
Идея алгоритма заключается в том, что в бесконечном цикле неявно конструируется матрица $A^TA$, выбирается пара индексов таких, что $i < j$, строится $2 \times 2$ матрица из элементов, стоящих на пересечении $i,j$-ых строк и столбцов матрицы $A^TA$, после чего она приводится к диагональному виду с помощью вращений Якоби. Итерации повторяются до тех пор, пока не достигается лимит итераций или достаточная точность.
С условием остановки возникает проблема: т.к. мы не конструируем матрицу $U^TU$ в явном виде, а лишь вычисляем 4 ее элемента, то нельзя ввести норму как сумму квадратов внедиагональных элементов. В \cite{citation3} в качестве нормы предлагают брать \[ ||U^TU|| = \dfrac{|c|}{\sqrt{\alpha \beta}}, \] поэтому критерием остановки является \[ \dfrac{|c|}{\sqrt{\alpha \beta}} \le \varepsilon. \]
После окончания итерирования у нас имеется 2 матрицы:
\begin{equation} \nonumber
\begin{aligned}
U_n &= U_0R_0^T...R_{n-1}^T \\
V_n &= R_{n-1}...R_0
\end{aligned}
\end{equation}
Чтобы из матриц $U_n,V_n$ получить сингулярное разложение, необходимо матрицу $U_n$ представить в виде $U_n = \bar{U}_n\Sigma$, где $\Sigma$ -- диагональная матрица, элементы которой заданы следующим образом: $\sigma_i = ||U_n^{(i)}||$, где $U_n^{(i)}$ -- $i$-ый столбец матрицы $U_n$. В таком случае $\bar{U}_n^{(i)} = \dfrac{U_n^{(i)}}{||U_n^{(i)}||}$.
Уточним, как вычисляются элементы матрицы $U^TU$:
\begin{equation} \nonumber
B = U^TU = \begin{pmatrix}
* & \cdots & * & \cdots & * & \cdots & * \\
\vdots & \ddots & \vdots & & \vdots & & \vdots \\
* & \cdots & b_{ii} & \cdots & b_{ij} & \cdots & * \\
\vdots & & \vdots & \ddots & \vdots & & \vdots \\
* & \cdots & b_{ji} & \cdots & b_{jj} & \cdots & * \\
\vdots & & \vdots & & \vdots & \ddots & \vdots \\
* & \cdots & * & \cdots & * & \cdots & *
\end{pmatrix}_{n \times n}
\end{equation}
$b_{ii}$ вычисляется как скалярное произведение $i$-ой строки матрицы $U^T$ и $i$-ого столбца матрицы $U$, но т.к. $i$-ый столбец матрицы $U$ после транспонирования становится $i$-ой строкой матрицы $U^T$, то вычисление элемента $b_{ii}$ можно представить в виде: \[ \alpha := b_{ii} = \sum_{k=1}^{m} u_{k,i}^2. \]
Аналогично рассуждая, получаем формулу для вычисления $b_{jj}$: \[ \beta := b_{jj} = \sum_{k=1}^{m} u_{k,j}^2. \]
Т.к. $b_{ij} = b_{ji}$, то для определенности будем вычислять $b_{ij}$. $b_{ij}$ вычисляется как скалярное произведение $i$-ой строки матрицы $U^T$ и $j$-ого столбца матрицы $U$, а т.к. $i$-я строка матрицы $U^T$ -- $i$-ый столбец матрицы $U$, то вычисление $b_{ij}$ можно представить в следующем виде: \[ \gamma := b_{ij} = b_{ji} = \sum_{k=1}^{m} u_{ki} u_{kj}. \]
\subsection{Исходный код}
В следующем листинге приведен код одностороннего алгоритма вращений Якоби для вычисления сингулярного разложения, написанный на языке Matlab.
\lstinputlisting{source_code/svd_jacobi.m}
\subsection{Анализ работы алгоритма}
\subsubsection{Анализ решения уравнения \ref{eq:quadratic}}
Во время вывода алгоритма первая произвольность в выборе какого-либо параметра возникает во время решения квадратного уравнения \[ t^2 + 2\zeta t - 1 = 0. \]
Как было показано, корней у уравнения 2, и каждое можно представить двумя способами:
\begin{equation} \nonumber
\left[
\begin{aligned}
t_{1,2} &= -\zeta \pm \sqrt{1 + \zeta^2}, \\
t_{1,2} &= \dfrac{1}{\zeta \pm \sqrt{1 + \zeta^2}}.
\end{aligned}
\right.
\end{equation}
Анализируя необходимое количество итераций для сходимости алгоритма и точность полученных вычислений, выберем наилучшее решение и его представление.
\begin{figure}[h!]
\centering
\includegraphics[height=0.89\textheight]{img/solution_analysis/barchart_rotated.png}
\caption{Диаграмма скорости сходимости в зависимости от размерности матрицы}
\label{barchartanalysis}
\end{figure}
На Рис.~\ref{barchartanalysis} представлен bar chart, с помощью которого можно проанализировать поведение алгоритма в зависимости от того, какой вариант корня был выбран. Для построения данной диаграммы создавались случайные квадратные гауссовы матрицы $A$ нужной размерности, после чего производилось сингулярное разложение с каждым из корней в обоих представлениях и составлялась статистика работы.
Вдоль координаты $x$ заданы размерности матрицы $A$, вдоль координаты $y$ -- среднее количество итераций, которое понадобилось для сходимости алгоритма. Для того, чтобы получить среднее количество итераций, которое требуется для сходимости алгоритма, необходимо произвести вычисления несколько раз. Число, стоящее над каждым набором столбцов, показывает, сколько вычислений было произведено перед тем, как составить результат.
Столбец, повернутый в отрицательную область показывает, что ни одно вычисление не сошлось. Если хотя бы одно вычисление сошлось, то столбец будет в положительной зоне.
Проанализировав данный график, приходим к выводу, что упорядоченный список (лучший $\rightarrow$ худший) выглядит следующим образом:
\begin{enumerate}
\item $t = \dfrac{\sign{(\zeta)}}{|\zeta| - \sqrt{1 + \zeta^2}}$. При таком выборе алгоритм сходится стандартно за 2-3 итерации.
\item $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$. В данном случае график сходимости больше напоминает логарифмическую кривую.
\item $t = -\sign{(\zeta)}(|\zeta| + \sqrt{1 + \zeta^2})$. Ярко выраженного шаблона при использовании такой формулы график сходимости не имеет. Однако график отдаленно напоминает зашумленную линейную функцию.
\item $t = -\sign{(\zeta)}(|\zeta| - \sqrt{1 + \zeta^2})$. Как видно из диаграммы, при использовании такой формулы алгоритм перестает сходиться, начиная с 7 размерности.
\end{enumerate}
Однако необходимо также учесть, насколько хороший результат получает алгоритм. Сходимость к какой-то матрице сама по себе ничего не значит, необходимо, чтобы алгоритм возвращал такие значения, которые соответствуют условиям.
Параллельно с вышеприведенной диаграммой на Рис.~\ref{barchartanalysis} строились и круговые диаграммы, которые показывают процент сходимости алгоритма с каждой используемой формулой и качество сходимости.
\begin{figure}[h!]
\centering
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/convergence_4.png}
\caption{Процент сходимости при\\ $t = \dfrac{\sign{(\zeta)}}{|\zeta| - \sqrt{1 + \zeta^2}}$}
\label{fig:t_4_analysis_a}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/correctness_4.png}
\caption{Процент правильности при\\ $t = \dfrac{\sign{(\zeta)}}{|\zeta| - \sqrt{1 + \zeta^2}}$}
\label{fig:t_4_analysis_b}
\end{subfigure}
\caption{Анализ сходимости при $t = \dfrac{\sign{(\zeta)}}{|\zeta| - \sqrt{1 + \zeta^2}}$}
\label{fig:t_4_analysis}
\end{figure}
Проанализируем сперва самое \quotes{лучшее} представление $t$. На Рис.~\ref{fig:t_4_analysis} представлены две круговые диаграммы, которые показывают процент сходимости алгоритма и качество сходимости при $t = \dfrac{\sign{(\zeta)}}{|\zeta| - \sqrt{1 + \zeta^2}}$. Из Рис.~\ref{fig:t_4_analysis_a} видно, что алгоритм сошелся в $100\%$ проведенных испытаниях, однако из Рис.~\ref{fig:t_4_analysis_b} становится ясно, что алгоритм выдает неправильные значения. Причиной тому может быть неподходящее условие выхода, которое указывает нам на сходимость алгоритма.
Продолжим аналогично анализировать и другие представления $t$.
\begin{figure}[h!]
\centering
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/convergence_3.png}
\caption{Процент сходимости при\\ $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$}
\label{fig:t_3_analysis_a}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/correctness_3.png}
\caption{Процент правильности при\\ $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$}
\label{fig:t_3_analysis_b}
\end{subfigure}
\caption{Анализ сходимости при $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$}
\label{fig:t_3_analysis}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/optimal.png}
\end{subfigure}
\caption{Оптимальный выбор $t$}
\label{fig:t_optimal}
\end{figure}
\begin{figure}[h!]
\centering
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/convergence_1.png}
\caption{Процент сходимости при\\ $t = -\sign{(\zeta)}(|\zeta| + \sqrt{1 + \zeta^2})$}
\label{fig:t_1_analysis_a}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/correctness_1.png}
\caption{Процент правильности при\\ $t = -\sign{(\zeta)}(|\zeta| + \sqrt{1 + \zeta^2})$}
\label{fig:t_1_analysis_b}
\end{subfigure}
\caption{Анализ сходимости при $t = -\sign{(\zeta)}(|\zeta| + \sqrt{1 + \zeta^2})$}
\label{fig:t_1_analysis}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/convergence_2.png}
\caption{Процент сходимости при\\ $t = -\sign{(\zeta)}(|\zeta| - \sqrt{1 + \zeta^2})$}
\label{fig:t_2_analysis_a}
\end{subfigure}
\begin{subfigure}{.495\textwidth}
\includegraphics[width=\textwidth]{img/solution_analysis/correctness_2.png}
\caption{Процент правильности при\\ $t = -\sign{(\zeta)}(|\zeta| - \sqrt{1 + \zeta^2})$}
\label{fig:t_2_analysis_b}
\end{subfigure}
\caption{Анализ сходимости при $t = -\sign{(\zeta)}(|\zeta| - \sqrt{1 + \zeta^2})$}
\label{fig:t_2_analysis}
\end{figure}
Из Рис.~\ref{fig:t_3_analysis_a} видно, что алгоритм сошелся в $100\%$ проведенных испытаниях, причем из Рис.~\ref{fig:t_3_analysis_b} становится ясно, что алгоритм выдает правильные значения в $>99\%$ случаев; $<1\%$ можно списать на вычислительную погрешность, т.к. вычисление признавалось ошибочным, если хотя бы один элемент полученной матрицы отличался от исходной более чем на $10^{-8}$. Т.о. мы приходим к выводу, что алгоритм при $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$ работает и быстро, и точно.
Продолжая исследовать остальные два варианта представления $t$ Рис.~\ref{fig:t_1_analysis}, Рис.~\ref{fig:t_2_analysis}, мы приходим к выводу, что $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$ является наиболее быстрым и точным.
Т.к. на Рис.~\ref{barchartanalysis} показаны средние значения, проанализируем Рис.~\ref{fig:t_optimal}, на котором выводится процент оптимальности каждого из представлений. Как можно сразу заметить, представление $t = \dfrac{\sign{(\zeta)}}{|\zeta| + \sqrt{1 + \zeta^2}}$ лидирует по всем критериям, поэтому и в выводе алгоритма, и в исходном коде используется именно такое представление.
\newpage % =======================================================================================================================================
\section{\quotes{Стандартный} алгоритм}
\quotes{Стандартный} алгоритм -- самый часто используемый алгоритм, который реализован в Matlab и библиотеке линейной алгебры Lapack. Данный алгоритм является улучшением алгоритма, разработанного математиками Golub и Van Loan. Улучшенный алгоритм представлен математиками J.Demmel и W.Kahan в \cite{citation6}.
\subsection{Описание алгоритма}
\quotes{Стандартный} алгоритм происходит в два этапа:
\begin{enumerate}
\item Сперва происходит бидиагонализация матрицы с помощью преобразований Хаусхолдера;
\item Затем с помощью итерационного алгоритма происходит нахождение сингулярного разложения бидиагональной матрицы.
\end{enumerate}
Сложность создания итерационного алгоритма заключается в том, что необходимо минимизировать одновременно и погрешность вычислений, и количество итераций. В оригинале данная задача решается тем, что на каждом шаге итерации выбирается, какой из двух алгоритмов: Golub \& Van Loan's algorithm или Demmel \& Kahan's algorithm подходит для данного случая, но в данной работе будет рассмотрен упрощенный алгоритм: Demmel \& Kahan's zero-shift algorithm.
\subsubsection{Преобразование Хаусхолдера}
Преобразование Хаусхолдера -- линейное преобразование $H_u$ векторного пространства $V$, которое описывает его отображение относительно гиперплоскости, которая проходит через начало координат.
\begin{wrapfigure}{l}{0pt}
\begin{tikzpicture}[scale=4]
\draw[<->] (0,1.5) -- (0,0) -- (1.5,0);
\draw[thick,->] (0,0) -- (0.638,1.183) node[anchor=base,black] at (0.27,0.73) {$\mathbf{x}$};
\draw[red,thick,dashed] (-0.1,-0.08333) -- (1.5,1.25) node[anchor=base] at (1.2,1.15) {$L$};
\draw[blue,thick,->] (0,0) -- (-0.1921,0.2305) node[anchor=base] at (-0.2,0.29) {$\mathbf{u}$};
\draw[green,thick,->] (0,0) -- (1.2787,0.4142) node[anchor=base,black] at (1.21,0.13) {$\mathbf{y}=H_{\mathbf{u}}(\mathbf{x})$};
\draw[blue,dashed] (1.2787,0.4142) -- (0.638,1.183);
\draw[blue,thick,->] (0.9583,0.7986) -- (0.638,1.183) node[anchor=base] at (0.85,1) {$\mathbf{w}$};
\end{tikzpicture}
\hspace{10pt}
\caption{Отражение}
\label{fig:hauseholder}
\end{wrapfigure}
Предположим, что имеется радиус-вектор $\mathbf{x}$ в $n$-мерном пространстве, который необходимо отразить относительно гиперплоскости $L$, которая задает $(n-1)$-мерное подпространство, проходящее через начало координат. Такую гиперплоскость можно задать с помощью вектора-нормали $\mathbf{u}$.
Т.е. задача состоит в том, чтобы составить линейное преобразование $H_\mathbf{u}$ векторного пространства $V$, которое бы отражало радиус-вектор $\mathbf{x}$ относительно гиперплоскости, заданной нормалью $\mathbf{u}$: \[\mathbf{y}=H_{\mathbf{u}}(\mathbf{x}).\]
Из рис.\ref{fig:hauseholder} понятно, что вектор $\mathbf{w}$ определяется следующим образом: \[\mathbf{w} = (\mathbf{x},\mathbf{u})\mathbf{u}.\]
Также можно заметить, что \[\mathbf{y} = \mathbf{x} - 2\mathbf{w} = \mathbf{x} - 2(\mathbf{x},\mathbf{u})\mathbf{u}.\]
Учитывая, что скалярное произведение двух столбцов $(a,b)$ можно представить в виде произведения строки и столбца $a^Tb$, перепишем полученное равенство: \[\mathbf{y}=\mathbf{x}-2(\mathbf{x}^T\mathbf{u})\mathbf{u}.\]
В данном выражении $(\mathbf{x},\mathbf{u})$ -- константа, поэтому можно записать равенство в следующем виде: \[ \mathbf{y} = \mathbf{x}-2\mathbf{u}(\mathbf{x}^T\mathbf{u}). \]
Воспользовавшись тем свойством, что $a^Tb \equiv b^Ta$, получим \[ \mathbf{y} = \mathbf{x}-2\mathbf{u}(\mathbf{x}^T\mathbf{u}) = \mathbf{x}-2\mathbf{u}(\mathbf{u}^T\mathbf{x}), \]
откуда, вынося $\mathbf{x}$, получим матрицу Хаусхолдера:
\[ H_{\mathbf{u}}\mathbf{x} = \mathbf{x}-2\mathbf{u}\mathbf{u}^T\mathbf{x} = \underbrace{(I-2\mathbf{u}\mathbf{u}^T)}_{H_{\mathbf{u}}}\mathbf{x}.\]
\[H_{\mathbf{u}} = I-2\mathbf{u}\mathbf{u}^T.\]
Важным свойством матрицы преобразования Хаусхолдера является тот факт, что она ортонормирована:
\begin{equation} \nonumber
\begin{aligned}
H_{\mathbf{u}}H^T_{\mathbf{u}} &= (I - 2\mathbf{u}\mathbf{u}^T)(I - 2\mathbf{u}\mathbf{u}^T)^T = I - 2\mathbf{u}\mathbf{u}^T - 2\mathbf{u}\mathbf{u}^T + 4\mathbf{u}\mathbf{u}^T = I \\
H_{\mathbf{u}}^TH_{\mathbf{u}} &= (I - 2\mathbf{u}\mathbf{u}^T)^T(I - 2\mathbf{u}\mathbf{u}^T) = I - 2\mathbf{u}\mathbf{u}^T - 2\mathbf{u}\mathbf{u}^T + 4\mathbf{u}\mathbf{u}^T = I
\end{aligned}
\end{equation}
\subsubsection{Бидиагонализация матрицы}
Бидиагонализация с помощью преобразований Хаусхолдера -- итерационный процесс с наперед заданным количеством шагов. Если матрица имеет размерность $m \times n$, то для бидиагонализации требуется $\min{(m,n)}$ шагов.
Покажем на примере итерационный процесс сведения матрицы $A_{4 \times 4}$ к бидиагональному виду:
\begin{equation} \nonumber
\begin{bmatrix}
\bullet & * & * & * \\
\bullet & * & * & * \\
\bullet & * & * & * \\
\bullet & * & * & *
\end{bmatrix} \rightarrow \begin{bmatrix}
\clubsuit & \bullet & \bullet & \bullet \\
& * & * & * \\
& * & * & * \\
& * & * & *
\end{bmatrix} \rightarrow \begin{bmatrix}
\clubsuit & \clubsuit & & \\
& \bullet & * & * \\
& \bullet & * & * \\
& \bullet & * & *
\end{bmatrix} \rightarrow ... \rightarrow \begin{bmatrix}
\clubsuit & \clubsuit & & \\
& \clubsuit & \clubsuit & \\
& & \clubsuit & \clubsuit \\
& & & \clubsuit
\end{bmatrix}
\end{equation}
В данном примере можно выделить 2 визуально разных преобразования:
\begin{enumerate}
\item обнуление элементов столбца, стоящих под главной диагональю;
\item обнуление элементов строки, стоящих правее бидиагонали.
\end{enumerate}
На практике эти преобразования аналогичные, т.е. если произвести транспонирование, то обнуление элементов строки матрицы $A$, стоящих правее бидиагонали, становится эквивалентным обнулению элементов стобца, стоящих под главной диагональю подматрицы матрицы $A^T$.
Представим матрицу
\begin{equation} \nonumber
A = \left[ \begin{array}{cccc}
\clubsuit & \bullet & \bullet & \bullet \\
& * & * & * \\
& * & * & * \\
& * & * & *
\end{array} \right]
\end{equation}
в виде блочной матрицы вида
\begin{equation} \nonumber
A = \left[ \begin{array}{c|ccc}
\clubsuit & \bullet & \bullet & \bullet \\
& * & * & * \\
& * & * & * \\
\hline
& * & * & *
\end{array} \right] = \left[ \begin{array}{c|ccc}
\clubsuit & \multicolumn{3}{c}{\multirow{3}{*}{\raisebox{-1mm}{\scalebox{1.5}{$\widetilde{A}$}}}} \\
& & & \\
& & & \\
\hline
& * & * & *
\end{array} \right],
\end{equation}
где матрица $\widetilde{A}$ записывается в виде:
\begin{equation} \nonumber
\widetilde{A} = \left[ \begin{array}{ccc}
\bullet & \bullet & \bullet \\
* & * & * \\
* & * & *
\end{array} \right].
\end{equation}
Тогда, транспонировав матрицу $\widetilde{A}$
\begin{equation} \nonumber
\widetilde{A}^T = \left[ \begin{array}{ccc}
\bullet & * & * \\
\bullet & * & * \\
\bullet & * & *
\end{array} \right],
\end{equation}
мы сведем задачу обнуления элементов строки матрицы $A$, стоящих правее бидиагонали, к задаче обнуления элементов столбца, стоящих под главной диагональю матрицы $\widetilde{A}^T$.
Покажем, как с помощью преобразований Хаусхолдера обнулить элементы столбца, стоящие под главной диагональю. Рассмотрим для простоты матрицу
\begin{equation} \nonumber
A_{2 \times 2} = \left[ \begin{array}{cc}
x_{1}^{(1)} & x_{1}^{(2)} \\
x_{2}^{(1)} & x_{2}^{(2)}
\end{array} \right],
\end{equation}
столбцы которой являются радиус-векторами в пространстве $\mathbb{R}^2$. Необходимо найти такое преобразование Хаусхолдера $H_{\mathbf{u}}$, которое бы выполняло следующее действие:
\begin{equation} \nonumber
H_{\mathbf{u}}\left[ \begin{array}{c}
x_{1}^{(1)} \\
x_{2}^{(1)}
\end{array} \right] = \left[ \begin{array}{c}
\widetilde{x}_{1}^{(1)} \\
0
\end{array} \right]
\end{equation}
\begin{wrapfigure}{l}{0pt}
\begin{tikzpicture}[scale=4]
\draw[<->] (0,1.5) -- (0,0) -- (1.5,0);
\draw[thick,->] (0,0) -- (0.638,1.183) node[anchor=base,black] at (0.27,0.73) {$x^{(1)}$};
\draw[red,thick,dashed] (-0.1000,-0.0597) -- (1.5000,0.8956) node[anchor=base] at (1.4,0.7) {$L$};
\draw[blue,thick,->] (0,0) -- (-0.1538,0.2576) node[anchor=base] at (-0.2,0.29) {$\mathbf{u}$};
\draw[green,thick,->] (0,0) -- (1.3441,0) node[anchor=base,black] at (0.71,0.07) {$\widetilde{x}^{(1)}=H_{\mathbf{u}}x^{(1)}$};
\draw[blue,thick,->] (1.3441,0) -- (0.638,1.183) node[anchor=base] at (1.15,1) {$\mathbf{v} = x^{(1)}-\widetilde{x}^{(1)}$};
\end{tikzpicture}
\hspace{20pt}
\caption{Обнуление $x_2$-координаты}
\label{fig:bidiagonalization}
\end{wrapfigure}
Зная исходный вектор и вид получаемого вектора, можно найти вектор $\mathbf{u}$, задающий матрицу преобразования Хаусхолдера. На рис.\ref{fig:bidiagonalization} представлена геометрическая интерпретация требуемого преобразования. Т.к. преобразование Хаусхолдера является изометрическим, то \[ ||x^{(1)}|| = ||\widetilde{x}^{(1)}||. \]
Это условие можно переписать в виде \[ H_{\mathbf{u}}x^{(1)} = \pm ||x^{(1)}|| e_1. \] Для того, чтобы найти $\mathbf{u}$, необходимо выразить $\mathbf{v}$, причем необходимо учитывать тот факт, что норма вектора $\mathbf{v}$ стоит в знаменателе, а значит, надо выбрать такое представление, чтобы избежать деления на близкое к нулю значение: \[ \mathbf{v} = x^{(1)}-\widetilde{x}^{(1)} = x^{(1)} - \pm ||x^{(1)}||e_1 = x^{(1)} \mp ||x^{(1)}||e_1 = x^{(1)} + \sign{(x^{(1)}_1)}||x^{(1)}||e_1. \]
Учитывая, что $\mathbf{u} = \dfrac{\mathbf{v}}{||\mathbf{v}||}$, получаем выражение: \[ \mathbf{u} = \dfrac{x^{(1)} + \sign{(x^{(1)}_1)}||x^{(1)}||e_1}{||x^{(1)} + \sign{(x^{(1)}_1)}||x^{(1)}||e_1||}. \]
Описание алгоритма: пусть дана матрица $M_{m \times n}$, которую необходимо привести к бидиагональному виду. \[ M = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & *
} \]
На каждом шаге обнуляются элементы одного столбца и одной строки. Рассмотрим первые 2 полных шага:
Пусть $A^{(1.1)} = M$. Обнулим элементы первого столбца матрицы \[A^{(1.1)} = \kbordermatrix{
\noindex & a_1 & & \noindex & \noindex & A_2 & \noindex & \noindex \\
\noindex & \bullet & \vrule & * & * & * & \hdots & * \\
\noindex & \bullet & \vrule & * & * & * & \hdots & * \\
\noindex & \bullet & \vrule & * & * & * & \hdots & * \\
\noindex & \vdots & \vrule & \vdots & \vdots & \vdots & & \vdots \\
\noindex & \bullet & \vrule & * & * & * & \hdots & * \\
\noindex & \bullet & \vrule & * & * & * & \hdots & *
} \rightarrow \kbordermatrix{
\noindex & \widetilde{a}_1 & & \noindex & \noindex & \widetilde{A}_2 & \noindex & \noindex \\
\noindex & \clubsuit & \vrule & \bullet & \bullet & \bullet & \hdots & \bullet \\
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\noindex & \vdots & \vrule & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\cline{2-8}
\noindex & 0 & \vrule & * & * & * & \hdots & *
} = H^1_1A^{(1.1)} = \widetilde{A}^{(1.1)} \]
Процесс конструирования матрицы преобразования Хаусхолдера был только что описан, поэтому воспользуемся готовой формулой. Учитывая, что $x^{(1)} := a_1$, найдем вектор $\mathbf{u}$, однозначно задающий матрицу преобразования:\[ \mathbf{u} = \dfrac{a_1 + \sign{(a_{11})}||a_1||e_1}{||a_1 + \sign{(a_{11})}||a_1||e_1||}; \quad H^{1}_{1} = I-2\mathbf{u}\mathbf{u}^T. \]
Т.о. применяя данное преобразование к матрице $A^{(1.1)}$, обнуляются элементы первого столбца под главной диагональю.
Сведем задачу обнуления элементов первой строки, стоящих правее бидиагонали, к задаче обнуления элементов первой строки, стоящих правее главной диагонали. Для этого достаточно рассмотреть подматрицу $\widetilde{A}_2$ матрицы $\widetilde{A}^{(1.1)}$. (Причем, если размерность матрицы $\widetilde{A}^{(1.1)}$ -- $m \times n$, то размерность подматрицы $\widetilde{A}_2$ -- $(m-1)\times(n-1)$).
\[ \widetilde{A}_2 = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \bullet & \bullet & \bullet & \hdots & \bullet \\
\cline{2-6}
\noindex & * & * & * & \hdots & * \\
\noindex & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & \hdots & *
} \]
Сведем задачу обнуления элементов строки матрицы $\widetilde{A}_2$, стоящих правее главной диагонали к задаче обнуления элементов столбца, стоящих под главной диагональю матрицы $A^{(1.2)} = \widetilde{A}^T_2$. \[ A^{(1.2)} = \kbordermatrix{
\noindex & a_1 & \noindex & \noindex & \noindex & A_2 & \noindex \\
\noindex & \bullet & \vrule & * & * & \hdots & * \\
\noindex & \bullet & \vrule & * & * & \hdots & * \\
\noindex & \bullet & \vrule & * & * & \hdots & * \\
\noindex & \vdots & \vrule & \vdots & \vdots & & \vdots \\
\noindex & \bullet & \vrule & * & * & \hdots & *
} \rightarrow \kbordermatrix{
\noindex & \widetilde{a}_1 & \noindex & \noindex & \noindex & \widetilde{A}_2 & \noindex \\
\noindex & \clubsuit & \vrule & \bullet & \bullet & \hdots & \bullet \\
\noindex & 0 & \vrule & * & * & \hdots & * \\
\noindex & 0 & \vrule & * & * & \hdots & * \\
\noindex & \vdots & \vrule & \vdots & \vdots & & \vdots \\
\noindex & 0 & \vrule & * & * & \hdots & *
} = H^1_2A^{(1.2)} = \widetilde{A}^{(1.2)} \]
Транспонировав обратно матрицу $\widetilde{A}^{(1.2)}$, мы получим матрицу с обнуленными элементами правее главной диагонали: \[ (\widetilde{A}^{(1.2)})^T = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & 0 & 0 & \hdots & 0 \\
\cline{2-6}
\noindex & \bullet & * & * & \hdots & * \\
\noindex & \bullet & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & & \vdots \\
\noindex & \bullet & * & * & \hdots & *
} \]
\begin{equation} \nonumber
\left.
\begin{aligned}
A^{(1.2)} = \widetilde{A}^T_2 \\
H^1_2A^{(1.2)} = \widetilde{A}^{(1.2)}
\end{aligned}
\right\} \Rightarrow \left. \begin{aligned}
(A^{(1.2)})^T = \widetilde{A}_2 \\
(A^{(1.2)})^T(H^1_2)^T = (\widetilde{A}^{(1.2)})^T
\end{aligned}
\right\} \Rightarrow \widetilde{A}_2 (H^1_2)^T = (\widetilde{A}^{(1.2)})^T
\end{equation}
Составим такие матрицы $H^{(1.1)}, H^{(1.2)}$, чтобы выполнялось следующее равенство: \[
H^{(1.1)} \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & *
} (H^{(1.2)})^T = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
}
\]
Мы нашли матрицу $H^{(1.1)} := H^1_1$, которая обнуляет элементы первого столбца, стоящие под главной диагональю: \[
H^{(1.1)} : \underbrace{\kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & *
}}_{m \times n} \rightarrow \underbrace{\kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
}}_{m \times n},
\]
а также матрицу $H^1_2$, которая обнуляет элементы первой строки, стоящих правее бидиагонали матрицы, полученной после применения преобразования $H^{(1.1)}$: \[
H^1_2 : \underbrace{\kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & * & * & * & \hdots & * \\
\noindex & * & * & * & \hdots & * \\
\noindex & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & \hdots & *
}}_{(m-1) \times (n-1)} \rightarrow \underbrace{\kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & * & * & * & \hdots & * \\
\noindex & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & \hdots & *
}}_{(m-1) \times (n-1)}
\]
Представим матрицу $H^{(1.1)}A^{(1.1)} = \widetilde{A}^{(1.1)}$ в виде блочной матрицы \[ \widetilde{A}^{(1.1)} = \left[ \begin{array}{c|c}
B_{(m-1)\times 1} & A^{(1.2)}_{(m-1)\times (n-1)} \\
\cline{1-2}
O_{1\times 1} & *_{1\times (n-1)}
\end{array} \right],
\]
для которой надо составить такую матрицу $H^{(1.2)}$, чтобы
\[ \left[ \begin{array}{c|c}
B_{(m-1)\times 1} & A^{(1.2)}_{(m-1)\times (n-1)} \\
\cline{1-2}
O_{1\times 1} & *_{1\times (n-1)}
\end{array} \right] (H^{(1.2)})^T = \left[ \begin{array}{c|c}
B_{(m-1)\times 1} & \widetilde{A}^{(1.2)}_{(m-1)\times (n-1)} \\
\cline{1-2}
O_{1\times 1} & *_{1\times (n-1)}
\end{array} \right]. \]
Подходящей матрицей является матрица вида \[
H^{(1.2)} = \left[ \begin{array}{c|c}
I_{1\times 1} & O_{1\times (n-1)} \\
\cline{1-2}
O_{(n-1)\times 1} & (H^1_2)_{(n-1)\times (n-1)}
\end{array} \right]
\]
Покажем, что $H^{(1.2)}$ -- именно та матрица, которую мы ищем: \[
(H^{(1.2)})^T = \left[ \begin{array}{c|c}
I_{1\times 1} & O_{1\times (n-1)} \\
\cline{1-2}
O_{(n-1)\times 1} & (H^1_2)_{(n-1)\times (n-1)}
\end{array} \right]^T = \left[ \begin{array}{c|c}
I_{1\times 1} & O_{1\times (n-1)} \\
\cline{1-2}
O_{(n-1)\times 1} & (H^1_2)^T_{(n-1)\times (n-1)}
\end{array} \right]
\]\[ \left[ \begin{array}{c|c}
B_{(m-1)\times 1} & A^{(1.2)}_{(m-1)\times (n-1)} \\
\cline{1-2}
O_{1\times 1} & *_{1\times (n-1)}
\end{array} \right] \left[ \begin{array}{c|c}
I_{1\times 1} & O_{1\times (n-1)} \\
\cline{1-2}
O_{(n-1)\times 1} & (H^1_2)^T_{(n-1)\times (n-1)}
\end{array} \right] = \left[ \begin{array}{c|c}
B_{(m-1)\times 1} & (\widetilde{A}^{(1.2)})^T_{(m-1)\times (n-1)} \\
\cline{1-2}
O_{1\times 1} & *_{1\times (n-1)}
\end{array} \right]. \]
\begin{equation} \nonumber
\left\{
\begin{aligned}
& B_{(m-1)\times 1} \cdot I_{1\times 1} + A^{(1.2)}_{(m-1)\times (n-1)} \cdot O_{(n-1)\times 1} = B_{(m-1)\times 1} \\
& B_{(m-1)\times 1} \cdot O_{1 \times (n-1)} + A^{(1.2)}_{(m-1)\times (n-1)} \cdot (H^1_2)^T_{(n-1)\times (n-1)} = (\widetilde{A}^{(1.2)})^T_{(m-1)\times (n-1)} \\
& O_{1 \times 1} \cdot I_{1 \times 1} + *_{1 \times (n-1)} \cdot O_{(n-1)\times 1} = O_{1\times 1} \\
& O_{1 \times 1} \cdot O_{1 \times (n-1)} + *_{1 \times (n-1)} \cdot (H^1_2)^T_{(n-1)\times (n-1)} = *_{1 \times (n-1)}
\end{aligned}
\right.
\end{equation}
Т.о. матрица, получающаяся в результате преобразования $H^{(1)}, H^{(2)}$ матрица, удовлетворяет требованиям: в первом столбце элементы, стоящие под главной диагональю, равны нулю; в первой строке элементы, стоящие правее бидиагонали, также равны нулю.
\[
H^{(1.1)}A^{(1.1)}(H^{(1.2)})^T = H^{(1.1)} \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & * & * & * & * & \hdots & * \\
\noindex & * & * & * & * & \hdots & *
}(H^{(1.2)})^T = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
} = A^{(2.1)}
\]
Первый шаг алгоритма полностью выполнен. На втором шаге алгоритм будет работать уже с матрицей $A^{(2.1)}$, для того, чтобы привести ее к виду \[
A^{(2.1)} = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
} \rightarrow \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & \clubsuit & \clubsuit & 0 & \hdots & 0 \\
\noindex & 0 & 0 & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & 0 & * & * & \hdots & * \\
\noindex & 0 & 0 & * & * & \hdots & *
} = A^{(3.1)}
\]
Матрицу $A^{(2.1)}$ можно представить в виде блочной матрицы
\[
A^{(2.1)} = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
} = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \vrule & \clubsuit & 0 & 0 & \hdots & 0 \\
\cline{2-8}
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\noindex & \vdots & \vrule & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & \vrule & * & * & * & \hdots & * \\
\noindex & 0 & \vrule & * & * & * & \hdots & *
} = \left[ \begin{array}{c|ccccc}
\clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\cline{1-6}
0 & \multicolumn{5}{c}{\multirow{5}{*}{\raisebox{-1mm}{\scalebox{1.8}{$M_2$}}}} \\
0 & & & & & \\
\vdots & & & & & \\
0 & & & & & \\
0 & & & & &
\end{array} \right]
\]
\[
A^{(2.1)} = \left[ \begin{array}{c|c}
D_{1\times 1} & R_{1\times (n-1)} \\
\cline{1-2}
O_{(m-1)\times 1} & (M_2)_{(m-1)\times (n-1)}
\end{array} \right]
\]
Повторяя первый шаг, находим матрицы $H^2_1, H^2_2$ для матрицы $M_2$. Однако искомые матрицы $H^{(2.1)}, H^{(2.2)}$ таковы, что \[
H^{(2.1)} \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & * & * & * & \hdots & * \\
\noindex & 0 & * & * & * & \hdots & *
} (H^{(2.2)})^T = \kbordermatrix{
\noindex & \noindex & \noindex & \noindex & \noindex & \noindex & \noindex \\
\noindex & \clubsuit & \clubsuit & 0 & 0 & \hdots & 0 \\
\noindex & 0 & \clubsuit & \clubsuit & 0 & \hdots & 0 \\
\noindex & 0 & 0 & * & * & \hdots & * \\
\noindex & \vdots & \vdots & \vdots & \vdots & & \vdots \\
\noindex & 0 & 0 & * & * & \hdots & * \\
\noindex & 0 & 0 & * & * & \hdots & *
}
\]
Т.е. они должны быть сконструированы таким образом, чтобы произвести необходимое преобразование подматрицы $M_2$, при этом не изменив $D$ и $R$. Представим матрицы $H^{(2.1)}, H^{(2.2)}$ в виде блочных матриц:
\[
H^{(2.1)} = \underbrace{\left[ \begin{array}{c|c}
H^{(2.1)}_1 & H^{(2.1)}_2 \\
\cline{1-2}
H^{(2.1)}_3 & H^{(2.1)}_4
\end{array} \right]}_{m \times m} \quad H^{(2.2)} = \underbrace{\left[ \begin{array}{c|c}
H^{(2.2)}_1 & H^{(2.2)}_2 \\
\cline{1-2}
H^{(2.2)}_3 & H^{(2.2)}_4
\end{array} \right]}_{n \times n}
\]
Найдем сперва элементы матрицы $H^{(2.1)}$:
\[
\left[ \begin{array}{c|c}
H^{(2.1)}_1 & H^{(2.1)}_2 \\
\cline{1-2}
H^{(2.1)}_3 & H^{(2.1)}_4
\end{array} \right] \left[ \begin{array}{c|c}
D_{1\times 1} & R_{1\times (n-1)} \\
\cline{1-2}
O_{(m-1)\times 1} & (M_2)_{(m-1)\times (n-1)}
\end{array} \right] = \left[ \begin{array}{c|c}
D_{1\times 1} & R_{1\times (n-1)} \\
\cline{1-2}
O_{(m-1)\times 1} & (\widetilde{M}_2)_{(m-1)\times (n-1)}
\end{array} \right]
\]
Составим систему уравнений, откуда найдем элементы блочной матрицы $H^{(2.1)}$ и их размерности.
\begin{equation} \nonumber
\left\{
\begin{aligned}
& H^{(2.1)}_1 \cdot D_{1 \times 1} + H^{(2.1)}_2 \cdot O_{(m-1) \times 1} = D_{1 \times 1} \\
& H^{(2.1)}_1 \cdot R_{1 \times (n-1)} + H^{(2.1)}_2 \cdot (M_2)_{(m-1)\times (n-1)} = R_{1 \times (n-1)} \\
& H^{(2.1)}_3 \cdot D_{1 \times 1} + H^{(2.1)}_4 \cdot O_{(m-1) \times 1} = O_{(m-1) \times 1} \\
& H^{(2.1)}_3 \cdot R_{1 \times (n-1)} + H^{(2.1)}_4 \cdot (M_2)_{(m-1) \times (n-1)} = (\widetilde{M}_2)_{(m-1)\times (n-1)}
\end{aligned}
\right. \Rightarrow \left\{
\begin{aligned}
& H^{(2.1)}_1 = I_{1 \times 1} \\
& H^{(2.1)}_2 = O_{1 \times (m-1)} \\
& H^{(2.1)}_3 = O_{(m-1) \times 1} \\
& H^{(2.1)}_4 = (H^2_1)_{(m-1) \times (m-1)}
\end{aligned}
\right.
\end{equation}
\[
H^{(2.1)} = \left[ \begin{array}{c|c}
I_{1 \times 1} & O_{1 \times (m-1)} \\
\cline{1-2}
O_{(m-1) \times 1} & (H^2_1)_{(m-1) \times (m-1)}
\end{array} \right] = \left[ \begin{array}{c|c}
I_1 & O \\
\cline{1-2}
O & H^2_1
\end{array} \right]
\]
Найдем элементы матрицы $H^{(2.2)}$ (необходимо помнить, что матрица транспонируется):
\[
\left[ \begin{array}{c|c}
D_{1\times 1} & R_{1\times (n-1)} \\
\cline{1-2}
O_{(m-1)\times 1} & (\widetilde{M}_2)_{(m-1)\times (n-1)}
\end{array} \right] \left[ \begin{array}{c|c}
H^{(2.2)}_1 & H^{(2.2)}_3 \\
\cline{1-2}
H^{(2.2)}_2 & H^{(2.2)}_4
\end{array} \right] = \left[ \begin{array}{c|c}
D_{1\times 1} & R_{1\times (n-1)} \\
\cline{1-2}
O_{(m-1)\times 1} & (\widetilde{M}_2)_{(m-1)\times (n-1)}
\end{array} \right]
\]
Составим вторую систему уравнений, откуда найдем элементы блочной матрицы $H^{(2.2)}$ и их размерности.
\begin{equation} \nonumber
\left\{
\begin{aligned}
& D_{1 \times 1} \cdot H^{(2.2)}_1 + R_{1 \times (n-1)} \cdot H^{(2.2)}_2 = D_{1 \times 1} \\
& D_{1 \times 1} \cdot H^{(2.2)}_3 + R_{1 \times (n-1)} \cdot H^{(2.2)}_4 = R_{1 \times (n-1)} \\
& O_{(m-1) \times 1} \cdot H^{(2.2)}_1 + (\widetilde{M}_2)_{(m-1)\times (n-1)} \cdot H^{(2.2)}_2 = O_{(m-1)\times 1} \\
& O_{(m-1) \times 1} \cdot H^{(2.2)}_3 + (\widetilde{M}_2)_{(m-1)\times (n-1)} \cdot H^{(2.2)}_4 = (\widetilde{A}^{(2.1)})_{(m-1)\times (n-1)}
\end{aligned}
\right. \Rightarrow \left\{
\begin{aligned}
& H^{(2.1)}_1 = I_{1 \times 1} \\
& H^{(2.1)}_2 = O_{(n-1) \times 1} \\
& H^{(2.1)}_3 = O_{1 \times (n-1)} \\
& H^{(2.1)}_4 = (H^2_1)_{(n-1) \times (n-1)}
\end{aligned}
\right.
\end{equation}
В данном случае следует отдельно разобрать второе уравнение, в котором утверждается, что
\[
R_{1 \times (n-1)} \cdot H^{(2.2)}_4 = R_{1 \times (n-1)}.
\]
Когда мы искали вид матрицы $H^{(1.2)}$, было показано, что она имеет вид:
\[
H^{(1.2)} = \left[ \begin{array}{c|c}
1 & O_{1\times (n-1)} \\
\cline{1-2}
O_{(n-1)\times 1} & (H^1_2)_{(n-1)\times (n-1)}
\end{array} \right].
\]
Подматрица $R_{1 \times (n-1)}$ имеет вид:
\[
R_{1 \times (n-1)} = \left[
\begin{array}{ccccc}
\clubsuit & 0 & 0 & \hdots & 0
\end{array}
\right],
\]
т.е. $R_{1 \times (n-1)}$ инвариантна относительно преобразования $H^{(2.2)}_4$.
Таким образом, матрица $H^{(2.2)}$ имеет вид:
\[
H^{(2.2)} = \left[ \begin{array}{c|c}
I_{1 \times 1} & O_{1 \times (n-1)} \\
\cline{1-2}
O_{(n-1) \times 1} & (H^2_1)_{(n-1) \times (n-1)}
\end{array} \right] = \left[ \begin{array}{c|c}
I_1 & O \\
\cline{1-2}
O & H^2_1
\end{array} \right]
\]
Только что было показано, что исходная матрица $M$ соотносится с матрицей, полученной после двух шагов $A^{(3.1)}$ следующим образом:
\[
H^{(2.1)}H^{(1.1)} \cdot M \cdot (H^{(1.1)})^T(H^{(2.1)})^T = (H^{(2.1)}H^{(1.1)}) \cdot M \cdot (H^{(2.1)}H^{(1.1)})^T = A^{(3.1)}
\]
Продолжая аналогичные рассуждения, можно написать формулы нахождения матриц $H^{(k.1)},H^{(k.2)}$ на $k$-ом шаге:
\begin{equation} \nonumber
\begin{aligned}
H^{(k.1)} = \left[ \begin{array}{c|c}
I_{(k-1) \times (k-1)} & O_{(n-k) \times (k-1)} \\
\cline{1-2}
O_{(k-1) \times (n-k)} & (H^k_1)_{(n-k) \times (n-k)}
\end{array} \right] = \left[ \begin{array}{c|c}
I_{k-1} & O \\
\cline{1-2}
O & H^k_1
\end{array} \right] \\
H^{(k.2)} = \left[ \begin{array}{c|c}
I_{(k-1) \times (k-1)} & O_{(k-1) \times (n-k)} \\
\cline{1-2}
O_{(n-k) \times (k-1)} & (H^k_2)_{(n-k) \times (n-k)}
\end{array} \right] = \left[ \begin{array}{c|c}
I_{k-1} & O \\
\cline{1-2}
O & H^k_2
\end{array} \right]
\end{aligned}
\end{equation}
Тогда матрица $A^{((k+1).1)}$ находится по формуле:
\[
A^{((k+1).1)} = \Bigg( H^{(k.1)} \cdot ... \cdot H^{(1.1)} \Bigg) \cdot M \cdot \Bigg( H^{(k.2)} \cdot ... \cdot H^{(1.2)} \Bigg)^T.
\]
Рассмотрим бидиагонализацию трех возможных видов матриц. Следует отдельно отметить тот факт, что ранг матрицы не влияет на работу алгоритма, т.е. вне зависимости от того $\rank{(A)} < \min{(m,n)}$ или $\rank{(A)} = \min{(m,n)}$, количество шагов работы алгоритма одинаково в обоих случаях (алгоритм инвариантен относительно ранга матрицы).
Первым делом рассмотрим самый простой случай, когда $m=n$, т.е. матрица квадратная.
\begin{wrapfigure}{l}{0pt}
\begin{tikzpicture}
\matrix [matrix of math nodes,left delimiter={[},right delimiter={]}] (m) {
* & \hdots & * \\
\vdots & \ddots & \vdots \\
* & \hdots & * \\};
\node[above=1pt of m-1-1] (top-1) {};
\node[above=1pt of m-1-3] (top-2) {};
\node[left=10pt of m-1-1] (left-1) {};
\node[left=10pt of m-3-3] (left-2) {};
\node[rectangle,above delimiter=\{] (del-top-1) at ($0.5*(top-1.south) +0.5*(top-2.south)$) {\tikz{\path (top-1.south west) rectangle (top-2.north east);}};
\node[above=10pt] at (del-top-1.north) {$n$ cols};
\node[rectangle,left delimiter=\{] (del-left-1) at ($0.5*(left-1.east) +0.5*(left-2.east)$) {\tikz{\path (left-1.north east) rectangle (left-2.south west);}};
\node[left=10pt] at (del-left-1.west) {$n$ rows};
\end{tikzpicture}
\end{wrapfigure}
На каждом шаге алгоритм работает с одним столбцом и одной строкой. У квадратной матрицы $n \times n$ необходимо обработать $(n-1)$ столбец и $(n-2)$ строки, следовательно бидиагональная матрица $B$ получается из исходной матрицы $A$ следующим образом: \[ B = \Bigg( H^{((n-1).1)} \cdot ... \cdot H^{(1.1)} \Bigg) \cdot A \cdot \Bigg( H^{((n-2).2)} \cdot ... \cdot H^{(1.2)} \Bigg)^T. \] Воспользовавшись условием ортонормированности матриц преобразований Хаусхолдера, выпишем разложение матрицы $A$ в виде произведения двух ортонормированных и одной бидиагональной матрицы:
\begin{equation} \nonumber
\begin{aligned}
A &= UBV^T, \\
U = \Big( H^{((n-1).1)} \cdot ... \cdot H^{(1.1)} \Big)^T; &\quad V = \Big( H^{((n-2).2)} \cdot ... \cdot H^{(1.2)} \Big)^T
\end{aligned}
\end{equation}
\begin{wrapfigure}{l}{0pt}
\begin{tikzpicture}
\matrix [matrix of math nodes,left delimiter={[},right delimiter={]}] (m) {
* & \hdots & * & * & \hdots & * \\
\vdots & \ddots & \vdots & * & \hdots & * \\
* & \hdots & * & * & \hdots & * \\};
\node[above=1pt of m-1-1] (top-1) {};
\node[above=1pt of m-1-3] (top-2) {};
\node[above=1pt of m-1-4] (top-3) {};
\node[above=1pt of m-1-6] (top-4) {};
\node[left=10pt of m-1-1] (left-1) {};
\node[left=10pt of m-3-3] (left-2) {};
\draw[dashed] ($0.5*(m-1-3.north east)+0.5*(m-1-4.north west)$) --
($0.5*(m-3-3.south east)+0.5*(m-3-4.south west)$);
\node[rectangle,above delimiter=\{] (del-top-1) at ($0.5*(top-1.south) +0.5*(top-2.south)$) {\tikz{\path (top-1.south east) rectangle (top-2.north west);}};
\node[above=10pt] at (del-top-1.north) {$m$ cols};
\node[rectangle,above delimiter=\{] (del-top-1) at ($0.5*(top-3.south) +0.5*(top-4.south)$) {\tikz{\path (top-3.south east) rectangle (top-4.north west);}};
\node[above=10pt] at (del-top-1.north) {$m-n$ cols};
\node[rectangle,left delimiter=\{] (del-left-1) at ($0.5*(left-1.east) +0.5*(left-2.east)$) {\tikz{\path (left-1.north east) rectangle (left-2.south west);}};
\node[left=10pt] at (del-left-1.west) {$m$ rows};
\end{tikzpicture}
\end{wrapfigure}
Рассмотрим случай, когда $m < n$. Из графического представления матрицы, у которой $m < n$ видно, что исходную матрицу $A_{m \times n}$ можно представить в виде блочной матрицы \[ A = \left[ \begin{array}{c|c}
S_{m \times m} & C
\end{array} \right]. \]
Т.к. на каждом шаге алгоритм работает с одним столбцом и одной строкой, а под главной диагональю находится $(m-1)$ столбец, то очевидно, что \[ U = \Big( H^{((m-1).1)} \cdot ... \cdot H^{(1.1)} \Big)^T. \] Строки матрицы $A$ имеют больше элементов по сравнению с квадратной матрицей, поэтому потребуется обрабатывать больше строк. Если подматрица $C$ является столбцом, то обрабатывается на одну строку больше по сравнению с квадратной матрицей, если $C$ -- матрица, имеющая $\ge 2$ количество столбцов, то обрабатываются все $m$ строк, и матрица $V$ принимает вид: \[ V = \Big( H^{(m.2)} \cdot ... \cdot H^{(1.2)} \Big)^T. \]
\begin{wrapfigure}{l}{0pt}
\begin{tikzpicture}
\matrix [matrix of math nodes,left delimiter={[},right delimiter={]}] (m) {
* & \hdots & * \\
\vdots & \ddots & \vdots \\
* & \hdots & * \\
* & \hdots & * \\
\vdots & & \vdots \\
* & \hdots & * \\};
\node[above=1pt of m-1-1] (top-1) {};
\node[above=1pt of m-1-3] (top-2) {};
\node[left=10pt of m-1-1] (left-1) {};
\node[left=10pt of m-3-1] (left-2) {};
\node[left=10pt of m-5-1] (left-3) {};
\node[left=10pt of m-6-1] (left-4) {};
\draw[dashed] ($0.5*(m-3-1.south west)+0.5*(m-4-1.north west)$) --
($0.5*(m-3-3.south east)+0.5*(m-4-3.north east)$);
\node[rectangle,above delimiter=\{] (del-top-1) at ($0.5*(top-1.south) +0.5*(top-2.south)$) {\tikz{\path (top-1.south east) rectangle (top-2.north west);}};
\node[above=10pt] at (del-top-1.north) {$n$ cols};
\node[rectangle,left delimiter=\{] (del-left-1) at ($0.5*(left-1.east) +0.5*(left-2.east)$) {\tikz{\path (left-1.north east) rectangle (left-2.north west);}};
\node[left=10pt] at (del-left-1.west) {$n$ rows};
\node[rectangle,left delimiter=\{] (del-left-2) at ($0.5*(left-3.north) +0.5*(left-4.north)$) {\tikz{\path (left-1.north east) rectangle (left-2.north west);}};
\node[left=10pt] at (del-left-2.west) {$m-n$ rows};
\end{tikzpicture}
\end{wrapfigure}
Последний случай, когда $m > n$ похож на предыдущий, за тем исключением, что теперь у столбцов больше элементов, которые необходимо обнулить. В данном случае бидиагональная матрица $B$ получается из исходной матрицы $A$ следующим образом: \[ B = \Bigg( H^{(n.1)} \cdot ... \cdot H^{(1.1)} \Bigg) \cdot A \cdot \Bigg( H^{((n-2).2)} \cdot ... \cdot H^{(1.2)} \Bigg)^T. \] Главное отличие от случая с квадратной матрицей -- обрабатывается последний столбец. Выпишем отдельно матрицы $U,V$:
\begin{equation} \nonumber
\begin{aligned}
A &= UBV^T, \\
U = \Big( H^{(n.1)} \cdot ... \cdot H^{(1.1)} \Big)^T; &\quad V = \Big( H^{((n-2).2)} \cdot ... \cdot H^{(1.2)} \Big)^T
\end{aligned}
\end{equation}
\subsubsection{Вращения Гивенса} \label{sec:givens_rot}
Очень важной частью \quotes{стандартного} алгоритма является эффективное вычисление элементов матриц вращений Гивенса. Более подробно описать задачу построения матриц Гивенса можно следующим образом: пусть имеется вектор $(f,g)^T$, необходимо образовать матрицу поворота, которая приведет вектор вида $(f,g)^T$ к виду $(r,0)^T$: \begin{equation} \label{eq:givens_rot}
\begin{bmatrix}
c & s \\ -s & c
\end{bmatrix} \begin{bmatrix}
f \\ g
\end{bmatrix} = \begin{bmatrix}
r \\ 0
\end{bmatrix},\ \ \bigg|\ \begin{aligned}
c &:= \cos{\theta} \\
s &:= \sin{\theta}
\end{aligned}
\end{equation}
Предположим, что $f > g$. В таком случае: \[t := \tan{\theta} = \dfrac{g}{f}.\]
Перепишем матричное уравнение (\ref{eq:givens_rot}) в виде системы алгебраических уравнений: \begin{equation}
\left\{ \begin{aligned}
cf - sg &= r \\
sf + cg &= 0 \Leftrightarrow g = -\dfrac{s}{c}f
\end{aligned} \right.
\end{equation}
Т.о. $r$ вычисляется по формуле: \[r = cf+\dfrac{s^2}{c}f = \dfrac{c^2f+s^2f}{c} = \dfrac{f}{c}.\]
Теперь остается вспомнить лишь тот факт, что $\cos{\theta} = \dfrac{1}{\sqrt{1+\tan^2{\theta}}}$ и подставить в предыдущую формулу:
\[ r = f \sqrt{1+t^2} \]
При условии, что $g > f$ все рассуждения аналогичны за тем исключением, что элементы векторов $(f,g)^T$ и $(r,0)^T$ записываются в обратном порядке: \[ \begin{bmatrix}
f \\ g
\end{bmatrix} \rightarrow \begin{bmatrix}
g \\ f
\end{bmatrix};\ \ \ \ \begin{bmatrix}
r \\ 0
\end{bmatrix} \rightarrow \begin{bmatrix}
0 \\ r
\end{bmatrix} \]
В работе авторов \cite{citation6} более подробно описывается данный алгоритм, а также изучается погрешность и накопление ошибки.
\subsubsection{\quotes{Стандартный} алгоритм} \label{sec:stand_alg}
Сам \quotes{стандартный} алгоритм собирается как конструктор из вышеописанных частей. Сингулярное разложение разбивается на две большие независимые части:
\begin{enumerate}
\item Бидиагонализация матрицы. \\\textit{\underline{Замечание 1:} \quotes{стандартный} алгоритм работает с вертикальными бидиагональными матрицами, т.е. матрицами $B_{m\times n}$, $m \ge n$. В противном случае перед бидиагонализацией исходную матрицу необходимо транспонировать.} \\\textit{\underline{Замечание 2:} возможна работа с горизонтальными матрицами, однако в таком случае необходимо привести матрицу к нижнему бидиагональному виду и соответствующим образом переопределить QR скачки.}
\item Использование QR скачков для диагонализации бидиагональной матрицы.
\end{enumerate}
Предположим, что $J^{(0)}$ -- бидиагональная матрица. Разновидностью QR-алгоритма можно итерационно диагонализировать бидиагональную матрицу $J^{(0)}$: \[J^{(0)} \rightarrow J^{(1)} \rightarrow ... \rightarrow J^{(n)} \rightarrow ... \rightarrow \Sigma,\ \ \ J^{(i+1)} = S^{(i)^T} J^{(i)} T^{(i)}, \] где преобразования $S^{(i)},T^{(i)}$ -- ортонормированные. Матрицы $T^{(i)}$ выбираются таким образом, чтобы последовательность матриц $M^{(i)} = J^{(i)^T}J^{(i)}$ сходилась к диагональной матрице, в то время как матрицы $S^{(i)}$ выбираются такими, чтобы все $J^{(i)}$ были бидиагональными. \\
Введем обозначения:
\[ J \equiv J^{(i)},\ \ \ \ \bar{J} \equiv J^{(i+1)},\ \ \ \ S \equiv S^{(i)},\ \ \ \ T \equiv T^{(i)},\ \ \ \ M \equiv J^TJ,\ \ \ \ \bar{M} \equiv \bar{J}^T\bar{J}. \]
Преобразование $J \rightarrow \bar{J}$ достигается с помощью применения поворотов Гивенса к матрице $J$ с обоих сторон:
\[ \bar{J} = \underbrace{S_nS_{(n-1)}...S_2}_{S}J\underbrace{T_2^TT_3^T...T_n^T}_{T^T}, \] где $S_i,T_i$ -- матрицы поворота, задающиеся по определенному правилу:
\begin{equation} \nonumber
\begin{aligned}
T_2^T &\text{ обнуляет элемент } J_{1,2} \text{; порождает элемент } J_{2,1} \\
S_2 &\text{ обнуляет элемент } J_{2,1} \text{; порождает элементы } J_{1,2}, J_{1,3} \\
T_3^T &\text{ обнуляет элементы } J_{1,3}, J_{2,3} \text{; порождает элемент } J_{3,2} \\
\hdots &\hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \hdots \\
T_n^T &\text{ обнуляет элементы } J_{n-2,n}, J_{n-1,n} \text{; порождает элемент } J_{n,n-1} \\
S_n &\text{ обнуляет элемент } J_{n,n-1}\text{; порождает элемент } J_{n-1,n}
\end{aligned}
\end{equation}
Рассмотрим качественно, что происходит в данных преобразованиях, на примере матрицы $4\times 4$: \[ J = \begin{bmatrix}
x & x & 0 & 0 \\
0 & x & x & 0 \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{bmatrix},\ \ \ \ J T_2^T = \begin{bmatrix}
x & 0 & 0 & 0 \\
+ & x & x & 0 \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{bmatrix},\ \ \ \ S_2 J T_2^T = \begin{bmatrix}
x & x & + & 0 \\
0 & x & x & 0 \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{bmatrix},
\] \[ S_2 J T_2^T T_3^T = \begin{bmatrix}
x & x & 0 & 0 \\
0 & x & 0 & 0 \\
0 & + & x & x \\
0 & 0 & 0 & x
\end{bmatrix},\ \ \ \ S_3 S_2 J T_2^T T_3^T = \begin{bmatrix}
x & x & 0 & 0 \\
0 & x & x & + \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{bmatrix}, \]
\[ S_3 S_2 J T_2^T T_3^T T_4^T = \begin{bmatrix}
x & x & 0 & 0 \\
0 & x & x & 0 \\
0 & 0 & x & 0 \\
0 & 0 & + & x
\end{bmatrix},\ \ \ \ \bar{J} = S_4 S_3 S_2 J T_2^T T_3^T T_4^T = \begin{bmatrix}
x & x & 0 & 0 \\
0 & x & x & 0 \\
0 & 0 & x & x \\
0 & 0 & 0 & x
\end{bmatrix} \]
В \cite{citation6} доказывается, что вышеприведенные преобразования итеративно приближают матрицу $J$ к диагональному виду: \[ \lim_{n \rightarrow \infty} J^{(n)} = \Sigma. \]
Таким образом набор преобразований, диагонализирующий бидиагональную матрицу $B$, выглядит следующим образом:
\[ \Sigma = \underbrace{ S_n^{(N)}\cdot ... \cdot S_2^{(N)}\cdot ... \cdot S_n^{(1)}\cdot ... \cdot S_2^{(1)}}_{S}B \underbrace{ T_2^{(1)^T}\cdot ... \cdot T_n^{(1)^T}\cdot ... \cdot T_2^{(N)^T}\cdot ... \cdot T_n^{(N)^T}}_{T^T},\ \ N \rightarrow \infty. \]
В (\ref{source_code:qr_sweeps}) приведен исходный код, реализующий одну итерацию QR скачка. Однако в статье \cite{citation6} приведен более эффективный метод реализации скачков, основанный на том, что работа происходит не со всей бидиагональной матрицей, а лишь с двумя ненулевыми диагоналями, которые задаются двумя векторами.
В данной работе рассматривался алгоритм с нулевым сдвигом. Алгоритм \quotes{стандартного} сингулярного разложения с ненулевым сдвигом подробно рассмотрен в работе \cite{citation9}.
Рассмотрим на примере построение матриц поворота $S_i,~ T_i$ с помощью разобранных в разделе (\ref{sec:givens_rot}) вращений Гивенса. Пусть дана исходная верхняя бидиагональная матрица $B_{0_{5\times 4}}$:
\[ B_0 = \begin{bmatrix}
d_1 & u_1 & 0 & 0 \\
0 & d_2 & u_2 & 0 \\
0 & 0 & d_3 & u_3 \\
0 & 0 & 0 & d_4 \\
0 & 0 & 0 & 0
\end{bmatrix},\ \ \ \ \ \begin{aligned}
d_i &\text{ -- $i$-ый диагональный элемент матрицы $B$;} \\
u_i &\text{ -- $i$-ый наддиагональный элемент матрицы $B$;} \\
l_i &\text{ -- $i$-ый поддиагональный элемент матрицы $B$.}
\end{aligned}
\]
Первой создается матрица Гивенса $T_2$. Матрица $T_2$ строится таким образом, чтобы обнулить элемент $u_1$ матрицы $B_0$. Для этого достаточно воспользоваться теорией из раздела (\ref{sec:givens_rot}): \begin{equation} \label{givens_rot_row_to_col}
\begin{bmatrix}
d_1 & u_1
\end{bmatrix} \rightarrow \begin{bmatrix}
d_1^* & 0
\end{bmatrix}. \end{equation}
Покажем, что задача (\ref{givens_rot_row_to_col}) эквивалентна задаче, рассмотренной в (\ref{sec:givens_rot}).
\[ \begin{bmatrix}
d_1 & u_1
\end{bmatrix} \underbrace{\begin{bmatrix}
c & -s \\ s & c
\end{bmatrix}}_{T_2^T} = \Bigg( \bigg( \begin{bmatrix}
d_1 & u_1
\end{bmatrix} \underbrace{\begin{bmatrix}
c & -s \\ s & c
\end{bmatrix}}_{T_2^T} \bigg)^T \Bigg)^T = \bigg( \underbrace{\begin{bmatrix}
c & s \\ -s & c
\end{bmatrix}}_{T_2} \begin{bmatrix}
d_1 \\ u_1
\end{bmatrix} \bigg)^T \rightarrow \begin{bmatrix}
d_1^* \\ 0
\end{bmatrix}^T = \begin{bmatrix}
d_1^* & 0
\end{bmatrix}. \]
Таким образом находим значения $c := \cos{\theta},~ s := \sin{\theta}$ и составляем матрицу поворота Гивенса $T_2^T$: \[ T_2^T = \begin{bmatrix}
c & -s & 0 & 0 \\
s & c & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix},\ \ \ \ B_1 = B_0T_2^T = \begin{bmatrix}
d_1^* & 0 & 0 & 0 \\
l_1 & d_2 & u_2 & 0 \\
0 & 0 & d_3 & u_3 \\
0 & 0 & 0 & d_4 \\
0 & 0 & 0 & 0
\end{bmatrix}. \]
Вторым шагом строится матрица Гивенса $S_2$. Матрица $S_2$ строится таким образом, чтобы обнулить элемент $l_1$ матрицы $B_1$. \[ \underbrace{\begin{bmatrix}
c & s \\ -s & c
\end{bmatrix}}_{S_2} \begin{bmatrix}
d_1^* \\ l_2
\end{bmatrix} \rightarrow \begin{bmatrix}
d_1^{**} \\ 0
\end{bmatrix}. \]
Составляем матрицу поворота Гивенса $S_2$: \[ S_2 = \begin{bmatrix}
c & s & 0 & 0 \\
-s & c & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix},\ \ \ \ B_2 = S_2B_1T_2^T = \begin{bmatrix}
d_1^{**} & u_1^* & t_1 & 0 \\
0 & d_2 & u_2 & 0 \\
0 & 0 & d_3 & u_3 \\
0 & 0 & 0 & d_4 \\
0 & 0 & 0 & 0
\end{bmatrix}. \]
На третьем шаге составляется матрица Гивенса $T_3^T$, которая обнуляет элемент $u_2$: \[ \begin{bmatrix}
d_2 & u_2
\end{bmatrix} \underbrace{\begin{bmatrix}
c & -s \\ s & c
\end{bmatrix}}_{T_3^T} \rightarrow \begin{bmatrix}
d_2^* & 0
\end{bmatrix}. \]
Можно показать, что помимо обнуления $u_2$, обнуляется также элемент $t_1$. Составляем матрицу поворота Гивенса $T_3^T$: \[ T_3^T = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & c & -s & 0 \\
0 & s & c & 0 \\
0 & 0 & 0 & 1
\end{bmatrix},\ \ \ \ B_3 = S_2B_2T_2^TT_3^T = \begin{bmatrix}
d_1^{**} & u_1^* & 0 & 0 \\
0 & d_2 & 0 & 0 \\
0 & l_2 & d_3 & u_3 \\
0 & 0 & 0 & d_4 \\
0 & 0 & 0 & 0
\end{bmatrix}. \]
Аналогично строятся последующие матрицы $S_i,~T_i$.
\subsection{Исходный код}
\subsubsection{Преобразование Хаусхолдера}
\lstinputlisting{source_code/householder.m}
\subsubsection{Бидиагонализация матрицы}
\lstinputlisting{source_code/bidiag_reduction.m}
\subsubsection{Вращения Гивенса}
\lstinputlisting{source_code/rot.m}
\subsubsection{QR скачки} \label{source_code:qr_sweeps}
\lstinputlisting{source_code/qr_sweep.m}
\subsubsection{\quotes{Стандартный} алгоритм}
\lstinputlisting{source_code/svd_standard.m}
\newpage % =======================================================================================================================================
\section{Приложения сингулярного разложения}
\subsection{Вычисление псевдообратной матрицы}
\subsubsection{Теория}
При вычислении обратной матрицы возникает 2 проблемы:
\begin{enumerate}
\item Вычисление обратной матрицы возможно только для квадратных матриц
\item Если определитель матрицы равен нулю, то при вычислении обратной матрицы происходит деление на ноль.
\end{enumerate}
Эти проблемы можно обойти, используя сингулярное разложение и вычисляя псевдообратную матрицу.
Представим, что дана матрица $A_{m \times n},\ m \ge n$, причем $r = \rank{(A)},\ r = n$, тогда можно утверждать, что $\exists\ (A^TA)^{-1}$. Пусть имеется система уравнений: \[ Ax = b, \] тогда, домножив слева на $A^T$, мы получим уравнение вида: \[ A^TAx=A^Tb. \] Т.к. $\exists\ (A^TA)^{-1}$, то можно выразить $x$ и тем самым определить псевдообратную матрицу: \begin{equation} \label{app:pseudo:def:A}
x = (A^TA)^{-1}A^Tb \Rightarrow A^+ = (A^TA)^{-1}A^T.
\end{equation}
В данном примере требовалось, чтобы $r = \rank{(A)},\ r = n$. Предположим теперь, что $r < n$. Вспомним, что \quotes{полное} сингулярное разложение имеет следующий вид: \begin{equation} \label{app:pseudo:SVD_full}
A_{m \times n} = U_{m \times m} \Sigma_{m \times n} V^T_{n \times n}.
\end{equation}
Подставим (\ref{app:pseudo:SVD_full}) в (\ref{app:pseudo:def:A}) и произведем необходимые сокращения: \[ A^+= (V \Sigma U^TU \Sigma^T V^T)^{-1}V \Sigma U^T = (V \Sigma^T\Sigma V^T)^{-1}V \Sigma^T U^T = V (\Sigma^T\Sigma)^{-1}\Sigma^T U^T = V \Sigma^{+} U^T. \]
\[ A^+ = V_{n \times n} (\Sigma_{m \times n})^{+} U^T_{m \times m} = V_{n \times n} (\Sigma^{+})_{n \times m} U^T_{m \times m}. \]
Основная проблема -- вычислить $\Sigma^{+}$. Эта операция выполняется в два шага:
\begin{enumerate}
\item Матрица транспонируется
\item Все ненулевые элементы главной диагонали возводятся в степень $-1$.
\end{enumerate}
Т.о. вычисление псевдообратной матрицы сводится к вычислению сингулярного разложения.
\subsubsection{Исходный код}
\lstinputlisting{source_code/p_inv.m}
\subsection{Метод главных компонент}% ==============================================================================
\label{applacation:pca}
Метод Главных Компонент -- один из основных способов уменьшить размерность данных, потеряв наименьшее количество информации. Изобретен К.Пирсоном в 1901 г. Применяется во многих областях, таких как распознавание образов, компьютерное зрение, сжатие данных и т.п.
\subsubsection{Теория}
Существует две основные матричные нормы в численных методах линейной алгебры: норма Фробениуса -- $||A||_F$ и 2-норма -- $||A||_2$. Норма Фробениуса была подробно рассмотрена в (\ref{theory:sing_vectors}), поэтому перейдем сразу к рассмотрению 2-нормы.
2-норма определяется следующим образом: \[ ||A||_2 = \max_{||\mathbf{v}||=1}{||A\mathbf{v}||} \Rightarrow ||A||_2 = \sigma_1. \]
Представим матрицу $А$ размеров $n \times d$ как $n$ точек $\{a^{(i)}\}_{i=1}^n$ в d-мерном пространстве. Норма Фробениуса матрицы $A$ эквивалентна квадратному корню суммы квадратов расстояний от точек $\{a^{(i)}\}_{i=1}^n$ до начала координат. 2-норма матрицы $A$ эквивалентна квадратному корню суммы квадратов расстояний точек $\{a^{(i)}\}_{i=1}^n$ до прямой, проходящей вдоль вектора $\mathbf{v_1}$.
Пусть \[A = \sum_{i=1}^r \sigma_i \mathbf{u_i} \mathbf{v_i^T}\] сингулярное разложение матрицы $A$. Пусть для $k \in \{1,2,...,r\}$ \[ A_k = \sum_{i=1}^k \sigma_i \mathbf{u_i} \mathbf{v_i^T} \] является $k$-ой частичной суммой. В таком случае матрица $A_k$ имеет ранг $k$, а также она является наилучшей аппроксимацией ранга $k$ матрицы $A$, если ошибка измеряется с помощью 2-нормы или нормы Фробениуса.
Докажем лемму, утверждающую, что строки матрицы $A_k$ являются проекциями строк матрицы $A$ на подпространство $V_k$, образованное первыми $k$ правыми сингулярными векторами $\{\mathbf{v_i}\}_{i=1}^k$ матрицы $A$.
Пусть $\mathbf{a}$ -- произвольная вектор-строка. Т.к. векторы $\mathbf{v_i}$ ортонормированы, то проекция вектора $\mathbf{a}$ на подпространство $V_k$ определяется следующим соотношением: $\sum_{i=1}^{k}(\mathbf{a} \cdot \mathbf{v_i})\mathbf{v_i^T}$. Таким образом, матрица, строки которой являются проекциями строк матрицы $A$ на подпространство $V_k$, задается выражением $\sum_{i=1}^{k}A\mathbf{v_i}\mathbf{v_i^T}$. \[ \sum_{i=1}^{k}A\mathbf{v_i}\mathbf{v_i^T} = \sum_{i=1}^{k}\sigma_i \mathbf{u_i} \mathbf{v_i^T} = A_k. \]
Т.к. строки матрицы $A_k$ -- проекции строк матрицы $A$ на $k$-размерное подпространство $V_k$, то $\rank{(A_k)} \le k$.
Докажем теорему, утверждающую, что для любой матрицы $B$ ранга не более, чем $k$ \[ ||A-A_k||_F \le ||A-B||_F. \]
Пусть $B$ минимизирует $||A-B||_F^2$ среди всех матриц ранга не больше, чем $k$. Пусть $V$ -- подпространство, порожденное векторами-строками матрицы $B$. Размерность подпространства $V$ не превосходит $k$. Т.к. $B$ минимизирует $||A-B||_F^2$, то отсюда следует, что строки матрицы $B$ -- проекции строк матрицы $A$ на подпространство $V$. Т.к. строки матрицы $B$ являются проекциями соответствующих строк матрицы $A$, то $||A-B||_F^2$ является суммой квадратов расстояний строк матрицы $A$ на подпространство $V$. Т.к. $A_k$ минимизирует сумму квадратов расстояний строк матрицы $A$ на подпространство размерности не больше $k$, то \[ ||A-A_k||_F^2 \le ||A-B||_F^2 \Rightarrow ||A-A_k||_F \le ||A-B||_F. \]
Докажем лемму, утверждающую, что \[ ||A-A_k||_2^2 = \sigma_{k+1}^2. \]
Пусть $A = \sum_{i=1}^{r}\sigma_i \mathbf{u_i} \mathbf{v_i^T}$ -- сингулярное разложение матрицы $A$. В таком случае $A_k = \sum_{i=1}^{k}\sigma_i \mathbf{u_i} \mathbf{v_i^T}$ и \[A-A_k = \sum_{i=k+1}^{r}\sigma_i \mathbf{u_i} \mathbf{v_i^T}\]
Вспомним, что \[ ||A|| = \sup_{||x||=1} ||Ax||. \]
В таком случае \[ ||A-A_k||_2 = \sup_{||\mathbf{x}||_2=1}||(A-A_k)\mathbf{x}||_2 = \bigg| \bigg| \sum_{i=k+1}^{r}\sigma_i \mathbf{u_i} \mathbf{v_i^T} \mathbf{x} \bigg| \bigg|_2. \]
Представим вектор $\mathbf{x}$ в виде линейной комбинации векторов $\mathbf{v_1},...,\mathbf{v_r}$: \[ \mathbf{x} = \sum_{i=1}^r \alpha_i \mathbf{v_i}. \]
В таком случае \[ ||(A-A_k)\mathbf{x}||_2 = \bigg| \bigg| \sum_{i=k+1}^{r}\sigma_i \mathbf{u_i} \mathbf{v_i^T} \sum_{i=1}^r \alpha_i \mathbf{v_i} \bigg| \bigg|_2 = \bigg| \bigg| \sum_{i=k+1}^{r} \alpha_i \sigma_i \mathbf{u_i} \mathbf{v_i^T} \mathbf{v_i} \bigg| \bigg|_2 = \bigg| \bigg| \sum_{i=k+1}^{r} \alpha_i \sigma_i \mathbf{u_i} \bigg| \bigg|_2. \]
Последние преобразования возможны благодаря тому, что
\begin{equation} \nonumber
\mathbf{v_i^Tv_j} = \left\{ \begin{aligned}
& 1,\ i = j \\
& 0,\ i \neq j
\end{aligned} \right.
\end{equation}
Таким образом получаем, что \[ ||(A-A_k)\mathbf{x}||_2 = \bigg| \bigg| \sum_{i=k+1}^{r} \alpha_i \sigma_i \mathbf{u_i} \bigg| \bigg|_2 = \sqrt{\sum_{i=k+1}^{r} \alpha_i^2 \sigma_i^2}, \]
т.к. получившееся выражение есть ни что иное, как разложение вектора по базису $\mathbf{u_i}$.
Учитывая тот факт, что последовательность $\{ \sigma_i \}_{i=1}^r$ не возрастает,
\begin{equation} \nonumber
\left. \begin{aligned}
& \mathbf{x} = (\underbrace{0,...,0}_{k},1,0,...,0)^T,\ ||\mathbf{x}||_2 = 1 \\
& \mathbf{x} = \argmax_{||\mathbf{x}||_2 = 1}{||(A-A_k)\mathbf{x}||_2}
\end{aligned} \right\} \Rightarrow ||A-A_k||_2^2 = \sigma_{k+1}.
\end{equation}
Докажем теорему, утверждающую, что для любой матрицы $B$ такой, что $\rank(B) \le k$ выполняется неравенство: \[||A-A_k||_2 \le ||A-B||_2.\]
Если $\rank(A) \le k$, то теорема, очевидно, выполняется, т.к. $||A-A_k||_2 = 0$. Таким образом, предположим, что $\rank(A) > k$. Используя вышедоказанную лемму, получим, что $||A-A_k||_2 = \sigma_{k+1}^2$. Теперь предположим, что $\exists B,\ \rank(B) \le k$, для которой \[ ||A-B||_2 \le ||A-A_k||_2\ \Rightarrow\ ||A-B||_2 \le \sigma_{k+1}^2. \]
Будем рассматривать матрицу $B$ как линейное преобразование. В таком случае размерность ядра преобразования $B$: \[\ker(B) = \{ \mathbf{x}\ |\ B\mathbf{x} = 0 \},\] будет не меньше, чем $(d-k)$\[\dim \ker(B) \ge (d-k).\]
Пусть векторы $\mathbf{v_1},...,\mathbf{v_{k+1}}$ -- первые $k+1$ сингулярные векторы матрицы $A$. Обозначим $d-k$ независимых векторов ядра преобразования $B$ как $\mathbf{u_1},...,\mathbf{u_{d-k}}$. В таком случае $\mathbf{u_1},...,\mathbf{u_{d-k}}, \mathbf{v_1},...,\mathbf{v_{k+1}}$ -- набор из $d+1$ векторов в $d$-мерном пространстве, следовательно пространства $V_k = \Span\{\mathbf{v_1},...,\mathbf{v_{k+1}}\}$ и $\ker(B) = \Span\{\mathbf{u_1},...,\mathbf{u_{d-k}}\}$ обязательно пересекаются: \[ \dim \bigg( \ker(B) \cap \Span\{\mathbf{v_1},...,\mathbf{v_{k+1}}\} \bigg) \neq 0\ \Rightarrow\ \exists \mathbf{z} \neq 0:\ \mathbf{z} \in \bigg( \ker(B) \cap \Span\{\mathbf{v_1},...,\mathbf{v_{k+1}}\} \bigg). \]
Выберем $\mathbf{z}$ таким образом, чтобы $||\mathbf{z}||_2 = 1$. Вспомнив определение нормы матрицы, получим: \[ ||A|| = \sup_{||\mathbf{x}||=1}||A\mathbf{x}|| \Rightarrow ||A-B||_2^2 \ge ||(A-B)\mathbf{z}||_2^2. \]
Учитывая тот факт, что $\mathbf{z} \in \ker(B):\ B\mathbf{z} = 0$, получаем, что \[||A-B||_2^2 \ge ||(A-B)\mathbf{z}||_2^2 = ||A\mathbf{z}||_2^2.\]
Учитывая тот факт, что $\mathbf{z} \in V_k = \Span\{\mathbf{v_1},...,\mathbf{v_{k+1}}\}$, получаем, что
\begin{equation} \nonumber
\begin{aligned}
||A\mathbf{z}||_2^2 &= \Bigg| \Bigg| \sum_{i=1}^{n}\sigma_i \mathbf{u_i} \mathbf{v_i^T}\mathbf{z} \Bigg| \Bigg|_2^2 = \Bigg| \Bigg| \sum_{i=1}^{n}\sigma_i \mathbf{u_i} \mathbf{v_i^T} \bigg( \sum_{j=1}^{k+1} \alpha_j \mathbf{v_j} \bigg) \Bigg| \Bigg|_2^2 = \Bigg| \Bigg| \sum_{i=1}^{k+1}\sigma_i \alpha_i \mathbf{u_i} \Bigg| \Bigg|_2^2 \\
&= \sum_{i=1}^{k+1}\sigma_i^2 \alpha_i^2 \ge \sigma_{k+1}^2 \sum_{i=1}^{k+1} \alpha_i^2 = \sigma_{k+1}^2.
\end{aligned}
\end{equation}
Таким образом получили выражение: \[ ||A-B||_2^2 \ge \sigma_{k+1}^2, \] что противоречит изначальному предположению о том, что $||A-B||_2 < \sigma_{k+1}$, следовательно \[\forall\ B,\ \rank(B) \le k:\ ||A-A_k||_2 \le ||A-B||_2.\]
\subsubsection{Исходный код}
\lstinputlisting{source_code/low_rank.m}
\subsubsection{Пример работы}
\begin{figure}[!h]
\centering
\includegraphics[width=0.65\textwidth]{low_rank.png}
\caption{Пример работы} \label{pic:pca:example}
\end{figure}
На (Рис.~\ref{pic:pca:example}) можно видеть, как точки, случайным образом разбросанные в трехмерном пространстве (синие) спроецировались на двумерное подпространство (красное).
\newpage
\subsection{Понижение размерности}
\subsubsection{Постановка задачи}
Ярким примером задачи понижения размерности является задача из области машинного обучения.
Будем рассматривать определенный класс объектов, каждый из которых можно представить с помощью конечного набора количественных и качественных характеристик, которые образуют вектор $(x_1,...,x_N)^T$. Рассмотрим набор, состоящий только из $n$ количественных характеристик $(n \le N)$. Т.о. каждому объекту ставится в соответствие вектор действительных чисел $(x_{i_1},...,x_{i_n})^T$, $\{i_1,...,i_n\} \subset \{1,...,N\}$. Без ограничения общности будем предполагать, что $(i_1,...,i_n) \equiv (1,...,n)$.
Предположим, что обучающая выборка содержит $M$ различных объектов, каждый из которых описывается вектором $n$ действительных чисел. В таком случае каждый объект можно представить радиус-вектором $n$-мерного пространства, а всю обучающую выборку -- матрицей $A_{n\times M}$.
Часто бывает так, что $n$ достаточно велико и работать с матрицей обучающей выборки либо неудобно, либо тяжело. В таком случае есть смысл попробовать понизить размерность матрицы, потеряв при этом наименьшее количество информации: рассматривать вместо вектора $(x_1,...,x_n)^T$, вектор $(y_1,...,y_k)^T$, причем $k < n$.
Однако проблема стоит не только в том, чтобы понизить размерность матрицы обучающей выборки, но также разработать алгоритм перевода любого вектора $(x_1,...,x_n)^T$ в соответствующий ему вектор $(y_1,...,y_k)^T$.
\subsubsection{Теория}
В данном разделе будет приведен алгоритм понижения размерности набора свойств класса объектов: \[ (x_1,...,x_n)^T \rightarrow (y_1,...,y_k)^T. \]
\begin{figure}[!h]
\centering
\includegraphics[width=0.95\textwidth]{img/dim_reduction.jpg}
\caption{Частный случай} \label{pic:dim_reduction:example}
\end{figure}
-- Первым делом необходимо понизить ранг матрицы $A_{n\times M}$ до $k < \min(n,M)$, воспользовавшись алгоритмом, описанным в (\ref{applacation:pca}).
После того, как матрица $A_{n\times M}$ понижена до ранга $k$, каждый столбец можно представить точкой на $k$-мерной гиперплоскости в $n$-мерном пространстве. Визуально представить эффект понижения порядка можно, посмотрев на (Рис.~\ref{pic:pca:example}).
Разложим матрицу $A_{n\times M}$ с помощью QR-разложения. Полученные матрицы $Q$ и $R$ являются матрицами размерностей $n \times n$ и $n \times M$ соответственно. Однако, учитывая то, что ранг матрицы $A_{n\times M}$ равен $k < n$, верхнетреугольная матрица $R$ имеет $n-k$ нулевых строк, т.е. $n-k$ характеристик обнулились. Необходимо запомнить преобразование $Q$, обнуляющее $n-k$ характеристик. Оно выглядит следующим образом: \[ Y := R = Q^TA. \]
-- На втором шаге необходимо выбрать \quotes{опорную} точку $P_0$, представляющую из себя рандомный вектор-столбец матрицы $A$, пониженного ранга.
-- Третим пунктом алгоритма является задание базиса $k$-мерной гиперплоскости $\Pi$. Для этого необходимо случайным образом выбрать $k$ точек $\{P_1,...,P_k\}$ таким образом, чтобы \[ \forall~ i \neq j:\ P_i \neq P_j,\ \ B_{n \times k} = \bigg[ P_1 - P_0\ \Big|\ ...\ \Big|\ P_k - P_0 \bigg],\ \ \rank{(B)} = k. \]
Фактически данный базис задает смещенную плоскость, проходящую через начало координат в точке $P_0$. Т.о. произошла неявная смена системы координат: $Oxyz \rightarrow P_0xyz$.
-- Далее рассмотрим произвольную точку $y$, которую необходимо спроецировать на гиперплоскость $\Pi$.
Изначально представим ее в новой, смещенной системе координат: \[\tilde{y} = y - P_0.\]
Обозначим базисные векторы \[P_i-P_0 = z_i = (z_i^1,...,z_i^n)^T.\]
Т.о. матрица $B_{n\times k}$ примет следующий вид: \[ B = \begin{bmatrix}
z_1 & \hdots & z_k
\end{bmatrix} =
\begin{bmatrix}
z_1^1 & \hdots & z_k^1 \\
\hdots & \hdots & \hdots \\
z_1^n & \hdots & z_k^n
\end{bmatrix}
\]
Если бы вектор $\tilde{y}$ лежал в плоскости $\Pi$, то $B\tilde{x} = \tilde{y}$. Однако в конкретном случае $\tilde{y}$ не обязан лежать в плоскости $\Pi$, поэтому сперва необходимо спроецировать его на эту плоскость с помощью псевдообратной матрицы $B^+$: \[ \tilde{x} = B^+\tilde{y} = B^+(y - P_0). \]
Полученный вектор $\tilde{x}$ задан в базисе плоскости $\{z_1,...,z_k\}$, поэтому необходимо перейти сперва к системе координат $P_0xyz$, а затем и в исходную $Oxyz$.
Поэтому искомый вектор $x$ примет окончательный вид: \[ x = B\tilde{x} + P_0 = BB^+\tilde{y} + P_0 = BB^+(y - P_0) + P_0 = BB^+y - BB^+P_0 + P_0 = BB^+y. \]
\[x = BB^+y.\]
Т.о. полностью проанализирована задача перехода от $y$ к $\tilde{x}$ во всех системах координат.
\newpage
\subsection{Сжатие изображения}
Сингулярное разложение оказывается очень удобным механизмом сжатия изображений с минимальной потерей информации. В данном разделе будет подробно разобран и запрограммирован простейший алгоритм сжатия серых изображений с использованием сингулярного разложения.
\subsubsection{Принцип работы} \label{compression}
Мы можем представить изображение как матрицу следующим образом. Предположим, что изображение имеет $n\times m$ пикселей. Оттенки серого кодируются одним числом, поэтому все изображение представляется обычной матрицей размерности $n\times m$.
У цветных изображений на каждый пиксель приходится три числовых значения, поэтому все изображение можно представить тремя матрицами размерностей $n \times m$, либо одной матрицей размерности $n \times 3m$, получающейся в результате конкатенации трех матриц, каждая из которых отвечает за свой цвет: $A = [A_{\text{red}}\ |\ A_{\text{green}}\ |\ A_{\text{blue}}]$.
Рассмотрим сперва случай квадратного изображения в градации серого, которое представляется матрицей $I_0$ размерности $n \times n$. В подавляющем большинстве случаев $\rank(I_0) = n$.
Матрицу $I_0$ можно разложить с помощью сингулярного разложения: \[I_0 = \sum_{i=1}^n \sigma_i \mathbf{u_i} \mathbf{v_i^T}. \]
Аппроксимируем матрицу $I_0$ матрицей \[I^{(k)} = \sum_{i=1}^k \sigma_i \mathbf{u_i} \mathbf{v_i^T}.\]
Матрица $I^{(k)}$ имеет ранг $k$, т.е. линейно независимыми являются лишь $k$ строк и столбцов.
Проанализируем полученный результат: мы аппроксимировали матрицу $I_0$ размерности $n \times n$ и ранга $n$ матрицей $I^{(k)}$ размерности $n \times n$ и рангом $k$, причем погрешность аппроксимации составляет \[ \delta_I = \dfrac{||I_0-I^{(k)}||_2}{||I_0||_2} \cdot 100\%\ \ \ \ \ \ \ \ \Delta_I = \dfrac{||I^{(k)}||_F}{||I_0||_F}. \]
Чем больше значение $k$, тем лучше $I^{(k)}$ аппроксимирует $I_0$. Из приведенных двух оценок аппроксимации принято использовать $\Delta_I$, т.к. она дает более точную и логичную оценку. На (Рис.~\ref{pic:compression:errs}) показаны изменения погрешностей в зависимости от выбранного ранга $k$ для изображения (Рис.~\ref{pic:compression:original}).
На данный момент изображение потеряло часть информации (возможно, пропала четкость и появилось размытие), однако само изображение как кодировалось матрицей $n \times n$, так и продолжает кодироваться матрицей такого же размера. Необходимо каким-то образом уменьшить размерность матрицы. Т.к. матрица $I^{(k)}$ имеет ранг $k$, то это значит, что она содержит $k$ линейно независимых столбцов. Следовательно, ортогонализовав матрицу $I^{(k)}$, мы получим матрицу с $n-k$ нулевыми строками или столбцами (в зависимости от того, как ортогонализуем).
Разложим матрицу $I^{(k)}$ с помощью QR-разложения. Полученные матрицы $Q$ и $R$ являются матрицами размерностей $n \times n$. Однако, учитывая то, что ранг матрицы $I^{(k)}$ равен $k < n$, верхнетреугольная матрица $R$ имеет $n-k$ нулевых строк, которые можно исключить из матрицы $R$. Исключив из $R$ $n-k$ последних строк, необходимо также исключить из матрицы $Q$ $n-k$ последних столбцов для согласованности матриц $Q,R$. Таким образом,
\begin{equation}
\left.
\begin{aligned}
I^{(k)}_{n\times n} = Q_{n\times n}R_{n\times n}\\
Q_{n\times n} \rightarrow \widetilde{Q}_{n\times k}\\
R_{n\times n} \rightarrow \widetilde{R}_{k\times n}
\end{aligned} \right\}\ \Rightarrow\ I^{(k)} = \widetilde{Q}\widetilde{R}.
\end{equation}
Таким образом, мы закодировали матрицу $n\times n$ двумя матрицами $n\times k$.
Абсолютно аналогичными рассуждениями прямоугольная матрица $n\times m$ сводится к двум матрицам с размерностями $n\times k$ и $k\times m$.
\subsubsection{Сжатие}
Рассмотрим случай квадратной матрицы $I_{n\times n}$. Применив алгоритм, описанный в (\ref{compression}), мы свели матрицу $I_{n\times n}$ к двум матрицам $\widetilde{Q}_{n\times k},\ \widetilde{R}_{k\times n}$. Т.о. $n^2$ элементов теперь кодируется $2nk$ элементами. Проанализируем, в каком случае происходит сжатие объема требуемой памяти.
\[ 2nk < n^2\ \Rightarrow\ 2k < n\ \Rightarrow\ k < \frac{n}{2}. \]
Т.е. сжатие изображения происходит, если выбираемый параметр $k < \dfrac{n}{2}$.
В таком случае целесообразно ввести коэффициент сжатия следующим образом: \[\alpha = \dfrac{n}{2k}.\]
Рассмотрим теперь случай прямоугольной матрицы $I_{n\times m}$. Применив алгоритм, описанный в (\ref{compression}), мы свели матрицу $I_{n\times m}$ к двум матрицам $\widetilde{Q}_{n\times k},\ \widetilde{R}_{k\times m}$. Т.о. $nm$ элементов теперь кодируется $(n+m)k$ элементами. Проанализируем, в каком случае происходит сжатие объема требуемой памяти.
\[ (n+m)k < nm\ \Rightarrow\ k < \frac{nm}{n+m}. \]
Т.е. сжатие изображения происходит, если выбираемый параметр $k < \dfrac{nm}{n+m}$.
Введем коэффициент сжатия аналогично тому, как он был введен в случае с квадратной матрицей: \[ \alpha = \frac{nm}{(n+m)k}. \]
\subsubsection{Исходный код}
\lstinputlisting{source_code/image_compression.m}
\subsubsection{Пример работы}
\begin{itemize}
\item Применим вышеприведенный код к изображению (Рис.~\ref{pic:compression:original}).
\item После аппроксимации матрицей меньшего ранга ($\alpha = 3$) получим изображение (Рис.~\ref{pic:compression:approximation}).
\item Сжав изображение, получим новое изображение (Рис.~\ref{pic:compression:compressed}).
\item На (Рис.~\ref{pic:compression:diff}) можно видеть разность оригинального изображения и его аппроксимации, к которой применен алгоритм \quotes{Histogram equalization}.
\end{itemize}
\begin{figure}[!h]
\centering
\includegraphics[width=0.65\textwidth]{img/img_compression/1.png}
\caption{Исходное изображение} \label{pic:compression:original}
\end{figure}
\begin{figure}[!h]
\centering
\includegraphics[width=0.75\textwidth]{img/img_compression/2.png}
\caption{Аппроксимация изображения} \label{pic:compression:approximation}
\end{figure}
\begin{figure}[!h]
\centering
\includegraphics[width=0.75\textwidth]{img/img_compression/3.png}
\caption{Сжатое изображение (Рис.~\ref{pic:compression:approximation})} \label{pic:compression:compressed}
\end{figure}
\begin{figure}[!h]
\centering
\includegraphics[width=0.7\textwidth]{img/img_compression/4.png}
\caption{Разность изображений} \label{pic:compression:diff}
\end{figure}
\begin{figure}[!h]
\centering
\includegraphics[width=0.7\textwidth]{img/img_compression/5.png}
\caption{Поведение ошибок изображений} \label{pic:compression:errs}
\end{figure}
\newpage % =======================================================================================================================================
\section{Дополнение}
\subsection{Определения}
Эрмитово-сопряжённая мтрица или сопряжённо-транспонрованная мтрица — это матрица $A^H$ с комплексными элементами, полученная из исходной матрицы $A$ транспонированием и заменой каждого элемента комплексно-сопряжённым ему. Эрмитово-сопряжённые матрицы во многом играют ту же роль при изучении комплексных векторных пространств, что и транспонированные матрицы в случае вещественных пространств.
{\vspace{0.5cm}}
Унитрная мтрица -- квадратная матрица с комплексными элементами, результат умножения которой на эрмитово-сопряжённую равен единичной матрице: $U^\dagger U = UU^\dagger = I$. Другими словами, матрица унитарна тогда и только тогда, когда существует обратная к ней матрица, удовлетворяющая условию $U^{-1}=U^\dagger$. Унитарная матрица, элементы которой вещественны, является ортогональной.
{\vspace{0.5cm}}
Норма Фробениуса -- субмультипликативная норма ($||AB||_F \leq ||A||_F||B||_F$), часто используемая в численных методах линейной алгебры из-за своих свойств: во-первых, она часто проще в вычислении, чем индуцированные нормы, а во-вторых, она инвариантна относительно поворотов $\forall R \in SO(n): ||A||_F = ||AR||_F = ||RA||_F$.
\subsection{Извлечение квадратного корня из матрицы}
Извлечение квадратного корня из матрицы возможно для любой матрицы, имеющей полный набор собственных векторов и собственных значений, и у которой все собственные значения положительны.
Корень квадратный из матрицы -- матрица, имеющая такой же набор собственных векторов, что и исходная матрица, а собственные значения которой являются квадратными корнями из собственных значений исходной матрицы.
Таким образом, чтобы посчитать $\sqrt{M}$, необходимо с помощью спектрального разложения представить матрицу M в виде $M = V \Lambda V^{-1}$ и тогда:
\[
\sqrt{M} = V \sqrt{\Lambda} V^{-1}.
\]
\newpage
\begin{thebibliography}{8}
\addcontentsline{toc}{section}{\refname}
\bibitem{citation1} J.~Hopcroft, R.~Kannan -- Computer Science Theory for the Information Age, 2012.
\bibitem{citation2} G.~W.~Stewart -- On the Early History of the Singular Value Decomposition.
\bibitem{citation3} J.~Demmel, K.~Veselic -- Jacobi's method is more accurate than QR.
\bibitem{citation4} J.~Krishnamoorthy, K.~Kocagoez -- Singular Values using Cholesky Decomposition.
\bibitem{citation5} N.~J.~Higham -- Analysis of the Cholesky Decomposition of a Semi-Definite Matrix.
\bibitem{citation6} J.~Demmel и W.~Kahan -- Accurate Singular Values of Bidiagonal Matrices.
\bibitem{citation7} M.~Plesinger, Z.~Strakos -- Some remarks on bidiagonalization and its implementation.
\bibitem{citation8} S.~Kahu, R.~Rahate -- Image Compression using Singular Value Decomposition, International Journal of Advancements in Research \& Technology, Volume 2, Issue 8, August-2013 244 ISSN 2278-7763.
\bibitem{citation9} G.~Golub и C.~Reinsch -- Singular Value Decomposition and Least Squares Solutions.
\end{thebibliography}
\end{document}