<jupyter_start><jupyter_text>Лекция 4 Часть 1. Умножение матрицы на вектор и умножение матриц, иерархия памяти и алгоритм Штрассена  На прошлой лекции- Линейные системы- <NAME> и число обусловленности- Линейная задача наименьших квадратов  Примеры пиковой производительности**Flops** –– число операций над числами с плавающей точкой в секунду (floating point operations per second).Гига = $2^{30} \approx 10^9$,  Тера = $2^{40} \approx 10^{12}$,      Пета = $2^{50} \approx 10^{15}$,       Экза = $2^{60} \approx 10^{18}$ Какова **пиковая производительность**:1. Cовременных CPU?2. Cовремененных GPU?3. Cамого большого суперкомпьютера?   Тактовая частота CPU vs. производительность в flopsFLOPS = sockets * (cores per socket) * (number of clock cycles per second) * (number of floating point operations per cycle).- Обычно sockets = 1- Число ядер обычно 2 или 4- Число тактов в секунду - это привычная частота- Число операций за один такт зависит от конкретного процессора 1. Современных CPU (Intel Core i7)   –– 400 Gflops2. Современных GPU (NVidia Titan X) –– 12 Tflops single precision 3. [Самый большой суперкомпьютер](https://www.top500.org/lists/top500/2020/06/) –– 513.85 Pflops   Умножение матрицы на вектор (matvec)Умножение матрицы $A$ размера $n\times n$ на вектор $x$ размера $n\times 1$ ($y=Ax$):$$y_{i} = \sum_{i=1}^n a_{ij} x_j$$требует $n^2$ умножений и $n(n-1)$ сложений. Таким образом, общая сложность равна $2n^2 - n =$  $\mathcal{O}(n^2)$   Насколько плоха сложность $\mathcal{O}(n^2)$?- Пусть $A$ – матрица попарных гравитационных взаимодействий между планетами в галактике.- Число планет в средней галактике $10^{11}$, поэтому размер этой матрицы $10^{11} \times 10^{11}$.- Для моделирования эволюции во времени нам нужно умножать эту матрицу на вектор для получения состояния в каждый момент времени.- Лучшие суперкомпьютеры делают примерно $10^{16}$ операций в секунду (flops), поэтому для умножения матрицы на вектор нам потребуется примерно$$\frac{(10^{11})^2 \text{ операций}}{10^{16} \text{ flops}} = 10^6 \text{ сек.} \approx 11.5 \text{ дней} $$для одного шага по времени. Если бы мы могли умножать матрицу на вектор за $\mathcal{O}(n)$ операций, мы бы получили$$\frac{10^{11} \text{ операций}}{10^{16} \text{ flops}} = 10^{-5} \text{ сек}.$$Далее посмотрим видео, моделирующее столкновение двух галактик, которое было получено с использованием алгоритма умножения матрицы на вектор со сложностью $\mathcal{O}(n \log n)$:<jupyter_code>from IPython.display import YouTubeVideo
YouTubeVideo("7HF5Oy8IMoM")<jupyter_output><empty_output><jupyter_text>Можно ли умножить матрицу на вектор быстрее $\mathcal{O}(n^2)$?- Вообще говоря **НЕТ**. - Это связано с тем, что мы имеем $\mathcal{O}(n^2)$ входных данных, поэтому для произвольной матрицы быстрее выполнить умножение нельзя.- К счастью, мы можем  умножать матрицу на вектор быстрее для определённого класса матриц.Например:    * Малоранговые матрицы (или матрицы с малоранговыми блоками)    * Разреженные матрицы с $\mathcal{O}(n)$ ненулевыми элементами    * Структурированные матрицы:      * Фурье      * Циркулянтные      * Тёплицевы       * Ганкелевы  Произведение матрицРассмотрим композицию линейных операторов:1. $y = Bx$2. $z = Ay$Тогда, $z = Ay =  A B x = C x$, где $C$ – **произведение двух матриц**.  Классическое определение**Определение**. Произведение матрицы $A$ размера $n \times k$ и матрицы $B$ размера $k \times m$ – это матрица $C$ размера $n \times m$ такая что её элементы записываются как $$   c_{ij} = \sum_{s=1}^k a_{is} b_{sj}, \quad i = 1, \ldots, n, \quad j = 1, \ldots, m $$Для $m=k=n$ сложность наивного алгоритма составляет $2n^3 - n^2 =$ $\mathcal{O}(n^3)$.   Обсуждение- Умножение матриц – это основная операция, необходимая для почти всех эффективных алгоритмов в вычислительной линейной алгебре.- Говоря нестрого, все алгоритмы для решения задач вычислительной линейной алгебры с плотными матрицами сводятся к последовательности умножений матриц.- Повышение эффективности умножения матриц, пропорционально увеличивает эффективность алгоритмов, которые используют эту операцию (то есть почти всех).- Однако реализация быстрых методов для умножения матриц – непростая задача!  Эффективная реализация алгоритмов вычисления произведения матриц**Вопрос**: Легко ли вычислить произведение матриц наиболее эфективным образом?  Ответ: нет, это не легко!Если вы хотите вычислить произведение как можно быстрее, используя обычный компьютер.  ДемоСравним время вычисления произведения матриц для следующих реализаций- стандартная функция `np.dot()`- рукописная реализация алгоритма - Numba версия рукописной реализации<jupyter_code>import numpy as np
def matmul(a, b):
    n = a.shape[0]
    k = a.shape[1]
    m = b.shape[1]  
    c = np.zeros((n, m))
    for i in range(n):
        for j in range(m):
            for s in range(k):
                c[i, j] += a[i, s] * b[s, j]
                
    return c
import numpy as np
from numba import jit # Just-in-time compiler for Python, see http://numba.pydata.org 

@jit(nopython=True)
def numba_matmul(a, b):
    n = a.shape[0]
    k = a.shape[1]
    m = b.shape[1]
    c = np.zeros((n, m))
    for i in range(n):
        for j in range(m):
            for s in range(k):
                c[i, j] += a[i, s] * b[s, j]
    return c<jupyter_output><empty_output><jupyter_text>Теперь давайте сравним время вычислений.Попробуйте угадать ответ :)<jupyter_code>n = 100
a = np.random.randn(n, n)
b = np.random.randn(n, n)

%timeit matmul(a, b)
%timeit numba_matmul(a, b)
%timeit np.dot(a, b)<jupyter_output>536 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.02 ms ± 12.7 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
18.4 µs ± 852 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)<jupyter_text>Для всех ли размерностей выполнены такие соотношения?<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

dim_range = [10*i for i in range(1, 11)]
time_range_matmul = []
time_range_numba_matmul = []
time_range_np = []
for n in dim_range:
    print("Dimension = {}".format(n))
    a = np.random.randn(n, n)
    b = np.random.randn(n, n)

    t = %timeit -o -q matmul(a, b)
    time_range_matmul.append(t.average)
    t = %timeit -o -q numba_matmul(a, b)
    time_range_numba_matmul.append(t.average)
    t = %timeit -o -q np.dot(a, b)
    time_range_np.append(t.average)
plt.plot(dim_range, time_range_matmul, label="Matmul")
plt.plot(dim_range, time_range_numba_matmul, label="Matmul Numba")
plt.plot(dim_range, time_range_np, label="Numpy")
plt.legend(fontsize=18)
plt.xlabel("Dimension", fontsize=18)
plt.ylabel("Time", fontsize=18)
plt.yscale("log")<jupyter_output><empty_output><jupyter_text>Почему наивная реализация такая медленная?Это так по двум причинам:- Она не использует преимуществ быстрой памяти (кэшей) и в целом архитектуру памяти- Она не использует параллелизм (особенно важно для GPU)  Архитектура памяти  - "Быстрой" памяти мало- "Медленной" памяти много  Реализация алгоритмов наиболее вычислительно интенсивным образом**Эффективные реализации для задач NLA** используют блочные версии алгоритмов. Эта идея лежит в основе **[BLAS (Basic Linear Algebra Subroutines)](http://www.netlib.org/blas/)**, написанной на Фортране много лет назад, однако до сих пор не потерявшей актуальность.  Разделим матрицу на блоки! Для примера рассмотрим разбиение на блочную матрицу размера $2 \times 2$:$$   A = \begin{bmatrix}         A_{11} & A_{12} \\         A_{21} & A_{22}        \end{bmatrix}, \quad B = \begin{bmatrix}         B_{11} & B_{12} \\         B_{21} & B_{22}        \end{bmatrix}$$Тогда  $$AB = \begin{bmatrix}A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22} \\            A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22}\end{bmatrix}.$$   Если $A_{11}, B_{11}$ и их произведение помещается в кэш-память (которая достигает 12 Мб (L3) в [последнем CPU от Intel](https://en.wikipedia.org/wiki/List_of_Intel_microprocessorsDesktop)), тогда мы только один раз загружаем их в память.    BLASBLAS имеет три уровня:1. BLAS-1,  операции типа $c = a + b$2. BLAS-2, операции типа умножения матрицы на вектор3. BLAS-3, произведение матриц**Q:** Какова принципиальная разница между ними? **A:** основная разница в соотношении между количеством операций и количеством данным!1. BLAS-1: $\mathcal{O}(n)$ данных, $\mathcal{O}(n)$ операций2. BLAS-2: $\mathcal{O}(n^2)$ данных, $\mathcal{O}(n^2)$ операций3. BLAS-3: $\mathcal{O}(n^2)$ данных, $\mathcal{O}(n^3)$ операций  Почему BLAS так важна и актуальна?1. Стандартные реализации базовых операций линейной алгебры2. Сохраняются стандартные имена для всех операций, сохраняющиеся в каждой новой реализации (например [ATLAS](https://www.netlib.org/atlas/), [OpenBLAS](https://www.openblas.net/), [MKL](https://software.intel.com/en-us/mkl)). Вы можете вызвать функцию умножения матриц (GEMM) в вашем коде, слинковать его с любой реализацией BLAS, и всё будет корректнро работать3. Позволяет формулировать новые алгоритмы в терминах BLAS операций – гарантии скорости4. Есть обёртки для большинства популярных языков  Реализации BLAS и их использование1. [ATLAS](http://math-atlas.sourceforge.net) - Automatic Tuned Linear Algebra Software. Она автоматически адаптируется под архитектуру конкретного компьютера.2. [LAPACK](http://www.netlib.org/lapack/) - Linear Algebra Package. Она выполняет высокоуровневые операции (например матричные разложения), которые опираются на вызовы функций из BLAS.3. [Intel MKL](https://software.intel.com/en-us/intel-mkl) - Math Kernel Library. Эта реализация является оптимизированной версией BLAS и LAPACK для процессоров Intel. Доступна в Anaconda: ```conda install mkl```MATLAB использует Intel MKL по умолчанию.4. OpenBLAS – это оптимизированная версия BLAS, основанная на [GotoBLAS](https://en.wikipedia.org/wiki/GotoBLAS). 5. PyTorch [поддерживает](https://pytorch.org/docs/stable/torch.htmlblas-and-lapack-operations) некоторые вызовы из BLAS и LAPACK6. Для GPU реализован [cuBLAS](https://docs.nvidia.com/cuda/cublas/index.html).Для сравнения производительности OpenBLAS и Intel MKL, см [этот обзор](https://software.intel.com/en-us/articles/performance-comparison-of-openblas-and-intel-math-kernel-library-in-r)  Более быстрые алгоритмы для умножения матрицНапомним, что сложность алгоритма матричного умножения $\mathcal{O}(n^3)$ операций.Однако, данных требуется только $\mathcal{O}(n^2)$.**Q:** возможно ли уменьшить число операций необходимое для умножения матриц до $\mathcal{O}(n^2)$? **A**: алгоритм умножения матриц со сложностью $\mathcal{O}(n^2)$ до сих пор не найден. * Метод Штрассена – $\mathcal{O}(n^{2.807\dots})$, иногда используется на практике* [Текущий мировой рекорд](http://arxiv.org/pdf/1401.7714v1.pdf) $\mathcal{O}(n^{2.37\dots})$ – большая константа, не практичный, основан на [алгоритме Coppersmith-Winograd'a](https://en.wikipedia.org/wiki/Coppersmith%E2%80%93Winograd_algorithm).Рассмотрим метод Штрассена более детально.  Наивный алгоритмПусть $A$ и $B$ две матрицы $2\times 2$. Наивное умножение $C = AB$$$\begin{bmatrix} c_{11} & c_{12} \\ c_{21} & c_{22}  \end{bmatrix}  =\begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22}  \end{bmatrix}\begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22}  \end{bmatrix} =\begin{bmatrix} a_{11}b_{11} + a_{12}b_{21} & a_{11}b_{21} + a_{12}b_{22} \\ a_{21}b_{11} + a_{22}b_{21} & a_{21}b_{21} + a_{22}b_{22} \end{bmatrix}$$состоит из $8$ умножений и $4$ сложений.  Алгоритм ШтрассенаВ работе [Gaussian elimination is not optimal](http://link.springer.com/article/10.1007%2FBF02165411?LI=true) (1969 г.) Штрассен обнаружил, что можно вычислить матрицу $C$, используя 18 сложений и только 7 умножений:$$\begin{split}c_{11} &= f_1 + f_4 - f_5 + f_7, \\c_{12} &= f_3 + f_5, \\c_{21} &= f_2 + f_4, \\c_{22} &= f_1 - f_2 + f_3 + f_6,\end{split}$$где$$\begin{split}f_1 &= (a_{11} + a_{22}) (b_{11} + b_{22}), \\f_2 &= (a_{21} + a_{22}) b_{11}, \\f_3 &= a_{11} (b_{12} - b_{22}), \\f_4 &= a_{22} (b_{21} - b_{11}), \\f_5 &= (a_{11} + a_{12}) b_{22}, \\f_6 &= (a_{21} - a_{11}) (b_{11} + b_{12}), \\f_7 &= (a_{12} - a_{22}) (b_{21} + b_{22}).\end{split}$$К счастью, эти формулы работают даже если $a_{ij}$ и $b_{ij}$, $i,j=1,2$ – блочные матрицы.Таким образом, схема алгоритма Штрассена такова - Сначала мы  разбиваем матрицы $A$ и $B$ размера $n\times n$, $n=2^d$  на 4 блока размера $\frac{n}{2}\times \frac{n}{2}$- После чего вычисляем произведения по указанным выше формулам рекурсивноЭто приводит нас к идее **разделяй и властвуй**.  Сложность метода Штрассена Число умноженийОбозначим $M(n)$ число умножений необходимых для вычисления произведения двух матриц $n\times n$ с помощью стратегии "разделяй и властвуй".Тогда для наивного алгоритма получаем следующее число умножений$$M_\text{naive}(n) = 8 M_\text{naive}\left(\frac{n}{2} \right) = 8^2 M_\text{naive}\left(\frac{n}{4} \right)= \dots = 8^{d-1} M(1) = 8^{d} = 8^{\log_2 n} = n^{\log_2 8} = n^3$$Итак, даже используя стратегию "разделяй и властвуй" мы не улучшаем сложность в $\mathcal{O}(n^3)$ операций.Вычислим число умножений в методе Штрассена:$$M_\text{strassen}(n) = 7 M_\text{strassen}\left(\frac{n}{2} \right) = 7^2 M_\text{strassen}\left(\frac{n}{4} \right)= \dots = 7^{d-1} M(1) = 7^{d} = 7^{\log_2 n} = n^{\log_2 7}$$  Число сложенийДля алгоритма Штрассена:$$A_\text{strassen}(n) = 7 A_\text{strassen}\left( \frac{n}{2} \right) + 18 \left( \frac{n}{2} \right)^2$$поскольку на первом уровне нам нужно сложить матрицы размера $\frac{n}{2}\times \frac{n}{2}$ 18 раз, а потом идти на следующий уровень рекурсии для каждого из 7 умножений. Таким образом,$$\begin{split}A_\text{strassen}(n) =& 7 A_\text{strassen}\left( \frac{n}{2} \right) + 18 \left( \frac{n}{2} \right)^2 = 7 \left(7 A_\text{strassen}\left( \frac{n}{4} \right) + 18 \left( \frac{n}{4} \right)^2 \right) + 18 \left( \frac{n}{2} \right)^2 =7^2 A_\text{strassen}\left( \frac{n}{4} \right) + 7\cdot 18 \left( \frac{n}{4} \right)^2 +  18 \left( \frac{n}{2} \right)^2 = \\=& \dots = 18 \sum_{k=1}^d 7^{k-1} \left( \frac{n}{2^k} \right)^2 = \frac{18}{4} n^2 \sum_{k=1}^d \left(\frac{7}{4} \right)^{k-1} = \frac{18}{4} n^2 \frac{\left(\frac{7}{4} \right)^d - 1}{\frac{7}{4} - 1} = 6 n^2 \left( \left(\frac{7}{4} \right)^d - 1\right) \leqslant 6 n^2 \left(\frac{7}{4} \right)^d = 6 n^{\log_2 7}\end{split}$$(так как $4^d = n^2$ и $7^d = n^{\log_2 7}$).Асимптотику для $A(n)$ также можно найти с помощью следующей [теоремы](https://en.wikipedia.org/wiki/Master_theorem).  Итоговая сложность метода ШтрассенаИтоговая сложность –  $M_\text{strassen}(n) + A_\text{strassen}(n)=$ $7 n^{\log_2 7}$. Метод Штрассена становится быстрее наивного алгоритма, если$$\begin{split}2n^3 &> 7 n^{\log_2 7}, \\n &> 667,\end{split}$$поэтому не самая хорошая идея спускаться на самый низкий уровень рекурсии.  Является ли метод Штрассена практически важным?- Долгое время считалось, что метод Штрассена имеет только теоретический интерес и плохо применим на практике- Статья [Strassen algorithm reloaded](http://jianyuhuang.com/papers/sc16.pdf)опровергает утверждение, что метод Штрассена не практичен, и предлагает реализацию, которая работает не только на больших матрицах.  Резюме первой части лекции- Умножение матриц – одна из основных операций в вычислительной линейной алгебре. Если необходима высокая эффективность, необходимо использовать блочные алгоритмы- Это связано с оптимальным использованием иерархии памяти- Концепция блочных алгоритмов- Алгоритм Штрассена  Часть 2. Матричные разложения и методы их вычисления  Основные матричные разложения- SVD- LU разложение- QR разложение- Разложение Шура  Общая концепция матричных разложенийВычислительная линейная алгебра занимается решением следующих задач:- Решение линейных систем $Ax = f$- Вычисление собственных значений и векторов- Вычисление сингулярных значений и векторов - Вычисление обратных матриц и иногда детерминантов- Вычисление **матричных функций** таких как $\exp(A), \cos(A)$ (это не поэлементные функции!)Для решения таких задач мы представляемм матрицу в виде суммы и/или произведения матриц **более простой структуры**, таких что мы можем решить эти задачи быстрее и/или более устойчивым образом.Что такое **более простая структура**?  Матрицы простой структурыМы уже упоминали некоторые классы структурированных матриц. Для плотных матриц это следующие матрицы- **унитарные матрицы**- **нижне- верхнетреугольные матрицы** - **диагональные матрицы**  PLU разложениеЛюбая невырожденная матрица моджет быть разложена как$$A = P L U,$$где $P$ – матрица перестановки, $L$ – **нижнетреугольная матрица**, $U$ – **верхнетреугольная матрица**.**Основная цель** вычисления LU разложения – это решение систем линейных уравнений, потому что$$    A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, $$что сводится к решению двух линейных систем $$     L y = f,  \quad U x = y$$с нижне- и верхнетреугольной матрицей, соответственно.   Положительно определённые матрицы и разложение ХолецкогоЕсли матрица эрмитова и положительно определена, то есть$$A = A^*, \quad (Ax, x) > 0, \quad x \ne 0,$$тогда она может быть представлена в виде$$A = RR^*,$$где $R$ – нижнетреугольная матрица.Это нам понадобится для QR разложения.  QR разложениеКак следует из названия, это представление матрицы в виде произведения$$    A = Q R, $$где $Q$ – матрица с **ортогональными столбцами** и $R$ – **верхнетреугольная**.  Размеры матриц: $Q$ – $n \times m$, $R$ – $m \times m$, если $n\geq m$.QR разложение определено для любой **прямоугольной матрицы**.  Приложения QR разложенияЭто разложение играет ключевую роль при решении многих задач, например:- Получение ортогонального базиса в линейном пространстве- Используется для препроцессинга при вычислении SVD- QR алгоритм для вычисления собственных векторов и собственных значений ([один из 10 самых важных алгоритмов ХХ века](https://archive.siam.org/pdf/news/637.pdf)) основан на вычислении QR разложения- Решение переопределённых систем линейных уравнений  QR разложение и линейная задача наименьших квадратов: напоминаниеПусть нам надо решить переопределённую систему $Ax = f$ в смысле наименьших квадратов невязок  $$\Vert A x - f \Vert_2 \rightarrow \min_x,$$где $A$ – $n \times m$ и $n \geq m$.Тогда мы разложим её $$A = Q R,$$и используем выражение для псевдообратной матрицы в случае матрицы полного ранга:$$x = A^{\dagger}b = (A^*A)^{-1}A^*b = ((QR)^*(QR))^{-1}(QR)^*b = (R^*Q^*QR)^{-1}R^*Q^*b = R^{-1}Q^*b. $$тогда $x$ может быть найдён из решения следующей системы$$R x = Q^*b$$Заметим, что это квадратная система линейных уравнений с верхнетреугольной матрицей. Какова сложность решения такой системы?  Существование QR разложения**Теорема.**Каждая матрица $n \times m$ может быть представлена в виде QR разложения. Существует несколько способов доказательства и вычисления:- Теоретический: используя матрицу Грама и разложение Холецкого- Геометрический: используя ортогонализацию Грама-Шмидта- Практический: используя преобразования Гивенса/Хаусхолдера  Доказательство с использованием разложения ХолецкогоЕсли у нас есть разложение$$A = QR,$$тогда $A^* A = ( Q R)^* (QR)  = R^* (Q^* Q) R = R^* R$, матрица $A^* A$ называется **матрицей Грама**, и её элементы – скалярные произведения столбцов матрицы $A$.    Случай матрицы полного рангаПусть $A$ имеет **полный столбцовый ранг**. Тогда легко показать, что $A^* A$ положительно определена:$$   (A^* A y, y) = (Ay, Ay) = \Vert Ay \Vert^2  > 0, \quad y\not = 0.$$Поэтому, $A^* A = R^* R$ всегда существует.Тогда матрица $A R^{-1}$ унитарна:  $$   (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I.$$  Случай матрицы неполного ранга- QR разложение по-прежнему существует.- Для любой матрицы неполного ранга существует последовательность матриц полного ранга $A_k$ такая что $A_k \rightarrow A$ (почему?).- Каждая $A_k$ может быть разложена $A_k = Q_k R_k$. - Множество унитарных матриц образует компакт, таким образом найдётся сходящаяся подпоследовательность $Q_{n_k} \rightarrow Q$ (почему?), и $Q^* A_k \rightarrow Q^* A = R$, которая верхнетреугольная.  Устойчивость вычисления QR разложения с помощью разложения Холецкого Итак, простейший способ вычисления QR разложения следующий$$A^* A = R^* R,$$и$$Q = A R^{-1}.$$Это **плохая идея** с точки зрения устойчивости. Покажем, почему это так.<jupyter_code>import numpy as np
n = 50
r = 4
a = [[1.0 / (i + j + 0.5) for i in range(r)] for j in range(n)]
a = np.array(a)
q, Rmat = np.linalg.qr(a)
e = np.eye(r)
print('Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e))
gram_matrix = a.T.dot(a)
Rmat1 = np.linalg.cholesky(gram_matrix)
q1 = np.dot(a, np.linalg.inv(Rmat1.T))
print('Via Cholesky:', np.linalg.norm(np.dot(q1.T, q1) - e))<jupyter_output>Built-in QR orth 5.061169633223855e-16
Via Cholesky: 3.2416882870748012e-12<jupyter_text>Второй способ: ортогонализация Грама-Шмидта- QR разложение – это способ записи процесса ортогонализации Грама-Шмидта- Для данного набора векторов $a_1, \ldots, a_m$ мы хотим найти ортонормированный базис $q_1, \ldots, q_m$ такой чтобы каждый вектор $a_i$ представлялся как линейная комбинация векторов из базиса.  **Метод Грама-Шмидта:**1. $q_1 := a_1/\Vert a_1 \Vert$2. $q_2 := a_2 - (a_2, q_1) q_1, \quad q_2 := q_2/\Vert q_2 \Vert$3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2, \quad q_3 := q_3/\Vert q_3 \Vert$4. И так далее  Заметим, что преобразование из $A$ в $Q$ имеет треугольную структуру, поскольку мы вычитаем из $k$-го вектора только предыдущие. Это следует из того, что произведение треугольных матриц – это треугольная матрица.  Модифицированный метод Грама-Шмидта- Метод Грама-Шмидта может быть очень неустойчивым (то есть генерировать векторы, которая не являются ортогональными, особенно если $q_k$ малой нормы).  Это называется **потеря ортогональности**.  - Этот недостаток метода Грама-Шмидта исправляется с помощью **модифицированного метода Грама-Шмидта**. Вместо вычисления$$q_k := a_k - (a_k, q_1) q_1 - \ldots - (a_k, q_{k-1}) q_{k-1}$$мы будем вычислять это выражение шаг за шагом. Сначала присвоим $q_k := a_k$, затем последовательно ортогонализуем:$$   q_k := q_k - (q_k, q_1)q_1, \quad q_k := q_{k} - (q_k,q_2)q_2, \ldots$$- В точной арифметике, это одинаковые алгоритмы. В неточной арифметике они абсолютно разные!- Заметим, что сложность по-прежнему $\mathcal{O}(nm^2)$ операций  Анонс следующей лекции- Собственные векторы и собственные значения- Какие задачи сводятся к их вычислению- Степенной метод- Теорема Гершгорина- Разложение Шура и алгоритм его вычисления<jupyter_code>from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13554003888987096, 'fraction_numerical': 0.030007070885628425, 'mean_word_length': 4.766819571865444, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 28, 'https://': 13, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '186020', 'license': 'MIT', 'n_tokens_mistral': 9951, 'n_tokens_neox': 11048, 'n_words': 2544, 'path': 'lectures/lecture4/lecture4.ipynb', 'repo_name': 'AlexanderSLevin/nla2020_ozon', 'size': '22633'}