Очередной семинар - ==Аппроксимация матрицы суммой одноранговых матриц== ~~TOC~~ ===Постановка задачи== Скалярное произведение в вещественном пространстве $ (m \times n) $-матриц будем считать заданным равенством $$ \langle A, B \rangle := \sum_{j,k} a_{jk} b_{jk} \ \mbox{при} \ \{A,B\} \subset \mathbb R^{m\times n}. $$ как если бы эти матрицы были предварительно ((:algebra2#vektorizacija векторизованы)), и затем перемножены скалярно как векторы. Эквивалентное определение --- через ((:algebra2#sled след)) матрицы: $$ \langle A, B \rangle = \operatorname{Sp} (A \cdot B^{\top}) = \operatorname{Sp} (A^{\top} \cdot B) \, . $$ Как следствие такого способа определения евклидова пространства, вводится понятие евклидовой нормы (или нормы Фробениуса) $$ \| B \|:=\sqrt{\sum_{j,k} b_{jk}^2} = \sqrt{\operatorname{Sp}(B^{\top} B)} = \sqrt{\operatorname{Sp}(B\cdot B^{\top})} \, . $$ Естественно, что расстояние между матрицами $ A $ и $ B $ это обобщение расстояния между точками в стандартном его определении: $$ \|A-B \|:= \sqrt{\sum_{j,k} (a_{jk}-b_{jk})^2} = \sqrt{\operatorname{Sp}((A-B)^{\top}(A- B))} \, . $$ **Задача.** Для заданной матрицы $ A \in \mathbb R^{m\times n} $ и натурального числа $ \mathfrak r < \operatorname{rank} (A) $ построить матрицу $ \widetilde A \in \mathbb R^{m\times n} $ ранга $ \mathfrak r $ наилучшим образом приближающую матрицу $ A $ в ((:norm_space#норма_матрицы евклидовой норме (норме Фробениуса))): $$ \min_{\widetilde A\in \mathbb R^{m\times n}, \operatorname{rank} (\widetilde A)= \mathfrak r} || A- \widetilde A|| \, . $$ Иными словами, среди всех матриц ранга $ \mathfrak r $ найти матрицу, ближайшую к матрице $ A $ (иными словами: найти ортогональную проекцию матрицы $ A $ на множество всех матриц ранга $ \mathfrak r $). В чем важность этой задачи? --- Дело в том, что если произвольная матрица $ A \in \mathbb R^{m\times n} $ задается $ m\cdot n $ элементами, то ((:algebra2/rank#matricy_ranga_1 одноранговая матрица)) того же порядка задается $ (m+n-1) $ элементами, двухранговая --- $ (m+n-1)+ (m+n-3) $ элементами, и т.д. Чем меньше $ \mathfrak r $ --- тем меньше требуется ресурсов для хранения. В чем сложности? --- 1. При фиксированном $ \mathfrak r $ найти матрицу $ \widetilde A $. 2. Определить оптимальную величину $ \mathfrak r $, т.е., с одной стороны, чтобы эта величина была поменьше, но, с другой стороны, чтобы соответствующая матрица $ \widetilde A $ давала хорошее приближение матрицы $ A $. Например, чтобы по построенным приближениям $ \widetilde A $ и $ \widetilde B $ матриц $ A $ и $ B $ можно было бы с большой вероятностью правильно идентифицировать эти матрицы. ===Рабочий аппарат== Представление матрицы $ A $ в виде суммы линейно независимых одноранговых матриц. Такое представление всегда существует в количестве слагаемых равном $ \operatorname{rank} (A) $. К примеру, см. ((:algebra2/rank#matricy_ranga_1 ЗДЕСЬ)). И подобных представлений можно привести бесконечно много. Наша цель: найти такое представление, где слагаемые будут ранжированы по "объему вклада" в матрицу $ A $: первая матрица ранга $1$ должна быть максимально близка к матрице $ A $, вторая матрица --- максимально близка к разности $ A $ и первой матрицы и т.д. Требуемое представление получается посредством **сингулярного разложения матрицы** $ A $. Формальный подход к изложению теории см. ☞ ((:algebra2/svd ЗДЕСЬ)). Это разложение используется не только для решения поставленной в настоящем разделе задачи. Если бегло посмотреть на ((:algebra2/svd#singuljarnoe_razlozhenie основной результат)) теории, то увидеть представление матрицы $ A $ в виде //суммы// каких-то матриц затруднительно: там стоит чистое //произведение//, т.е. разложение матрицы $ A $ на три сомножителя, из которых один (средний) --- простой структуры (диагональная матрица). На самом деле, можно показать (см. ☟ ((#singuljarnoe_razlozhenie_v_vide_proizvedenija_matric НИЖЕ)) )эквивалентность этого произведения именно сумме одноранговых матриц. Обратим внимание на похожесть упоминаемого результата на уже пройденный в курсе алгебры результат о представлении квадратной матрицы $ A $ также в виде произведения трех матриц $$ A=C^{-1} A_{diag} C \, , $$ где матрица $ A_{diag} $ --- диагональная, с собственными числами матрицы $ A $ на главной диагонали. Не акцентируем пока внимания на то, что такое представление не для всех матриц возможно; считаем, что случайным образом выбранная матрица $ A $ ((:mapping/operator#diagonalizuemost_matricy_operatora диагонализуема)). Нас, в первую очередь, интересует другое: является ли этот результат частью общего результата о сингулярном разложении? Для некоторого класса матриц ответ оказывается положительным. Вот именно на основе этих матриц мы и сделаем сначала "заготовку" для общей теории. ===Квадратные матрицы== Как произведении $ A=C^{-1} A_{diag} C $ "перегнать" в сумму одноранговых матриц? !!Т!! **Теорема 1.** //Если все собственные числа// $ \lambda_{1},\dots,\lambda_n $ //матрицы// $ A \in \mathbb C^{n\times n} $ //различны, то для произвольного полинома// $ g_{}(x) $ //справедлива формула// $$ g(A)= \sum_{k=1}^n \frac{g(\lambda_k)}{f_k(\lambda_k)}f_k(A) = \sum_{k=1}^n \frac{g(\lambda_k)}{f'(\lambda_k)}f_k(A) $$ //где// $ f_k(x) $ //означает частное от деления характеристического полинома// $ f_{}(x):=\det (A-xE_{n\times n}) $ //на// $ x- \lambda_k $. Поскольку этот результат справедлив для произвольного полинома $ g(x) $, то он выполняется и для частного случая $ g(x) \equiv x $: $$ A = \sum_{k=1}^n \lambda_k \frac{f_k(A)}{f_k(\lambda_k)} = \sum_{k=1}^n \lambda_k \frac{f_k(A)}{f^{\prime}(\lambda_k)} \, . $$ То есть мы получили представление матрицы $ A $ в виде линейной комбинации матриц $$ \left\{ M_k(A):=\frac{f_k(A)}{f_k(\lambda_k)} \right\}_{k=1}^n \, . $$ Теперь исследуем свойства этих матриц --- они называются **ковариантами Фробениуса матрицы** $ A $. !!Т!! **Теорема 2.** //Если алгебраическая кратность собственного числа// $ \lambda_k $ //равна// $ 1 $, //то// $$ \operatorname{rank} (M_k(A)) = 1 \, . $$ Любая одноранговая $m\times n $-матрица ((:algebra2/rank#matricy_ranga_1 представима в виде произведения)) столбца $ U \in \mathbb C^{m\times 1} $ на строчку $ V\in \mathbb C^{1\times n} $: $$ \left(\begin{array}{c} u_1\\ u_2\\ \vdots \\ u_m \end{array} \right) (v_1,\dots,v_n)=\left[u_jv_k\right]_{j=1,\dots,m;k=1,\dots,n} $$ при хотя бы одной паре индексов $ j_{} $ и $ k_{} $ таких, что $ u_j\ne 0, v_k \ne 0 $. Для матрицы $ f_k(A) $ такое представление обеспечивается собственными векторами матрицы $ A $. !!Т!!** Теорема 3.** //Матрица// $ M_k(A) $ //представима в виде произведения// $$ M_k(A) = X_k \cdot Y_k \, , $$ где столбец $ X_k $ ---// ((:algebra2/charpoly#sobstvennyj_vektor правый собственный вектор)), а строка $ Y_k $ --- ((:algebra2/charpoly#sobstvennyj_vektor левый собственный вектор)) матрицы// $A $, //принадлежащие// $ \lambda_k $. **Доказательство** следует из теоремы Гамильтона-Кэли: $$ f(A)=\mathbb O_{n\times n} \quad \Rightarrow \ (A-\lambda_k E) f_k(A)= \mathbb O_{n\times n} \, . $$ То есть любой ненулевой столбец матрицы $ f_k(A) $ должен быть правым собственным вектором матрицы $ A $, принадлежащим $ \lambda_k $. Аналогично --- для любой строки матрицы $ f_k(A) $. !!Т!!** Теорема 4.** //Если// $ \lambda_{j} \ne \lambda_k $, то $$ M_j(A)M_k(A) = M_k(A)M_j(A) = \mathbb O_{n\times n} \, . $$ **Доказательство.** $$ M_j(A) = X_j Y_j \, , M_k(A) = X_k Y_k \, . $$ Результат теоремы следует из равенства $ Y_j X_k = 0 $. Оно справедливо поскольку из определений собственных векторов $$ Y_jA=\lambda_j Y_j,\ AX_k=\lambda_k X_k \, , $$ вытекает $$ Y_jAX_k = \left\{\begin{array}{ll} \lambda_j Y_jX_k ;\\ \lambda_k Y_jX_k \end{array} \right. $$ При условии $ \lambda_{j} \ne \lambda_k $ должно выполняться заявленное равенство. !!Т!! **Теорема 5.** //Матрицы// $ M_k(A) $ ((:mapping/operator#osnovnye_opredelenija идемпотентны)), т.е. $$ [M_k(A)]^2=M_k(A) \, . $$ //Кроме того, система матриц// $ \{M_k(A)\}_{k=1}^n $ //линейно независима и// $$ \sum_{k=1}^n M_k(A) =E_{n\times n} \, . $$ Таким образом, представление $$ A=\sum_{k=1}^n \lambda_k M_k(A) $$ можно интерпретировать как разложение матрицы $ A $ по базису $ \{M_k(A)\}_{k=1}^n $, где собственные числа играют роль координат. Причем, ввиду теоремы 4, можно предполагать ортогональность этого базиса... На самом деле, в смысле формального ((#postanovka_zadachi определения)) скалярного произведения $$ \langle M_j(A), M_k(A) \rangle = \operatorname{Sp} (M_j(A) \cdot [M_k(A)]^{\top}) $$ здесь нет ортогональности. Имеется разве что //похожесть на ортогональность//. Теперь проанализируем полученное разложение на предмет ответа на вопрос: какие именно слагаемые вносят в него наибольший вклад? Какие из них наилучшим образом аппроксимируют матрицу $ A $? !!П!! **Пример 1.** $$ A= \left[ \begin {array}{rrrr} -19873&40167&14897&28232 \\ -11466&23174&8594&16288\\ 14522 &-29351&-10886&-20632\\ -5349&10812&4011&7601 \end {array} \right] $$ **Решение.** $ f(\lambda):=\det (A-\lambda E) = {\lambda}^{4}-16\,{\lambda}^{3}-27\,{\lambda}^{2}+178\,\lambda-136 $. Спектр матрицы: $$ \lambda_1 =1,\ \lambda_2=2,\ \lambda_3=17,\ \lambda_4=-4 \, . $$ Не вдаваясь пока в вычислительные детали, приведем выражения для матриц $ M_k(A) $ $$ M_1(A) \approx \left[ \begin {array}{rrrr} 356.3&- 757.4&- 370.3& - 705.6\\ 203.6&- 432.8&- 211.6&- 403.2\\ - 254.5& 541& 264.5& 504\\ 95.4375&- 202.875&- 99.1875&- 189 \end {array} \right] , $$ $$ \ M_2(A) \approx \left[ \begin {array}{rrrr} - 751.66& 1469.16& 587.66& 1239.11\\ - 432.66& 845.66& 338.26& 713.24\\ 539.0&- 1053.5&- 421.4&- 888.53\\ - 198.0& 387.0& 154.8& 326.4 \end {array} \right] , \ \dots $$ Теперь оценим разности $ A-\lambda_k M_k $ для каждого $ k \in \{1,2,3,4\} $ вычислив ((:norm_space#norma_matricy евклидовы нормы)) (нормы Фробениуса): $$ \| A-\lambda_1 M_1 \| \approx 75041.7, \ \| A-\lambda_2 M_2 \| \approx 82542.5,\ \| A-\lambda_3 M_3 \| \approx 6632.4,\ \| A-\lambda_4 M_4 \| \approx 65668.2 $$ Видим, что наилучшее приближение дает слагаемое, соответствующее максимальному собственному числу. А следующим за ним, в плане качества приближения, оказывается слагаемое, соответствующее собственному числу, абсолютная величина (модуль) которого является максимальной из модулей оставшихся собственных чисел. Возникает гипотеза, что за качество приближения отвечают именно модули собственных чисел. Эта гипотеза тут же опровергается следующим по убыванию модуля собственным числом, но... Делаем скидку на то, что $ \lambda_1 $ и $ \lambda_2 $ близки. А если взять пару наилучших слагаемых, то ошибка аппроксимации уменьшается еще значительнeй $$ \| A-\lambda_3 M_3 -\lambda_4 M_4 \| \approx 4348.4 \, . $$ ♦ А сколько таких слагаемых можно оставить в сумме $ \sum_{k=1}^n \lambda_k M_k(A) $ чтобы матрица $ A $ стала распознаваемой? Именно эту задачу мы будем решать в дальнейшем. Поиск решения будет производиться по уже намеченному выше пути. Но придется принять во внимание несколько нетривиальных аспектов. 1. Как правило, возникающие в приложениях матрицы вещественны. И аппроксимировать их требуется также вещественными матрицами. А собственные числа вещественных квадратных матриц не обязательно вещественны (как случилось в предыдущем примере). 2. Собственные числа квадратных матриц не обязательно различны. А все предшествующие рассуждения были справедливы именно в предположении простоты всех собственных чисел. 3. Возникающие в приложениях матрицы не обязательно квадратные. А для неквадратных матриц отсутствует понятие собственного числа. 4. Задача для таких матриц ставится в постановке: аппроксимировать их наилучшим образом матрицами меньшего ранга. Можно ли как-то предсказать ошибку этого наилучшего приближения (без явного построения самого этого приближения)? ===Симметричные матрицы== $$ A= A^{\top} $$ Среди общего случая квадратных матриц, симметричные матрицы выделяются рядом хороших свойств для собственных чисел и собственных векторов. 1. Все собственные числа $ \lambda_1,\dots, \lambda_n $ вещественны ( ☞ ((:algebra2/charpoly/symm доказательство))). 2. Если $ X_k $ --- правый собственный вектор матрицы, принадлежащий собственному числу $ \lambda_k $, то $ X_k^{\top} $ --- левый собственный вектор, принадлежащий тому же собственному числу: $$ A X_k = \lambda_k X_k \quad \Rightarrow \quad (A X_k)^{\top} = (\lambda_k X_k)^{\top} \quad \Leftrightarrow \quad X_k^{\top} A^{\top} = \lambda_k X_k^{\top} \quad \Leftrightarrow \quad X_k^{\top} A = \lambda_k X_k^{\top} \, . $$ 3. Если $ \lambda_j \ne \lambda_k $, то $ X_j^{\top} X_k = 0 $; говорят, что эти собственные векторы **ортогональны** ( ☞ ((algebra2/symmetric#xarakteristicheskij_polinom_sobstvennye_chisla_sobstvennye_vektory доказательство))). 4. Симметричная матрица всегда диагонализуема, т.е. ее ЖНФ диагональна даже в случае наличия кратных собственных чисел ( ☞ ((algebra2/symmetric/vspom1 доказательство))). Более того, матрица, приводящая к диагональному виду может быть выбрана во множестве ортогональных матриц: $$ \exists P \in \mathbb R^{n\times n} \quad \mbox{такая, что} \ P^{\top}P=E,\ P^{\top}A P= \left( \begin{array}{cccc} \lambda_1 & & & \mathbb O \\ & \lambda_2 & & \\ && \ddots & \\ \mathbb O&& & \lambda_n \end{array} \right) \, . $$ !!Т!! **Теорема 6.** //Симметричную матрицу// $ A \ne \mathbb O_{n\times n} $ //всегда можно представить в виде линейной комбинации//: $$ A=\sum_{k=1}^{\operatorname{rank}(A)} \lambda_k M_k(A) $$ //одноранговых симметричных матриц// $$ \{M_k(A):=X_k \cdot X_k^{\top}\}_{k=1}^{\operatorname{rank}(A)} $$ //Здесь// $ X_k $ --- //нормированный собственный вектор//: $$ AX_k=\lambda_k X_k,\ X_k^{\top}X_k=1 \, . $$ // Матрицы// $ \{M_k(A)\}_{k=1}^{\operatorname{rank}(A)} $ //удовлетворяют условиям// $$ \left[ M_k(A) \right]^2= M_k(A),\ \operatorname{Sp} ( M_k(A) ) =1, \ \operatorname{Sp} ( M_j(A) M_k(A) ) = 0 \quad \mbox{при} \ \lambda_j \ne \lambda_k \, . $$ !!П!! **Пример 2.** Представить матрицу $$ A= \left[ \begin {array}{rrr} 0&-1&2\\ -1&3&-1 \\ 2&-1&0\end {array} \right] $$ в виде комбинации матриц из теоремы. **Решение.** Собственные числа и нормированные собственные векторы матрицы: $$ \lambda_1=4, X_{1} = \left( \begin{array}{r} 1/\sqrt{6} \\ -2/\sqrt{6} \\ 1/\sqrt{6} \end{array} \right) \ , \ \lambda_2=-2, X_{2} = \left( \begin{array}{r} -1/\sqrt{2} \\ 0 \\ 1/\sqrt{2} \end{array} \right) \ , \ \lambda_3=1, X_{3}= \left( \begin{array}{r} 1/\sqrt{3} \\ 1/\sqrt{3} \\ 1/\sqrt{3} \end{array} \right) \, . $$ Следовательно $$ A= \lambda_1 X_{1} X_{1}^{\top}+ \lambda_2 X_{2} X_{2}^{\top}+ \lambda_3 X_{3} X_{3}^{\top}= $$ $$ = 4 \, \left[ \begin {array}{rrr} 1/6&-1/3&1/6\\ -1/3&2/3& -1/3\\ 1/6&-1/3&1/6 \end {array} \right] - 2\, \left[ \begin {array}{rrr} 1/2&0&-1/2\\ 0&0&0 \\ -1/2&0&1/2\end {array} \right] + 1\, \left[ \begin {array}{rrr} 1/3&1/3&1/3\\ 1/3&1/3&1/3 \\ 1/3&1/3&1/3 \end {array} \right] \, . $$ ♦ Условия на матрицы $ \{M_k(A)\}_{k=1}^n $ из теоремы 6 означают, что эти матрицы образуют ортонормированный "базис". Почему ортонормированный --- понятно. А вот почему слово базис приведено в кавычках? Можно интерпретировать результат теоремы как разложение симметричной матрицы по ортогональному базису подпространства симметричных матриц? Формально это неверно: количество одноранговых матриц в теореме меньше размерности подпространства симметричных матриц ($=n(n-1)/2$). Тем не менее, для конкретной матрицы $ A $ система матриц $ \{M_k(A)\}_{k=1}^n $ может считаться базисной, а собственные числа $\{\lambda_k\}_{k=1}^n $ --- координатами матрицы в этом ортонормированном базисе. И в этом базисе мы будем искать решение для следующей задачи. **Задача.** Для натурального числа $ \mathfrak r ♦ Количество всевозможных комбинаций из $ n $ собственных чисел по $ \mathfrak r $ равно $ C_n^{\mathfrak r} $. Среди них наименьшее значение для $ || A- A_{\mathfrak r}|| $ достигается когда собственные числа пронумерованы в порядке невозрастания: $$ |\lambda_1| \ge |\lambda_2| \ge \dots \ge |\lambda_{\mathfrak r}| \ge |\lambda_{{\mathfrak r}+1}| \ge \dots \ge |\lambda_n | \, . $$ Оказывается, что соответствующая матрица $ A_{\mathfrak r} $ обеспечивает глобальный минимум в поставленной задаче аппроксимации матрицы $ A$ матрицей ранга $ {\mathfrak r} $. Фактически формула $$ A=\sum_{k=1}^n \lambda_k M_k(A) $$ представляет собой разбиение матрицы $ A $ на слагаемые, расположенные по убыванию их "вклада" в матрицу. Если матрица отражает своими элементами какую-то информационную сущность, Лишь бы только эта матрица не была абсолютно случайной! то можно попробовать заменить ее хранение (пересылку) на матрицу $ A_{\mathfrak r} $ при подходящем выборе $ {\mathfrak r} < \operatorname{rank} (A) $. !!П!! **Пример 3.** Для матрицы примера 2 $$ \|A-\lambda_1M_1(A) \| =\sqrt{5},\ \|A-\lambda_2M_2(A) \|=\sqrt{17},\ \|A-\lambda_3M_3(A) \|= \sqrt{20} \ ; $$ $$ \|A-\lambda_1M_1(A) -\lambda_2M_2(A) \| =1, \ \|A-\lambda_1M_1(A) -\lambda_3M_3(A) \| = 2, \ \|A-\lambda_2M_2(A) -\lambda_3M_3(A) \| = 4 \, . $$ Каков геометрический смысл этой теоремы? ===Метод главных компонент== **Задача.** Найти координаты точки на плоскости, сумма квадратов расстояний до которой от точек $ P_1=(x_{1},y_1),\dots, P_n=(x_n,y_n) $ минимальна. Можно дать этой задаче следующее наглядное толкование: предположим, что на лист бумаги прикрепляется точечная мишень, в которую стрелок стреляет из пистолета. После стрельбы лист бумаги с $ n $ пулевыми отверстиями, но со снятой мишенью, передается эксперту с требованием установить вероятное положение мишени. Предположения, в которых эта задача решается, стандартные: пистолет пристрелян и не дает сносов, а сам стрелок не страдает косоглазием. !!Т!! **Теорема 1.** //Координаты искомой точки определяются средними значениями координат точек// $ \{ P_j\}_{j=1}^n $: $$ \overline{x}=\frac{x_1+x_2+ \dots+x_n}{n},\quad \overline{y}=\frac{y_1+y_2+ \dots+y_n}{n}\ . $$ Обозначение среднего значения числовой последовательности через $ \overline{\color{White}{x}} $ является традиционным; не следует его путать с комплексным сопряжением. Точку $ (\overline{x},\overline{y}) $ иногда называют **центроидом** для системы точек $ \{ P_j\}_{j=1}^n $. !!=>!! Координаты точки, для которой величина $$ \sum_{j=1}^n m_j \left[(x-x_j)^2+(y-y_j)^2 \right] \ , $$ ($ m_{1},\dots,m_n $ --- положительные константы) будет минимальной: $$ \hat x=\frac{m_1x_1+m_2x_2+ \dots+m_nx_n}{m_1+m_2+ \dots+m_n},\quad \hat y=\frac{m_1y_1+m_2y_2+ \dots+m_ny_n}{m_1+m_2+ \dots+m_n}\ . $$ Иными словами: искомая точка совпадает с **центром масс** (**барицентром**) системы масс $ \{m_j\}_{j=1}^n $, расположенных в точках $ \{P_j\}_{j=1}^n $. **Задача.** Найти прямую на плоскости, сумма квадратов расстояний до которой от точек $ \{P_j\}_{j=1}^n $ минимальна. В постановке задачи стрельбы по мишени, источником этой задачи можно считать стрельбу по прямолинейному отрезку (скажем, по участку дороги). Требуется установить наиболее правдоподобное расположение дороги исходя из результатов попаданий. Если обозначить $ \delta_j $ расстояние от $ P_j $ до прямой $$ ax+by+c=0 \, , $$ то речь о поиске $$ \min_{(a,b,c)\in \mathbb R^3} (\delta_1^2+\delta_2^2+\dots + \delta_n^2) \, . $$ {{ pearson2.png |}} !!Т!! **Теорема 2 [Пирсон].** //Искомый минимум равен минимальному положительному корню// $ \lambda_{min} $ //полинома// $$ f(\lambda):=\left| \begin{array}{cc} \displaystyle \sum_{j=1}^n (x_j-\overline{x})^2 - \lambda & \displaystyle \sum_{j=1}^n (x_j-\overline{x})(y_j-\overline{y}) \\ \displaystyle \sum_{j=1}^n (x_j-\overline{x})(y_j-\overline{y}) & \displaystyle \sum_{j=1}^n (y_j-\overline{y})^2 - \lambda \end{array} \right| \, . $$ **Доказательство** ☞ ((:interpolation:pearson ЗДЕСЬ)). Корни квадратного полинома $ f(\lambda) $ из теоремы Пирсона вещественны и неотрицательны. Можно доказать это применением ((:algebra/dets/binet_cauchy неравенства Коши-Буняковского)). С другой стороны, для доказательства можно применить и более общий теоретический результат, который позволит распространить результат теоремы в $\mathbb R^3 $ (и даже в пространства бóльших размерностей). Дело в том, что полином $ f(\lambda) $ является ((:algebra2/charpoly#xarakteristicheskij_polinom характеристическим полиномом матрицы)) $$ G=\left(\begin{array}{cc} \displaystyle \sum_{j=1}^n (x_j-\overline{x})^2 & \displaystyle \sum_{j=1}^n (x_j-\overline{x})(y_j-\overline{y}) \\ \displaystyle \sum_{j=1}^n (x_j-\overline{x})(y_j-\overline{y}) & \displaystyle \sum_{j=1}^n (y_j-\overline{y})^2 \end{array} \right) \, , $$ которую можно интерпретировать как ((:dets:gram матрицу Грама)) системы векторов $$ \left\{ (x_1-\overline{x},\dots, x_n-\overline{x}), (y_1-\overline{y},\dots, y_n-\overline{y}) \right\} \, . $$ при ((:euclid_space#opredelenija стандартном способе задания скалярного произведения)) в $\mathbb R^n $. Доказывается (см. ☞ ((:dets/gram#sobstvennye_chisla_sobstvennye_vektory ЗДЕСЬ))), что корни характеристического полинома матрицы Грама (т.е. ее ((:algebra2/charpoly#sobstvennoe_chislo собственные числа))) всегда вещественны и неотрицательны. !!=>!! Если корни полинома $ f(\lambda) $ различны, то прямая, минимизирующая сумму квадратов расстояний до нее от точек $ \{P_j\} $ единственна, она проходит через точку $ (\overline{x},\overline{y}) $, и ее направляющим вектором является ((:algebra2:charpoly#sobstvennyj_vektor собственный вектор матрицы)) $ G $, принадлежащий __максимальному__ собственному числу этой матрицы. По аналогии хочется утверждать, что собственный вектор, соответствующий числу $ \lambda_{min} $ будет направляющим вектором прямой, обеспечивающей __максимум__ сумме квадратов. Однако, это утверждение неверно. Оно будет верно, если мы добавим условие "...среди всех прямых, проходящих через точку $ (\overline{x},\overline{y}) $. Обе прямые, порождаемые этими "оптимизационными" задачами, перпендикулярны. !!П!! **Пример [Пирсон].** Построить прямую из теоремы $ 2 $ для точек $$ \begin{array}{c|cccccccccc} x & 0.0 & 0.9 & 1.8 & 2.6 & 3.3 & 4.4 & 5.2 & 6.1 & 6.5 & 7.4 \\ \hline y & 5.9 & 5.4 & 4.4 & 4.6 & 3.5 & 3.7 & 2.8 & 2.8 & 2.4 & 1.5 \end{array} $$ **Решение.** Здесь $ n=10 $ и $$ \sum_{j=1}^{10} x_j= 38.2,\ \sum_{j=1}^{10} y_j= 37.0 \quad \Rightarrow \quad \overline x =3.82, \overline y = 3.70 \, . $$ $$ G=\left[ \begin {array}{rr} 56.396&- 30.430 \\ - 30.430& 17.220 \end {array} \right] $$ Характеристический полином $$ f(\lambda)=\lambda^2-73.616\,\lambda+45.15422 $$ имеет корни $$ \lambda_1 \approx 0.618573,\ \lambda_2 \approx 72.997427 \, . $$ Собственные векторы, принадлежащие этим собственным числам: $$ \approx \left[ 0.478924, 0.877856\right]^{\top} \quad \mbox{и} \quad \approx \left[ 0.877856, - 0.478924\right]^{\top} $$ соответственно; оба нормированы. Сумма квадратов расстояний от заданных точек до ближайшей к ним прямой равна $ \lambda_1 $. Из следствия к теореме Пирсона получаем уравнение этой прямой $$ 0.478924\, x+0.877856\,y-5.077556=0 \quad \Leftrightarrow \quad y = -0.545561\, x + 5.784082 \, . $$ {{ :pearson1.png?400 |}} Эта прямая совпадает с решением Пирсона, но вот значение для минимума у него приведено другое: $ 0.2484 $. Непосредственная проверка подтверждает значение $ \lambda_1 $. ♦ В статистике матрица $$ \frac{1}{n} G $$ составленная для векторов $ \{X,Y\} \subset \mathbb R^n $ экспериментальных данных, называется **ковариационной матрицей** этих векторов. В примере Пирсона исходные данные взяты абстрактными, не имеющими прототипа в реальном мире. ==Вычислительные аспекты== Для симметричных матриц больших порядков проблема вычисления полного разложения в сумму одноранговых симметричных матриц $$ A=\sum_{k=1}^n \lambda_k M_k(A) \quad \mbox{при} \ M_k(A)=X_k X_k^{\top}, \ X_k^{\top} X_k =1 $$ требует нахождения ее спектра, т.е. набора всех ее собственных чисел, а также соответствующих собственных векторов. Однако уже проблема вычисления характеристического полинома матриц порядка $ 100 $ является неприятной. Теоретически ((algebra2/charpoly#metody_vychislenija_xarakteristicheskogo_polinoma алгоритмы)) существуют, но на практике их не используют. Применяют подход, решающий ((:algebra2/charpoly#chastichnaja_problema_sobstvennyx_chisel частичную проблему собственных чисел)), т.е. проблему ограничивают вычислением максимальных по модулю собственных чисел и соответствующих собственных векторов. Пусть собственные числа симметричной матрицы $ A \in \mathbb R^{n\times n} $ занумерованы в порядке убывания модулей: $$ |\lambda_1|>|\lambda_2|> \dots > |\lambda_n | $$ и среди этих модулей нет одинаковых. Тогда $ \lambda_1 $ может быть вычислен ((:algebra2/charpoly#chastichnaja_problema_sobstvennyx_chisel методом степенных итераций)). Находим его с необходимой точностью (пока не известно какой), а также вычисляем соответствующий ему (правый) собственный вектор $ X_1 \in \mathbb R^{n} $ (столбец), нормируем его $ \| X_1 \|=1 $. И выбрасываем из разложения матрицы $ A $ первое слагаемое: $$ A_2:=A- \lambda_1 X_1 \cdot X_1^{\top} = \sum_{k=2}^n \lambda_k M_k(A) $$ Если все вычисления производятся точно, то $ \operatorname{rank} (A_2)= \operatorname{rank} (A)-1 $, спектром матрицы $ A_2 $ является набор $ \{\lambda_2,\lambda_3, \dots, \lambda_n, 0 \} $, а ее доминирующим (по модулю) собственным числом становится $ \lambda_2 $. Для его поиска снова применяем метод степенных итераций. ==Сингулярное разложение произвольной матрицы== ===Некоторые предварительные сведения== о матрицах $ A\cdot A^{\top} $ и $ A^{\top} A $ для матрицы $ A\in \mathbb R^{m\times n} $. 1. Обе матрицы --- симметричны (хотя и не обязательно одинаковых порядков). 2. Матрица $ A\cdot A^{\top} $ --- ((:dets:gram#матрица_и_определитель_грама матрица Грама)) системы строк, а $ A^{\top} A $ --- матрица Грама системы столбцов матрицы $ A $. Так если представить матрицу $ A $ как результат конкатенации ее столбцов $ A=[A_{[1]}|A_{[2]}|\dots |A_{[n]} ] $, то $$ A^{\top} A = G(A_{[1]},A_{[2]},\dots ,A_{[n]}) \, ; $$ здесь скалярное произведение задано стандартным способом. Отсюда, в частности, следует: $$ \operatorname{rank} ( A\cdot A^{\top}) = \operatorname{rank} (A^{\top} A) = \operatorname{rank}(A) \, . $$ 3. ((:algebra2:charpoly Характеристические полиномы)) этих матриц могут различаться только лишь на степень переменной: $$ (-\lambda)^n \det (A \cdot A^{\top} - \lambda E_{m\times m})\equiv (-\lambda)^m \det ( A^{\top} A - \lambda E_{n\times n}) \ . $$ 4. Любое ненулевое ((:algebra2:charpoly#собственное_число собственное число)) одной матрицы вещественно, положительно и является собственным числом другой матрицы. Количество ненулевых собственных чисел совпадает с $ \operatorname{rank} (A) $. 5. Если обозначить $ \{ \sigma_1^2,\dots, \sigma_{\operatorname{rank} (A)}^2 \} $ ненулевые собственные числа матрицы $ A\cdot A^{\top} $ (или $ A^{\top} A $), а соответствующие им ортонормированные собственные векторы $ \{ U_{[1]},\dots U_{[\operatorname{rank} (A)]} \} \subset \mathbb R^m $ (соответственно, $ \{ V_{[1]},\dots V_{[\operatorname{rank} (A)]} \} \subset \mathbb R^n $), то справедливы равенства $$ A\cdot A^{\top} = \sum_{k=1}^{\operatorname{rank} (A)} \sigma_k^2 U_{[k]}U_{[k]}^{\top} \quad \mbox{и} \quad A^{\top} A = \sum_{k=1}^{\operatorname{rank} (A)} \sigma_k^2 V_{[k]}V_{[k]}^{\top} \, . $$ 6. Существуют ((:algebra2:funmatrix#квадратный_корень_из_матрицы положительно полуопределенные квадратные корни)) из матриц $ A\cdot A^{\top} $ и $ A^{\top} A $, т.е. вещественные симметричные матрицы $ P \in \mathbb R^{m\times m} $ и $ Q \in \mathbb R^{n\times n}$ такие, что $$ P^2= A\cdot A^{\top},\ Q^2=A^{\top} A \, . $$ ===Сингулярное разложение== !!Т!! **Теорема ??.** //Матрицу// $ A \in \mathbb R^{m\times n}, A\ne \mathbb O_{m\times n} $ //всегда можно представить в виде линейной комбинации//: $$ A=\sum_{k=1}^{\operatorname{rank}(A)} \sigma_k M_k(A) \qquad \qquad \qquad {\color{Red}{ (1)} } $$ //одноранговых матриц// $$ \left\{M_k(A):=U_{[k]} V_{[k]}^{\top}\right\}_{k=1}^{\operatorname{rank}(A)} \, . $$ //Здесь// $ \{\sigma_k>0\} $ --- //арифметические корни квадратные из ненулевых собственных чисел матрицы// $ A \cdot A^{\top} $ . //Система столбцов// $ \{U_{[k]}\} \subset \mathbb R^{m} $ //состоит из ортонормированных собственных векторов, принадлежащих// $ \{\sigma_k^2\} $, //а система столбцов// $ \{V_{[k]}\} \subset \mathbb R^{n} $ //состоит из ортонормированных собственных векторов матрицы// $ A^{\top} A $,// принадлежащих// $ \{\sigma_k^2\} $. //При этом столбцы// $ U_{[k]} $ и $ V_{[k]} $ //связаны дополнительными соотношениями// $$ A V_{[k]} = \sigma_k U_{[k]}, \ A^{\top} U_{[k]} = \sigma_k V_{[k]} \, \qquad \qquad \qquad {\color{Red}{ (2)} } . $$ Числа $ \{ \sigma_k \} $ называются **сингулярными числами матрицы** $ A $; столбец $ U_{[k]} $ называют **левым сингулярным вектором**, а столбец $ V_{[k]} $ --- **правым сингулярным вектором** матрицы $ A $, принадлежащими сингулярному числу $ \sigma_k $. Представление матрицы $ A $ по формуле $ {\color{Red}{ (1)} } $ называется **сингулярным разложением**[[(//англ.//) Singular Value Decomposition, SVD]] матрицы $ A $. Единственности сингулярного разложения не гарантируется. Условия из теоремы для определения левых и правых собственных векторов избыточны. Действительно, вычислив, к примеру, правый сингулярный вектор $ V_{[1]} $ как __нормированный__ собственный вектор матрицы $ A^{\top} A $, принадлежащий $ \sigma_1 $, мы немедленно найдем и левый сингулярный вектор $ U_{[1]} $ по первой из формул $ {\color{Red}{ (2)} } $: $ U_{[1]}= \frac{1}{\sigma_1} A V_{[1]}$. И он действительно будет собственным вектором матрицы $ A\cdot A^{\top} $: $$ A\cdot A^{\top} U_{[1]} = \frac{1}{\sigma_1} A\cdot A^{\top} A V_{[1]} = \frac{1}{\sigma_1} A (\sigma_1 V_{[1]}) = \sigma_1 U_{[1]} \, . $$ Кроме того, $$ U_{[1]}^{\top} U_{[1]}=\frac{1}{\sigma_1^2} V_{[1]}^{\top} A^{\top} A V_{[1]} = \frac{1}{\sigma_1^2} V_{[1]}^{\top} (\sigma_1^2 V_{[1]})=V_{[1]}^{\top} V_{[1]}=1 \ , $$ т.е. если $ V_{[1]} $ --- нормированный, то и $ U_{[1]} $ --- нормированный. Имеет смысл начинать вычисления с тех сингулярных векторов, которые являются собственными для матрицы наименьшего порядка из $ m\times m $ или $ n \times n $ (да и собственные числа в таком варианте вычислять проще!). !!П!! **Пример 1.** Построить сингулярное разложение матрицы $$ A=\left( \begin{array}{rrr} 1 & -1 & -2 \\ -7/3 & 1/3 & 2/3\\ 1/3 & -7/3 & -2/3 \\ -5/3 & 5/3 & -2/3 \end{array} \right) \, . $$ **Решение.** Из двух вариантов $ A \cdot A^{\top} $ и $ A^{\top} A $ выбираем второй: $$ A^{\top} A = \left( \begin{array}{rrr} 28/3 & -16/3 & -8/3 \\ -16/3 & 28/3 & 8/3\\ -8/3 & 8/3 & 16/3 \end{array} \right) \, . $$ Характеристический полином $$ \det (A^{\top} A - \lambda E)= -(\lambda-16)(\lambda-4)^2 \, . $$ Следовательно сингулярными числами матрицы $ A $ являются $ \sigma_1=4,\sigma_2=2,\sigma_3=2 $. Ищем ортонормированные собственные векторы матрицы $ A^{\top} A $. Если бы все собственные числа матрицы $ A^{\top}A $ были простыми, то можно было бы воспользоваться универсальным представлением собственного вектора посредством матрицы, взаимной матрице $ A \cdot A^{\top} - \lambda E $ (см. ((/algebra2/charpoly#sobstvennyj_vektor ЗДЕСЬ))). Достаточно найти ее первый столбец: $$ \left( \begin{array}{crr} (\lambda-32/3)(\lambda-4) & * & * \\ 16/3(4-\lambda) & * & * \\ 8/3(4-\lambda) & * & * \end{array} \right) $$ Подстановка $ \lambda=16 $ в него дает собственный вектор в виде $ [64,-64,-32]^{\top} $; после нормирования получим $ V_{[1]}=[2/3,-2/3,-1/3]^{\top} $ . Однако при значении кратного корня $ \lambda = 4 $ столбец вырождается в нулевой, т.е. поиск двух линейно независимых собственных векторов посредством взаимной матрицы принципиально неосуществим. На помощь приходит ((:mapping:operator:jordan#аннулирующий_полином минимальный аннулирующий полином)) матрицы $ A^{\top} A $. Таким является $ (\lambda-16)(\lambda-4) $. Матрица $$ A^{\top} A - 16 E= \left( \begin{array}{rrr} -20/3 & -16/3 & -8/3 \\ -16/3 & -20/3 & 8/3 \\ -8/3 & 8/3 & -32/3 \end{array} \right) $$ имеет ранг равный $ 2 $, и любая пара ее столбцов будет линейно независимой системой собственных векторов матрицы, принадлежащих собственному числу $ \lambda=4 $. Возьмем, например, второй и третий столбцы и проведем с ними процедуру ((:euclid_space#ортогонализация ортогонализации Грама-Шмидта)). Получим столбцы $ [-16/3,-20/3, 8/3]^{\top} $, $ [-12,-0,-24]^{\top} $, которые после нормирования дают $$ V_{[2]}=\left[ -4/(3\sqrt{5}),-\sqrt{5}/3, 2/\sqrt{5} \right]^{\top}, V_{[3]}=\left[1/\sqrt{5},0, 2/\sqrt{5} \right]^{\top} \, . $$ Нахождение левых сингулярных векторов производится с помощью формулы $ {\color{Red}{ (2)} } $ как $ U_{[j]}= 1/\sigma_j A V_{[j]} $: $$ U_{[1]}=\left[\frac{1}{2},-\frac{1}{2},\frac{1}{2},-\frac{1}{2}\right]^{\top},\ U_{[2]}=\left[-\frac{1}{2\sqrt{5}},\frac{3}{2\sqrt{5}},\frac{3}{2\sqrt{5}},-\frac{1}{2\sqrt{5}}\right]^{\top},\ $$ $$ U_{[3]}=\left[-\frac{3}{2\sqrt{5}},-\frac{1}{2\sqrt{5}},-\frac{1}{2\sqrt{5}},-\frac{3}{2\sqrt{5}}\right]^{\top} \, . $$ Формула $ {\color{Red}{ (1)} } $ дает представление матрицы $ A $ в виде $$ 4 \left( \begin{array}{rrr} 1/3 & -1/3 & -1/6 \\ -1/3 & 1/3 & 1/6 \\ 1/3 & -1/3 & -1/6 \\ -1/3 & 1/3 & 1/6 \end{array} \right)+2 \left( \begin{array}{rrr} -1/3 & -1/6 & -1/3 \\ -1/3 & -1/6 & -1/3 \\ -1/3 & -1/6 & -1/3 \\ -1/3 & -1/6 & -1/3 \end{array} \right) + 2 \left( \begin{array}{rrr} 1/6 & 1/3 & -1/3 \\ -1/6 & -1/3 & 1/3 \\ -1/6 & -1/3 & 1/3 \\ 1/6 & 1/3 & -1/3 \end{array} \right) \, . $$ Впрочем, на практике экономнее хранить не матрицы с фактически одинаковыми строками (столбцами), а именно их представления через $ U_{[j]} $ и $ V_{[j]} $. !!=>!! Система матриц $\{M_k(A)\} $ ортонормирована. $$ \operatorname{Sp} \left( (U_{[j]} V_{[j]}^{\top})^{\top} U_{[k]} V_{[k]}^{\top} \right)= \operatorname{Sp} \left( V_{[j]}U_{[j]}^{\top} U_{[k]} V_{[k]}^{\top} \right) = 0 $$ поскольку $ U_{[j]}^{\top} U_{[k]} =0 $. ====Представление сингулярного разложения в виде произведения матриц== Часто формула сингулярного разложения матрицы записывается в форме матричного произведения. !!Т!! **Теорема.** //Составим матрицы из сингулярных чисел и сингулярных векторов матрицы// $ A\in \mathbb R^{m\times n} $: $$ {\color{Red}{\mathbf{\Sigma}}}:= \operatorname{diag}(\sigma_1,\dots, \sigma_{\operatorname{rank} (A)}) = \left( \begin{array}{cccc} \sigma_1 & 0 & \dots & 0 \\ 0 & \sigma_2 & \dots & 0 \\ \vdots && \ddots & \vdots \\ 0 & 0 & \dots & \sigma_{\operatorname{rank} (A)} \end{array} \right) $$ $$ U:=\left[U_{[1]}\mid \dots \mid U_{[\operatorname{rank} (A)]} \right], \ V:=\left[V_{[1]}\mid \dots \mid V_{[\operatorname{rank} (A)]} \right] \ , $$ //Справедливо равенство// $$ A=U \cdot {\color{Red}{\mathbf{\Sigma}}} \cdot V^{\top} \qquad \qquad {\color{Red}{ (1)^{\prime}} }\, . $$ Эквивалентность этого представления формуле $ {\color{Red}{ (1)} } $ вытекает из следующих соображений: $$ {\color{Red}{\mathbf{\Sigma}}}= \operatorname{diag}(\sigma_1, 0,\dots, 0) + \operatorname{diag}(0, \sigma_2,\dots, 0)+\dots $$ $$ U {\color{Red}{\mathbf{\Sigma}}} V^{\top}=U(\operatorname{diag}(\sigma_1, 0,\dots, 0) + \operatorname{diag}(0, \sigma_2,\dots, 0)+\dots) V^{\top} = $$ $$ =U \operatorname{diag}(\sigma_1, 0,\dots, 0) V^{\top} + \dots = \sigma_1 \left[U_{[1]}\mid \dots \mid U_{[\operatorname{rank} (A)]} \right] \operatorname{diag}(1, 0,\dots, 0) V^{\top} + \dots = $$ $$ = \sigma_1 \left[U_{[1]}\mid \mathbb O \mid \dots \mid \mathbb O \right] V^{\top} + \dots = \sigma_1 \left[U_{[1]}\mid \mathbb O \mid \dots \mid \mathbb O \right] \left[ \begin{array}{c} V_{[1]}^{\top} \\ V_{[2]}^{\top} \\ \vdots \end{array} \right] + \dots = \sigma_1 U_{[1]} V_{[1]}^{\top} + \dots $$ То есть мы получили первое слагаемое формулы $ {\color{Red}{ (1)} } $. Для остальных слагаемых рассуждения аналогичны. Традиционно диагональная матрица в произведении $ {\color{Red}{ (1)^{\prime} }} $ обозначается $ {\mathbf{\Sigma}} $. И это приводит к коллизии с обозначением ((:interpolation/pearson#kovariacionnaja_matrica ковариационной матрицы))[[Не говоря уж об обозначении суммирования!]]. Поэтому в настоящем разделе будем использовать красный цвет: ${\color{Red}{\mathbf{\Sigma}}} $. !!=>!! Для симметричной ((:2form#znakoopredelennost положительно полуопределенной матрицы)) $ A\in \mathbb R^{n\times n} $ формула сингулярного разложения эквивалентна ((:algebra2/symmetric/vspom1 формуле приведения этой матрицы к диагональному виду с помощью ортогонального преобразования)). Действительно, в качестве сингулярных чисел выступают ненулевые собственные числа матрицы $ A $. Матрицы $ U $ и $ V $ одинаковы, и состоят из ортонормированной системы собственных векторов, принадлежащих этим собственным числам. В частности, при $ \operatorname{rank}(A)=n $ матрица $ U=V $ ортогональна: $ U^{\top} U = E $. ===Аппроксимация матриц== ====Симметричная матрица== ====Произвольная матрица== !!Т!! **Теорема ??.** Упорядочим сингулярные числа матрицы $ A\in \mathbb R^{m\times n} $ по невозрастанию: $$ \sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_{\operatorname{rank}(A)} >0 \, . $$ //Матрицей// $ \widetilde A \in \mathbb R^{m\times n} $ ранга $ \mathfrak r , 0<\mathfrak r \le \operatorname{rank}(A) $, //наилучшим образом приближающей// $ A $ //в ((:norm_space#норма_матрицы евклидовой норме (норме Фробениуса)))//: $$ \min_{\widetilde A\in \mathbb R^{m\times n}, \operatorname{rank} (\widetilde A)= \mathfrak r} \| A- \widetilde A\| \, , $$ //является матрица// $$ A_{\mathfrak r}=\sum_{j=1}^{\mathfrak r} \sigma_j U_{[j]} V_{[j]}^{\top} \, . $$ //При такой матрице// $ \widetilde A $ //величина указанного минимума равна// $$ \sqrt{\sum_{j=\mathfrak r+1}^{\operatorname{rank}(A)} \sigma_j^2} \, . $$ **Доказательство** ((:algebra2:sumrankone:approx ЗДЕСЬ)) !!=>!! Обозначим $$ U_{\mathfrak r}:= [U_{[1]}|\dots | U_{[\mathfrak r]}] \in \mathbb R^{n\times \mathfrak r } \, .$$ Тогда матрица $ A_{\mathfrak r} $ из теоремы представима в виде $$ A_{\mathfrak r} = U_{\mathfrak r} \cdot U_{\mathfrak r}^{\top} A \, . $$ Аналогично, если $$ V_{\mathfrak r}:= [V_{[1]}|\dots | V_{[\mathfrak r]}] \in \mathbb R^{\mathfrak r \times m } \, .$$ то $$ A_{\mathfrak r} = A V_{\mathfrak r} \cdot V_{\mathfrak r}^{\top} \, . $$ ==Вычислительные аспекты==