\documentclass[11pt,english,russian]{article} \usepackage[english,russian]{babel} \usepackage{fontspec} \usepackage[a4paper, bindingoffset=0.0in, left=25mm, right=25mm, top=25mm, bottom=25mm, footskip=10mm]{geometry} \usepackage{enumitem} \usepackage{epsfig} \usepackage{epsf} \usepackage{amsmath,amssymb,amsfonts,amsthm,bm} \usepackage{mathtools} \usepackage{mathrsfs} \usepackage{array} \usepackage{tabularx} \usepackage{multirow} \usepackage{longtable} \usepackage{pstricks} \usepackage{lipsum} \usepackage{ulem} \usepackage{blindtext} \usepackage{setspace} \usepackage{graphicx} \usepackage{algorithm} \usepackage{algpseudocode} \usepackage{booktabs} \usepackage{geometry} \usepackage{float} \setromanfont{Times New Roman} \setsansfont{Arial} % --- Геометрия и оформление --- % \usepackage{geometry} % \geometry{left=2.5cm, right=2.5cm, top=2.5cm, bottom=2.5cm} \usepackage{titlesec} \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} \fancyhead[L]{\leftmark} \fancyfoot[C]{\thepage} % --- Гиперссылки --- \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=tumblue, urlcolor=tumblue, citecolor=tumblue } % --- Цвета и графика --- \usepackage{xcolor} \definecolor{tumblue}{RGB}{0, 101, 165} \definecolor{alertred}{RGB}{180, 30, 30} \definecolor{successgreen}{RGB}{30, 130, 30} \definecolor{warmorange}{RGB}{200, 120, 0} \usepackage{tikz} \usetikzlibrary{arrows.meta, angles, quotes, calc} % --- Пользовательские команды --- \newcommand{\N}{\mathbb{N}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \renewcommand{\vec}[1]{\mathbf{#1}} \newcommand{\E}{\mathbb{E}} \renewcommand{\i}{\mathrm{i}} \newcommand{\cond}{\operatorname{cond}} \newcommand{\rank}{\operatorname{rank}} \newcommand{\tr}{\operatorname{tr}} \newcommand{\diag}{\operatorname{diag}} \newcommand{\sign}{\operatorname{sign}} \newcommand{\argmin}{\operatorname{arg\,min}} \newcommand{\diff}{\mathrm{d}} \newcommand{\re}{\operatorname{Re}} \newcommand{\im}{\operatorname{Im}} \newcommand{\argg}{\operatorname{arg}} \titleformat{\section}{\Large\bfseries\color{tumblue}}{\thesection}{1em}{} \titleformat{\subsection}{\large\bfseries\color{tumblue!80!black}}{\thesubsection}{1em}{} \titleformat{\subsubsection}{\normalsize\bfseries\color{tumblue!60!black}}{\thesubsubsection}{1em}{} % Окружение для важных замечаний \usepackage{tcolorbox} \tcbuselibrary{skins, breakable} \newtcolorbox{warningbox}[1][]{ colback=alertred!5!white, colframe=alertred, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{tipbox}[1][]{ colback=tumblue!5!white, colframe=tumblue, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{successbox}[1][]{ colback=successgreen!5!white, colframe=successgreen, fonttitle=\bfseries, title=#1, breakable } \newtcolorbox{orangebox}[1][]{ colback=warmorange!5!white, colframe=warmorange, fonttitle=\bfseries, title=#1, breakable } \begin{document} \title{\vspace{-2cm}\textbf{\color{tumblue} Математика, Вычислительные методы и \\ Обработка сигналов}\\[0.5cm] \large Скрипт для студентов химического факультета TUM} \author{Ильгиз Ибрагимов \& Елена Ибрагимова \& ИИ\footnote{ассистент по правкам и переводу}} \date{Сентябрь 2026 \\ DOI: 10.5281/zenodo.23002397} \maketitle \thispagestyle{empty} \vfill \begin{center} \textit{Посвящается нашим дочерям Марие и Алевтине. \\ Пусть математика станет для вас не набором формул, \\ а языком, на котором написана Вселенная.} \end{center} \newpage \tableofcontents \newpage \section{Предисловие} Дорогие наши дочери! Нам так хотелось, чтобы в вашем жизненном пути всё было понятно и просто. Мы знаем, что часто математику объясняют, специально запутывая человека терминами и сложными формулами --- так, что за деревьями не видно леса. Этот скрипт --- наша попытка дать обзор математических и вычислительных идей, особенно полезных современному химику, который работает с экспериментальными данными, моделированием и вычислениями. Наша цель --- чтобы у тебя была \textit{понятная картина} современного уровня математических знаний. Не всех знаний --- это невозможно в одной книге --- а именно тех, которые применяются в химии, биохимии, спектроскопии, молекулярном моделировании и обработке данных. Мы прошли вместе длинный путь: \begin{itemize} \item От самых основ --- множеств чисел, комплексных чисел, векторов и матриц, \item Через линейную алгебру --- SVD, обусловленность, методы решения систем, \item К нелинейным задачам --- оптимизация, градиенты, метод Ньютона, \item Через обработку сигналов --- Фурье, Прони, сжатые измерения, \item К численным методам --- интегрирование, конечные элементы, базисные функции, \item И наконец --- к машинному обучению, нейросетям и современной вычислительной технике. \end{itemize} Всё это --- не просто набор разрозненных тем. Это --- \textbf{единый язык}, на котором говорит современная наука. И если вы овладеете этим языком --- перед вами откроются двери, о которых вы сейчас даже не подозреваете. \subsection{Как читать эту книгу} Хотя текст написан так, чтобы сделать его максимально просто читаемым, без детальных математических доказательств, некоторые места могут быть слишком заумно сформулированными. Это нормально --- математика не даётся сразу. Поэтому к тексту прилагается исходник этой книги в \LaTeX. Если что-то будет не сильно понятно: \begin{enumerate} \item Найдите соответствующий текст в исходнике, \item Скопируйте его в какой-нибудь чат (Qwen, DeepSeek, Kimi, ChatGPT, GROK, Gemini), \item Попросите его понятнее рассказать про это, или задайте вопрос, который вы тут не понимаете. \end{enumerate} Также с помощью этих чатов можно перевести весь текст на английский или немецкий --- для этого надо просто попросить перевести с сохранением \LaTeX-разметки и полученный текст скомпилировать командой \texttt{xelatex NumMathForChemists.tex} в PDF-файл. \begin{tipbox}[Совет] Не пытайтесь прочитать всё за один присест. Математика --- это как музыка: её надо ``играть'', то есть решать задачи, пробовать код, экспериментировать. Если какая-то глава кажется сложной --- вернитесь к ней через неделю, и вы удивитесь, насколько проще она станет. \end{tipbox} % \subsection{Зачем химику математика?} % % Прежде чем мы начнем разбирать числа, матрицы, операторы и преобразования Фурье, полезно ответить на самый естественный вопрос: % \textit{Зачем всё это вообще нужно химику?} % Ответ очень простой: потому что современная химия всё чаще работает не непосредственно с веществами, а с \textbf{данными, моделями и вычислениями}. % % Представь совершенно обычный эксперимент. Мы помещаем образец в прибор и получаем некоторый сигнал. Например, спектрометр может вернуть несколько тысяч чисел. Что мы хотим сделать с этими числами? % % Мы хотим удалить шум, найти пики, определить их положение и интенсивность, сравнить полученный спектр с известными спектрами, определить концентрацию вещества или восстановить параметры молекулы. То есть эксперимент можно представить в очень упрощенном виде: % % $$ % \boxed{ % \text{вещество} % \longrightarrow % \text{измерение} % \longrightarrow % \text{данные} % \longrightarrow % \text{математическая модель} % \longrightarrow % \text{химический вывод} % } % $$ % % И почти каждый переход в этой цепочке требует математики. Поэтому весь дальнейший текст можно воспринимать как путешествие по одной цепочке: % % $$ % \boxed{ % \begin{array}{c} % \text{числа}\\ % \downarrow\\ % \text{векторы и матрицы}\\ % \downarrow\\ % \text{операторы}\\ % \downarrow\\ % \text{линейные и нелинейные задачи}\\ % \downarrow\\ % \text{статистика и неопределенность}\\ % \downarrow\\ % \text{сигналы и Фурье}\\ % \downarrow\\ % \text{численные методы}\\ % \downarrow\\ % \text{сжатие и извлечение информации}\\ % \downarrow\\ % \text{машинное обучение}\\ % \downarrow\\ % \text{современная вычислительная химия} % \end{array} % } % $$ % Не обязательно запоминать все формулы. Гораздо важнее постепенно научиться узнавать математическую структуру задачи. % % \begin{itemize} % \item Если появляется набор чисел, спроси: \textit{Что это за математический объект?} % \item Если появляется множество измерений: \textit{Можно ли представить их как вектор или матрицу?} % \item Если есть неизвестные параметры: \textit{Как построить модель и найти эти параметры?} % \item Если есть шум: \textit{Как отделить информацию от случайной ошибки?} % \item Если данных слишком много: \textit{Можно ли найти скрытую структуру или уменьшить размерность?} % \end{itemize} % % % Именно такой способ мышления, а не запоминание большого количества формул, является главной целью этого скрипта. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Введение: Обозначения и природа чисел} \subsection{Иерархия числовых множеств} Прежде чем мы начнем говорить о сложных вещах, давай договоримся о языке. В математике мы группируем числа в множества, обладающие определенными свойствами. Мы будем использовать следующие стандартные обозначения: \begin{itemize} \item $\N = \{1, 2, 3, \dots\}$ --- \textbf{натуральные числа}. Используются для счета. \item $\Z = \{\dots, -2, -1, 0, 1, 2, \dots\}$ --- \textbf{целые числа}. Добавляют ноль и отрицательные числа, позволяя описывать долг и имущество, температуру ниже нуля и т.д. \item $\Q$ --- \textbf{рациональные числа}. Любое число, которое можно представить в виде дроби $\frac{p}{q}$, где $p \in \Z$, $q \in \Z \setminus \{0\}$. \item $\R$ --- \textbf{действительные (вещественные) числа}. Включают в себя все рациональные, а также иррациональные числа (например, $\sqrt{2}, \pi, e$), которые нельзя представить в виде обычной дроби. Они непрерывно заполняют всю числовую прямую. \item $\C$ --- \textbf{комплексные числа}. О них мы подробно поговорим ниже. \end{itemize} \textit{Вложенность множеств:} $\N \subset \Z \subset \Q \subset \R \subset \C$. \subsection{Числа в компьютере} Мы привыкли, что в математике числа бесконечно точны. Но когда мы переходим к \textbf{вычислительной математике} и программированию, компьютер вынужден хранить числа в ячейках памяти фиксированного размера (в битах и байтах). \begin{itemize} \item \textbf{Целые типы} (int8, int16, int32, int64) хранят точные значения из $\Z$, но ограничены диапазоном (например, int32 может хранить числа от $-2^{31}$ до $2^{31}-1$). \item \textbf{Вещественные типы} (float32, float64 / double) хранят числа из $\R$ в экспоненциальном формате (мантисса и порядок). Из-за этого возникают ошибки округления: компьютерное $0.1 + 0.2$ не всегда в точности равно $0.3$. Понимание этого критически важно для численных методов! \end{itemize} \section{Когда одного числа мало: Комплексные числа} Наш мир настолько сложен, что не всегда все можно описать одним числом. В физике и химии нам часто приходится работать с сущностями, которые требуют \textit{двух} чисел для своего описания (например, амплитуда и фаза волны, или действительная и мнимая часть волновой функции). \subsection{Откуда они берутся?} Самый простой и исторический способ познакомиться с ``числами-парами'' --- это попытаться решить уравнение: \[ x^2 + 1 = 0 \quad \Rightarrow \quad x^2 = -1 \] В множестве действительных чисел $\R$ квадрат любого числа неотрицателен. Корня из $-1$ не существует. Но математики (и физики вслед за ними) сказали: ``А что, если мы просто \textit{придумаем} такое число, которое при умножении само на себя дает $-1$?''. Так появилась \textbf{мнимая единица} $\mathrm{i}$. Мы не ``извлекаем корень'' из $-1$, мы \textit{определяем} новую сущность: \[ \mathrm{i}^2 = -1 \] \textit{Важное замечание: правило $\sqrt{a}\sqrt{b} = \sqrt{ab}$ работает только для $a,b \ge 0$. Поэтому мы не пишем $\sqrt{-1}\sqrt{-1} = \sqrt{(-1)(-1)} = 1$. Мы просто принимаем $\mathrm{i}^2 = -1$ как факт.} \subsection{Алгебраическая форма и комплексная плоскость} Любое комплексное число $z \in \C$ можно записать в виде: \[ z = x + \mathrm{i}y \] где $x = \re(z)$ --- \textbf{действительная часть}, а $y = \im(z)$ --- \textbf{мнимая часть}. Обе части $x, y \in \R$. Геометрически комплексное число --- это точка на \textbf{комплексной плоскости}. \begin{itemize} \item По горизонтальной оси (ось $\re$) откладывается действительная часть $x$. \item По вертикальной оси (ось $\im$) откладывается мнимая часть $y$. \end{itemize} Комплексное число также можно рассматривать как \textbf{вектор}, идущий из начала координат $(0,0)$ в точку $(x,y)$. \subsection{Арифметика комплексных чисел} \textbf{Сложение} интуитивно понятно. Если у нас есть два числа $z_1 = x_1 + \mathrm{i}y_1$ и $z_2 = x_2 + \mathrm{i}y_2$, то мы просто складываем их действительные и мнимые части отдельно: \[ z_1 + z_2 = (x_1 + x_2) + \mathrm{i}(y_1 + y_2) \] Геометрически это работает как \textbf{правило треугольника} для векторов: мы строим два вектора от начала координат, переносим начало второго вектора в конец первого, и сумма --- это вектор от начала первого до конца второго. \textbf{Умножение} в алгебраической форме выглядит чуть сложнее. Мы раскрываем скобки как в обычной алгебре, помня, что $\mathrm{i}^2 = -1$: \begin{align*} z_1 \cdot z_2 &= (x_1 + \mathrm{i}y_1)(x_2 + \mathrm{i}y_2) \\ &= x_1 x_2 + \mathrm{i}x_1 y_2 + \mathrm{i}y_1 x_2 + \mathrm{i}^2 y_1 y_2 \\ &= (x_1 x_2 - y_1 y_2) + \mathrm{i}(x_1 y_2 + y_1 x_2) \end{align*} Но чтобы понять \textit{истинную магию} умножения комплексных чисел, нам нужно перейти к другой форме записи. \subsection{Тригонометрическая и экспоненциальная формы} Вместо координат $(x,y)$ точку на плоскости можно задать с помощью полярных координат: расстояния от начала координат (модуля) $r$ и угла $\varphi$ относительно положительного направления оси $\re$. Связь с алгебраической формой очевидна из тригонометрии: \[ x = r \cos \varphi, \quad y = r \sin \varphi \] Тогда комплексное число принимает \textbf{тригонометрическую форму}: \[ z = r(\cos \varphi + \mathrm{i} \sin \varphi) \] \subsubsection{Формула Эйлера и экспоненциальная форма} Здесь на сцену выходит величайшая формула математики --- \textbf{формула Эйлера}: \[ e^{\mathrm{i}\varphi} = \cos \varphi + \mathrm{i} \sin \varphi \] Благодаря ей, любое комплексное число можно записать в невероятно удобной \textbf{экспоненциальной форме}: \[ z = r e^{\mathrm{i}\varphi} \] \subsubsection{Геометрический смысл умножения} Давай посмотрим, что происходит при умножении двух чисел в экспоненциальной форме: \[ z_1 = r_1 e^{\mathrm{i}\varphi_1}, \quad z_2 = r_2 e^{\mathrm{i}\varphi_2} \] \[ z_1 \cdot z_2 = (r_1 e^{\mathrm{i}\varphi_1}) \cdot (r_2 e^{\mathrm{i}\varphi_2}) = (r_1 r_2) e^{\mathrm{i}(\varphi_1 + \varphi_2)} \] \textbf{Вывод, который нужно запомнить навсегда:} Умножение комплексных чисел --- это \textbf{умножение их длин (модулей)} и \textbf{сложение их углов (аргументов)}! Геометрически: умножение на комплексное число --- это поворот вектора на угол $\varphi$ и его растяжение/сжатие в $r$ раз. \subsection{Комплексные экспоненты и логарифмы} \subsubsection{Свойства экспоненты} Поскольку $e^{\mathrm{i}\varphi}$ ведет себя как обычная экспонента, для нее работают все стандартные правила: \[ e^{z_1} \cdot e^{z_2} = e^{z_1 + z_2}, \quad \frac{e^{z_1}}{e^{z_2}} = e^{z_1 - z_2}, \quad (e^z)^n = e^{nz} \] Для комплексного числа $z = x + \mathrm{i}y$ экспонента раскладывается на действительную и мнимую части: \[ e^z = e^{x + \mathrm{i}y} = e^x \cdot e^{\mathrm{i}y} = e^x (\cos y + \mathrm{i} \sin y) \] Модуль этого числа равен $e^x$, а аргумент --- $y$. \subsubsection{Комплексный логарифм} Логарифм --- это операция, обратная экспоненте. Если $e^w = z$, то $w = \ln z$. Пусть $z = r e^{\mathrm{i}\varphi}$ и $w = u + \mathrm{i}v$. Тогда: \[ e^{u + \mathrm{i}v} = r e^{\mathrm{i}\varphi} \quad \Rightarrow \quad e^u e^{\mathrm{i}v} = r e^{\mathrm{i}\varphi} \] Приравнивая модули и аргументы, получаем: \[ e^u = r \quad \Rightarrow \quad u = \ln r \] \[ v = \varphi + 2\pi k, \quad k \in \Z \] Таким образом, \textbf{комплексный логарифм} многозначен: \[ \ln z = \ln|z| + \mathrm{i}(\argg z + 2\pi k) \] Главное значение логарифма (при $k=0$ и $\argg z \in (-\pi, \pi]$) обозначается как $\text{Ln } z$. Многозначность возникает из-за того, что поворот на $2\pi$ возвращает нас в ту же точку на комплексной плоскости. \subsection{Шпаргалка: Тригонометрия и Гиперболические функции} Гиперболические функции часто пугают студентов-химиков, но на самом деле они --- близкие родственники обычных синусов и косинусов, просто ``живущие'' в комплексной плоскости. \subsubsection{Определения через экспоненту} \begin{align*} \cos z &= \frac{e^{\mathrm{i}z} + e^{-\mathrm{i}z}}{2}, & \sin z &= \frac{e^{\mathrm{i}z} - e^{-\mathrm{i}z}}{2\mathrm{i}} \\ \cosh z &= \frac{e^{z} + e^{-z}}{2}, & \sinh z &= \frac{e^{z} - e^{-z}}{2} \end{align*} \subsubsection{Связь тригонометрических и гиперболических функций} Если мы подставим $\mathrm{i}z$ вместо $z$ в определения, то увидим удивительную симметрию (формулы Осборна): \begin{align*} \cos(\mathrm{i}z) &= \cosh z, & \cosh(\mathrm{i}z) &= \cos z \\ \sin(\mathrm{i}z) &= \mathrm{i}\sinh z, & \sinh(\mathrm{i}z) &= \mathrm{i}\sin z \end{align*} \textit{Запоминалка:} При переходе от тригонометрии к гиперболике аргумент умножается на $\mathrm{i}$, а перед синусом/косинусом могут появляться мнимые единицы. Гиперболические функции описывают не колебания (как синус), а экспоненциальный рост/спад (как $e^x$), что критически важно для затухающих процессов. \subsection{Зачем это нужно в реальной химии и физике?} Может показаться, что мнимые числа --- это просто красивая математическая абстракция. Но наш мир устроен так, что \textit{без} мнимого пространства мы не смогли бы описать реальные вещи. \subsubsection{Квантовая механика и орбитали} Уравнение Шрёдингера, которое описывает поведение электронов в атомах, содержит мнимую единицу $\mathrm{i}$ в явном виде: \[ \mathrm{i}\hbar \frac{\partial \Psi}{\partial t} = \hat{H}\Psi \] Волновая функция $\Psi$ --- комплекснозначная. Для основного состояния атома водорода (1s-орбиталь) волновая функция действительная, и ее можно нарисовать в реальном пространстве. Но как только электрон переходит в возбужденное состояние (например, 2p-орбиталь с магнитным квантовым числом $m \neq 0$), его волновая функция содержит фазовый множитель $e^{\mathrm{i}m\varphi}$. Без комплексных чисел описать угловой момент электрона и тонкую структуру спектров просто невозможно. \subsubsection{Ядерный Магнитный Резонанс (ЯМР) и Фурье-анализ} Вы будете много работать с ЯМР-спектроскопией. Когда вы помещаете образец в магнитное поле и даете радиочастотный импульс, ядра начинают прецессировать. Детектор регистрирует сигнал, затухающий со временем --- это называется \textbf{FID} (Free Induction Decay). Сигнал от одного типа ядер выглядит как затухающая осцилляция: \[ S(t) = A e^{-t/T_2} \cos(\omega t) \] Если в образце много разных ядер, сигнал превращается в кашу из множества наложенных друг на друга косинусов. Как понять, какие частоты $\omega$ там спрятаны? Здесь на помощь приходит \textbf{Преобразование Фурье}. Но работать с косинусами сложно. Гораздо проще перейти к комплексным экспонентам, используя формулу Эйлера: \[ \cos(\omega t) = \frac{e^{\mathrm{i}\omega t} + e^{-\mathrm{i}\omega t}}{2} \] В ЯМР-спектрометрах используют квадратурное детектирование, которое позволяет измерять сразу комплексный сигнал: \[ S_{complex}(t) = A e^{-t/T_2} e^{\mathrm{i}\omega t} \] Когда мы применяем к этому сигналу быстрое преобразование Фурье (БПФ), время $t$ переходит в частоту $\omega$. Длинная осциллирующая функция во временной области превращается в \textbf{узкий пик (Лоренцеву кривую)} в частотной области! Каждому типу ядер соответствует своя частота $\omega$, и на итоговом спектре ЯМР мы видим отдельные пики. Вся современная обработка сигналов (от МРТ в медицине до аудио MP3) работает именно благодаря магии комплексных чисел и экспонент Эйлера. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Когда одного измерения мало: Векторы и Нормы} \subsection{От плоскости к N-мерному пространству} Что же, если у нас есть комплексные числа (по сути, пары чисел), то, наверное, есть и что-то более сложное? Кватернионы? Октонионы? Да, действительно, математики обобщили числа с плоскости на большие размерности. Но за пределами комплексных чисел и кватернионов (которые, кстати, отлично прижились в робототехнике и 3D-графике для описания вращений в нашем трехмерном пространстве) специальные ``типы чисел'' практически не используются. Вместо этого люди просто группируют числа в наборы и называют их \textbf{N-мерными векторами}. Это просто N-мерное пространство. $N$ может быть равно двум (тогда это плоскость), трем (наше обычное физическое пространство) или очень-очень большим. Такие векторы мы будем обозначать жирным шрифтом, например $\vec{v} \in \R^N$. Это вектор: \[ \vec{v} = (v_1, v_2, \dots, v_N)^T \] \textit{Здесь $T$ означает транспонирование, чтобы записать столбец в одну строку текста}. Мы сами вольны придумывать, что делать с этим пространством, но первое, что нам нужно --- это понять, как измерять ``размер'' или ``длину'' такого вектора. \subsection{Зачем нам векторы? Примеры из химии} На первый взгляд, вектор --- это просто набор чисел. Но числа-то откуда-то пришли! Допустим, это данные с хроматографа или ЯМР-спектрометра. Представь, что мы хотим узнать, какой же тут самый ``яркий'' пик. Так это же просто максимум среди этих чисел, верно? А теперь давай возьмем хроматограмму какой-нибудь сложной грязи и хроматограмму одного чистого вещества. Мы вколем в хроматограф и то, и другое в одинаковом молярном объеме (пусть условимся, что детектор одинаково чувствует все молекулы). \begin{itemize} \item В случае чистого вещества мы получим один красивый, очень острый и высокий пик. \item В случае грязи мы получим гору пиков, может быть, даже слившихся в один сплошной пологий горб. \end{itemize} Но вот что интересно: \textbf{площадь} (интегральная интенсивность) обеих хроматограмм будет одинакова! Сам пик чистого вещества будет очень высоким по сравнению с горой, но сумма всех откликов равна количеству вещества. Из этой задачи естественно возникают три разных способа оценить наш вектор данных: \begin{enumerate} \item \textbf{Максимум.} Нам важна максимальная концентрация (или максимальное поглощение, если пик ``смотрит'' вниз). Мы берем самое большое значение. \item \textbf{Сумма.} Нам важно общее количество вещества. Мы складываем все значения (обычно по модулю, так как базовая линия может быть не идеальной). \item \textbf{Корень из суммы квадратов.} А это зачем? Представь, что мы описываем молекулу не хроматограммой, а тремя параметрами: \textit{полярность, молярная масса, температура кипения}. Это вектор в $\R^3$. Чтобы понять, насколько две молекулы \textit{похожи} друг на друга (например, для дизайна лекарств), мы не можем просто сложить разности параметров или взять максимальную разность. Нам нужно именно \textbf{евклидово расстояние} --- кратчайший путь в этом ``химическом пространстве''. Или другой пример: в физике квадратичная норма вектора (спектра) пропорциональна \textbf{полной энергии} этого сигнала. \end{enumerate} \subsection{Математический язык норм} Все эти три интуитивные понятия в математике называются \textbf{нормами векторов} и обозначаются двойными вертикальными чертами $\| \vec{v} \|$. \begin{itemize} \item \textbf{Максимальная норма ($L_\infty$):} \[ \| \vec{v} \|_\infty = \max_{1 \le k \le N} |v_k| \] \item \textbf{Манхэттенское расстояние / Сумма модулей ($L_1$):} \[ \| \vec{v} \|_1 = \sum_{k=1}^N |v_k| \] \item \textbf{Евклидова норма / Длина вектора ($L_2$):} \[ \| \vec{v} \|_2 = \sqrt{\sum_{k=1}^N |v_k|^2} \] \end{itemize} Но математики любят порядок, поэтому они объединили всё это в одну общую формулу для \textbf{$L_p$-нормы}: \[ \| \vec{v} \|_p = \left( \sum_{k=1}^N |v_k|^p \right)^{1/p} \] В этой формуле $L_1$ --- это случай, когда $p=1$. $L_2$ --- когда $p=2$. А что же будет в качестве $p$ для максимума? Да всё верно, $p \to \infty$. Если взять предел этой формулы при $p \to \infty$, то математически строго получится именно максимальный элемент вектора. \subsection{Непрерывный случай: Функции как векторы} Более того, наш вектор может быть совсем непрерывным. Тогда это на самом деле --- \textbf{функция} $f(x)$. Сейчас мы только запомним, что функцию можно рассматривать как ``бесконечномерный вектор'', значения которого заданы в каждой точке $x$. И для функций нормы работают точно так же, только сумму мы заменяем на интеграл: \[ \| f \|_p = \left( \int |f(x)|^p \, dx \right)^{1/p} \] Мы будем иногда оперировать то массивом (вектором), то функцией, а позже, в курсе обработки сигналов, увидим, что это две стороны одной медали. \subsection{Неравенство Минковского} Есть один критически важный момент, про который нужно помнить. Мы часто будем сравнивать $\| \vec{x} + \vec{y} \|$ с $\| \vec{x} \|$ и $\| \vec{y} \|$. Интуитивно понятно, что если мы пройдем из точки A в точку B, а потом из B в C, то общий путь будет не меньше, чем если бы мы пошли напрямую из A в C. Это геометрическое свойство (длина стороны треугольника не больше суммы двух других) для любых $L_p$-норм формализуется в \textbf{неравенстве Минковского}: \[ \| \vec{x} + \vec{y} \|_p \le \| \vec{x} \|_p + \| \vec{y} \|_p \] \textit{Как это использовать?} Это позволяет нам оценивать ошибки. Если $\vec{x}$ --- это истинный сигнал, а $\vec{y}$ --- это шум (ошибка измерения), то неравенство Минковского гарантирует, что норма (размер) нашего измеренного сигнала $(\vec{x}+\vec{y})$ не превысит сумму нормы истинного сигнала и нормы шума. \textit{Строгое и красивое доказательство этого неравенства можно посмотреть в \href{https://ru.wikipedia.org/wiki/Неравенство_Минковского}{Википедии}, чтобы не перегружать этот скрипт.} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Два измерения и больше: Матрицы и Тензоры} \subsection{От векторов к матрицам} Раз у нас есть вектор $\vec{a} \in \R^N$ (одномерная сущность, дискретный аналог функции $f(x)$), то логично спросить: а что если у нас функция двух переменных $f(x,y)$? Каков её дискретный аналог? А для $f(x,y,z)$? Да, двухмерные объекты в дискретном пространстве --- это \textbf{матрицы}. И как и в случае с комплексными числами, всё тут очень-очень хорошо изучено. А всё, что имеет три и более измерений, обычно называется \textbf{тензорами}. Кстати, если ты когда-нибудь столкнешься со старой научной литературой по хемометрике или анализу данных, ты можешь удивиться зоопарку названий для тензоров. Их там называют \textit{multilinear, multidimensional, procrustes, multimodal, three-way, multi-way arrays}. Не пугайся, под всеми этими красивыми словами скрывается просто многомерный массив чисел. Пока давай остановимся на двумерных объектах: функция $f(x,y)$ и матрица $\mathbf{A} \in \R^{N \times M}$. \subsection{Комплексные векторы и матрицы} Почему я до сих пор писал только $\R$? И вектор, и матрица, могут быть, конечно же, комплексными! Мы можем рассмотреть пространства $\C^N$ и $\C^{N \times M}$. Тут всё просто: вместо одного действительного числа в каждой ячейке вектора или матрицы стоит комплексная пара. Но как считать норму для комплексного числа? Да всё то же самое! Модуль комплексного числа $z = x + \i y$ вычисляется через сопряженное число $\bar{z} = x - \i y$: \[ |z| = \sqrt{x^2 + y^2} = \sqrt{z \cdot \bar{z}} \] Соответственно, $L_2$-норма комплексного вектора $\vec{z} \in \C^N$ будет выглядеть так: \[ \| \vec{z} \|_2 = \sqrt{\sum_{k=1}^N |z_k|^2} = \sqrt{\sum_{k=1}^N z_k \bar{z}_k} \] В матричной нотации это записывается через эрмитово сопряжение (транспонирование + комплексное сопряжение) $\vec{z}^H$: \[ \| \vec{z} \|_2 = \sqrt{\vec{z}^H \vec{z}} \] \subsection{Нормы матриц} А для матриц мы тоже можем задать нормы! Но тут есть важный нюанс. Большинство простых норм для матриц считаются так: мы мысленно ``разрезаем'' матрицу на строки, выстраиваем все её $N \times M$ ячеек в один длинный вектор-столбец и берем норму уже этого вектора. Самая популярная из таких норм --- \textbf{норма Фробениуса} (или евклидова норма для матриц), которая является аналогом $L_2$-нормы: \[ \| \mathbf{A} \|_F = \sqrt{\sum_{i=1}^N \sum_{j=1}^M |a_{ij}|^2} \] Для неё, как и для обычной $L_2$, есть масса важных и красивых свойств, которые мы скоро рассмотрим. Однако, кроме таких ``поэлементных'' норм, в линейной алгебре есть специальные \textbf{операторные (согласованные) нормы}. Они отвечают на вопрос: ``Насколько сильно матрица может растянуть вектор, если мы умножим её на него?''. Но это уже история для следующей главы, где мы вплотную займемся линейными операторами. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Операторы: когда математика начинает действовать} \subsection{Что такое оператор?} Мы подошли к одной очень важной концепции, которая сплошь и рядом встречается и в математике, и в химии. Это --- \textbf{оператор}. Оператор --- это некоторое \textit{действие} над объектом. Но тут сразу становится как-то запутанно: действие над чем? Над числом? Над вектором? Над функцией? Давай вместе рассмотрим несколько примеров, и всё встанет на свои места. \subsection{Матрица как оператор} Пусть у нас есть квадратная матрица $\mathbf{A} \in \C^{N \times N}$ (я буду далее почти всегда использовать $\C$, так как $\R$ --- это всего лишь подмножество $\C$, и всё, что работает для действительных чисел, работает и для комплексных). Рассмотрим действие $\mathbf{A}\vec{b}$, то есть умножение матрицы на вектор. Результат --- снова вектор из $\C^N$. Матрица \textit{преобразует} один вектор в другой. Это и есть простейший пример линейного оператора. \subsection{Аналогия с функциональным оператором} А теперь --- самое интересное. Аналогом матричного оператора в мире функций является \textbf{интегральный оператор}: \begin{equation} (\hat{K}g)(x) = \int_a^b K(x,y)\, g(y)\, \diff y \label{IntegralOperator} \end{equation} Здесь $K(x,y)$ --- это \textit{ядро} оператора (аналог матрицы $\mathbf{A}$), $g(y)$ --- входная функция (аналог вектора $\vec{b}$), а результат --- новая функция от $x$. Интересно, где же такая абстракция встречается в реальной жизни? Оказывается, сплошь и рядом! Например, в квантовой механике оператор Гамильтона $\hat{H}$ действует на волновую функцию $\Psi$, и это как раз интегрально-дифференциальный оператор. В спектроскопии отклик прибора на входной сигнал часто описывается интегралом свёртки --- это тоже оператор. \subsection{Возвращаемся в школу: система уравнений} Но давай начнем с самого простого. Мы еще в школе решали системы уравнений, например: \[ \begin{cases} 2x - 3y = -4 \\ 3x + 2y = 9 \end{cases} \] Ты, наверное, думаешь: ``Ну и что? Мы же умели это решать подстановкой или методом сложения''. Да, умели. Но давай перепишем это в матричном виде $\mathbf{A}\vec{x} = \vec{b}$: \[ \underbrace{ \begin{pmatrix} 2 & -3 \\ 3 & 2 \end{pmatrix} }_{\mathbf{A}} \underbrace{ \begin{pmatrix} x \\ y \end{pmatrix} }_{\vec{x}} = \underbrace{ \begin{pmatrix} -4 \\ 9 \end{pmatrix} }_{\vec{b}} \] Ты резонно спросишь меня: ``Зачем так усложнять?''. А я отвечу: нет, мы не усложняем --- мы \textit{приводим в порядок} наши знания. Мы выделяем \textit{структуру} задачи. И теперь мы можем сказать: если у нас есть \textbf{линейный оператор} $\mathbf{A}$, то он действует на вектор $\vec{x}$ и результатом является вектор $\vec{b}$. \subsection{Напоминание: базовые операции} Прежде чем идти дальше, давай зафиксируем в явном виде те операции, которыми мы будем постоянно пользоваться. \textbf{Вектор как матрица-столбец.} Вектор $\vec{x} \in \C^N$ --- это на самом деле матрица размера $N \times 1$. Поэтому все правила матричного умножения автоматически применяются и к векторам. \textbf{Транспонирование.} Операция $\mathbf{A}^T$ меняет строки и столбцы местами: $(\mathbf{A}^T)_{ij} = A_{ji}$. Для комплексных матриц чаще используют \textbf{эрмитово сопряжение} $\mathbf{A}^H = \overline{\mathbf{A}^T}$ (транспонирование + комплексное сопряжение всех элементов). \textbf{Скалярное произведение.} Для двух векторов $\vec{x}, \vec{y} \in \C^N$ скалярное произведение определяется как: \[ \langle \vec{x}, \vec{y} \rangle = \vec{x}^H \vec{y} = \sum_{k=1}^N \bar{x}_k y_k \] Для действительных векторов это просто $\vec{x}^T \vec{y} = \sum x_k y_k$. \textbf{Важнейшее свойство:} скалярное произведение вектора самого на себя равно \textbf{квадрату его евклидовой нормы}: \[ \langle \vec{x}, \vec{x} \rangle = \|\vec{x}\|_2^2 = \sum_{k=1}^N |x_k|^2 \] Это мостик между алгеброй и геометрией: длина вектора --- это корень из его ``скалярного квадрата''. \textbf{Умножение матрицы на вектор.} Если $\mathbf{A} \in \C^{M \times N}$ и $\vec{x} \in \C^N$, то результат $\vec{y} = \mathbf{A}\vec{x} \in \C^M$ вычисляется по правилу: \[ y_i = \sum_{j=1}^N A_{ij} x_j \] То есть $i$-я компонента результата --- это скалярное произведение $i$-й строки матрицы на вектор $\vec{x}$. \textbf{Умножение матриц.} Если $\mathbf{A} \in \C^{M \times K}$ и $\mathbf{B} \in \C^{K \times N}$, то произведение $\mathbf{C} = \mathbf{A}\mathbf{B} \in \C^{M \times N}$: \[ C_{ij} = \sum_{k=1}^K A_{ik} B_{kj} \] Важно: умножение матриц, вообще говоря, \textbf{не коммутативно} --- $\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}$. Это одно из самых важных отличий от умножения обычных чисел! Все эти операции --- скалярное произведение, умножение матрицы на вектор, умножение матрицы на матрицу --- имеют свои полные аналоги в мире функций и интегральных операторов. Только везде, где у нас была сумма $\sum$, появляется интеграл $\int$, а транспонирование заменяется на более сложный оператор сопряжения. Но это --- уже дело техники. \subsection{Четыре главные задачи линейной алгебры} На самом деле, пока наши операторы линейные (то есть их можно записать в виде (\ref{IntegralOperator}) или $\mathbf{A}\vec{x}$), весь мир возможных задач крутится вокруг нескольких типовых вопросов. Давай их перечислим. \textbf{Задача 1. Посчитать скалярное произведение.} Даны два вектора --- найти число $\langle \vec{x}, \vec{y} \rangle$. Это базовая операция, из которой, как из кирпичиков, строится всё остальное. \textbf{Задача 2. Посчитать действие оператора.} Даны $\mathbf{A}$ и $\vec{x}$ --- найти $\mathbf{A}\vec{x}$. Либо даны $\mathbf{A}$ и $\mathbf{B}$ --- найти $\mathbf{A}\mathbf{B}$. Это прямое вычисление. \textbf{Задача 3. Решить линейную систему уравнений.} Даны $\mathbf{A}$ и $\vec{b}$ --- найти $\vec{x}$ такой, что $\mathbf{A}\vec{x} = \vec{b}$. Это \textit{прямая} задача: известен оператор и результат, надо найти, на что он действовал. \textbf{Задача 4. Минимизировать невязку (метод наименьших квадратов).} А что, если система $\mathbf{A}\vec{x} = \vec{b}$ не имеет точного решения? (Например, уравнений больше, чем неизвестных --- переопределенная система, что типично для обработки экспериментальных данных). Тогда мы ищем такое $\vec{x}$, которое \textit{минимизирует} невязку: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_p \] Почти всегда в качестве нормы используется $p = 2$ --- это и есть знаменитый \textbf{метод наименьших квадратов (МНК)}. Он имеет глубокую геометрическую интерпретацию: мы проецируем вектор $\vec{b}$ на подпространство, натянутое на столбцы матрицы $\mathbf{A}$. \subsection{А что, если правая часть --- ноль?} А теперь --- хитрый поворот. Представим, что $\vec{b} = \vec{0}$. Тогда задача $\min_{\vec{x}} \|\mathbf{A}\vec{x}\|_2$ имеет тривиальное решение $\vec{x} = \vec{0}$. Но это же неинтересно! Давай поставим дополнительное условие: мы ищем \textbf{ненулевое} решение, например, с ограничением $\|\vec{x}\|_2 = 1$. То есть мы ищем такое направление, в котором матрица $\mathbf{A}$ сильнее всего ``сжимает'' пространство (или, наоборот, сильнее всего растягивает --- это уже вопрос знака). И вот тут на сцену выходят \textbf{собственные значения и собственные векторы}. Мы ищем такие специальные векторы $\vec{v}$ и числа $\lambda$, при которых действие матрицы сводится просто к растяжению: \[ \mathbf{A}\vec{v} = \lambda \vec{v} \] Геометрически: собственный вектор --- это такое направление, которое матрица \textit{не поворачивает}, а только растягивает или сжимает в $\lambda$ раз. Задача минимизации $\|\mathbf{A}\vec{x}\|_2$ при $\|\vec{x}\|_2 = 1$ решается именно через собственные значения $\mathbf{A}^H \mathbf{A}$ или сингулярные значения $\mathbf{A}$: минимум достигается на собственном векторе матрицы $\mathbf{A}^H \mathbf{A}$, соответствующем \textit{минимальному по модулю} собственному значению. А максимум --- на векторе, соответствующем максимальному. \subsection{А если неизвестен сам оператор?} Ты спросишь: ``А что, если нам оператор неизвестен, а мы знаем только входные функции и результаты?'' Тут ответ немного неоднозначный. Подумай: в матричном виде у нас есть только $N$ входных параметров (вектор $\vec{x}$), а мы хотим найти $N^2$ параметров матрицы $\mathbf{A}$. Природа редко устроена так, что малое число входов хорошо определяет большее число параметров. Это классическая \textbf{недоопределенная} задача. Но и тут решение есть! Если мы скажем, что одна и та же матрица $\mathbf{A}$ действует на несколько разных векторов $\vec{x}_1, \dots, \vec{x}_K$ с известными результатами $\vec{b}_1, \dots, \vec{b}_K$, то мы можем сформулировать задачу так: \begin{equation} \forall k = 1, \dots, K: \quad \mathbf{A}\vec{x}_k = \vec{b}_k, \quad \text{причём } \mathbf{A} \text{ --- неизвестна.} \label{MatApprox} \end{equation} Если собрать векторы в матрицы $\mathbf{X} = (\vec{x}_1, \dots, \vec{x}_K)$ и $\mathbf{B} = (\vec{b}_1, \dots, \vec{b}_K)$, то задача переписывается как: \[ \mathbf{A}\mathbf{X} \approx \mathbf{B} \] А это, внимание, \textbf{та же самая задача минимизации}, только теперь мы минимизируем не по $\vec{x}$, а по $\mathbf{A}$: \[ \min_{\mathbf{A}} \|\mathbf{A}\mathbf{X} - \mathbf{B}\|_F \] (Здесь $\|\cdot\|_F$ --- норма Фробениуса для матриц, о которой мы говорили в прошлой главе.) Транспонируем для удобства: $\mathbf{X}^T \mathbf{A}^T \approx \mathbf{B}^T$. И мы снова получили задачу вида ``минимизировать $\|\mathbf{M}\vec{z} - \vec{c}\|_2$'', только для каждого столбца $\mathbf{A}^T$ отдельно. Это и есть классическая \textbf{линейная регрессия} --- основа хемометрики, QSAR, калибровки приборов и многого другого. \subsection{Сингулярное разложение: мостик в большие миры} И вот тут на сцену выходит одна из самых красивых конструкций во всей линейной алгебре --- \textbf{сингулярное разложение (SVD, Singular Value Decomposition)}. Любую матрицу $\mathbf{A} \in \C^{M \times N}$ можно представить в виде: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] где $\mathbf{U} \in \C^{M \times M}$ и $\mathbf{V} \in \C^{N \times N}$ --- унитарные матрицы (их столбцы --- ортонормированные векторы), а $\mathbf{\Sigma}$ --- диагональная матрица с неотрицательными числами $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ на диагонали. Эти числа называются \textbf{сингулярными значениями}. Геометрический смысл SVD потрясающ: любая матрица --- это просто поворот ($\mathbf{V}^H$), потом растяжение вдоль осей ($\mathbf{\Sigma}$), потом еще один поворот ($\mathbf{U}$). Всё! Больше матрицы ничего делать не умеют. SVD --- это универсальный инструмент. Через него решаются: \begin{itemize} \item задачи наименьших квадратов, \item поиск псевдообратной матрицы, \item сжатие данных и выделение главных компонент (PCA --- Principal Component Analysis, основа хемометрики!), \item регуляризация плохо обусловленных задач. \end{itemize} Здесь есть важный момент, унитарные матрицы, имеют одно очень важное свойство, что всегда $\mathbf{V}^H \mathbf{V} = \mathbf{I}$, где $\mathbf{I}$ --- единичная матрица, то есть такая, у которой на главной диагонали стоят единички, а все остальное заполнено нулями. \subsection{Мостик в функциональный анализ: теория Фредгольма} А теперь --- самое интересное. Помнишь интегральный оператор (\ref{MatApprox})? Так вот, для него тоже есть аналог SVD! Это называется \textbf{сингулярное разложение компактного оператора} или, в более классической формулировке, \textbf{теория Фредгольма}. Суть в следующем: ядро интегрального оператора $K(x,y)$ можно разложить в ряд по ``сингулярным функциям'' $u_k(x)$ и $v_k(y)$: \[ K(x,y) = \sum_{k=1}^\infty \sigma_k \, u_k(x) \overline{v_k(y)} \] где $\sigma_k$ --- сингулярные значения (убывающие к нулю), а $u_k, v_k$ --- ортонормированные системы функций. Это в точности аналог SVD для матриц, только в бесконечномерном пространстве! И все идеи, которые мы поняли на матрицах, переносятся сюда почти дословно. Но не будем перегружать сознание Фредгольмом и функциональным анализом прямо сейчас. Хорошая новость в том, что \textbf{большую часть всего, что нам нужно в химии и обработке сигналов, можно сделать на матричном уровне}. А там, где уже не получается (например, в квантовой механике или в строгой теории интегральных уравнений), мы просто вспомним, что ``где-то там есть Фредгольм'', и при необходимости спросим детали у ИИ --- он с радостью расскажет про теорию Фредгольма, про ядра Гильберта--Шмидта и про спектры компактных операторов. Главное --- понять \textit{структуру}. А структура везде одна и та же: оператор, пространство, действие, разложение по ``базисным направлениям''. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Обусловленность: когда матрица ``хочет'' нас обмануть} \subsection{Когда решение есть, но его как будто нет} Итак, пусть у нас есть квадратная матрица $\mathbf{A}$ и линейная система $\mathbf{A}\vec{x} = \vec{b}$. Что мы можем про неё сказать? Мы помним из школьного курса, что иногда система уравнений бывает такой, что когда мы начинаем подстановку, вдруг у нас либо нет решений, либо мы получаем в конце два и более одинаковых уравнения --- и тогда решений получается бесконечно много. В химии, да и в любой индустриальной задаче, часто исходные данные для матричных уравнений приходят из измерений с приборов. И тут может не повезти сразу тремя способами: \begin{enumerate} \item Наша система будет вырожденной (нет единственного решения). \item Наша система будет иметь множество решений. \item \textbf{Самый коварный случай:} мы, решая в машинной арифметике, окажемся \textit{очень близко} к плохому случаю, но этого не заметим. И получим ответ, который выглядит разумно, но на самом деле --- полная ерунда. \end{enumerate} Когда же это может произойти? Давай разберёмся на живом примере. \subsection{Пример из хроматографии: ловушка большого и малого} Пусть у нас есть две хроматограммы, в которых были сняты два вещества на фоне растворителя. Пик растворителя вышел очень широким и ``заполз'' на пики искомых веществ, а отклики искомых веществ оказались очень слабенькими. Фактически, в пике искомых веществ мы имеем: \begin{itemize} \item Первая хроматограмма: $(a + b_1)$, где $a$ --- огромный отклик от растворителя, а $b_1$ --- очень маленькое число от искомого вещества. \item Вторая хроматограмма: $(a + b_2)$, но так как температура чуть-чуть стала выше, то это измерение слегка неточно, скажем, повысилось на какое-то значение $\varepsilon$ \textit{относительно пика растворителя}, то есть $(1+\varepsilon)(a+b_2)$. \end{itemize} Если мы захотим вычесть из первого второе, чтобы посчитать, чему равно $b_1 - b_2$, то вместо этого мы получим: \[ (a+b_1) - (1+\varepsilon)(a+b_2) = b_1 - b_2 - \varepsilon(a+b_2) \approx b_1 - b_2 - \varepsilon a \] То есть если $\varepsilon a$ сравнимо с $b_1 - b_2$ по порядку, мы можем получить не просто большую ошибку, но и \textbf{противоположный знак}! Это и есть суть плохой обусловленности: когда в одном числе ``живут'' одновременно очень большая и очень маленькая компоненты, и мы пытаемся извлечь маленькую через вычитание. \begin{warningbox}[Правило вычитания близких чисел] Никогда не вычитай два близких по величине числа, если нужна относительная точность. Потеря значащих цифр --- это не баг компьютера, это фундаментальное свойство арифметики. \end{warningbox} \subsection{Как SVD раскрывает механизм катастрофы} То же самое происходит и при решении систем уравнений. И у нас есть классный способ \textit{заранее} оценить, что же нам делать. Пусть у нас есть система $\mathbf{A}\vec{x} = \vec{b}$, где мы ищем $\vec{x}$, а всё остальное дано. Вспомним сингулярное разложение $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$. Тогда решение можно переписать в несколько стадий: \begin{align*} \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \vec{x} &= \vec{b} \\ \mathbf{\Sigma} \mathbf{V}^H \vec{x} &= \mathbf{U}^H \vec{b} \equiv \vec{t}_1 \\ \mathbf{V}^H \vec{x} &= \mathbf{\Sigma}^{-1} \vec{t}_1 \equiv \vec{t}_2 \\ \vec{x} &= \mathbf{V} \vec{t}_2 \end{align*} Здесь мы используем свойства того, что решать систему с унитарными матрицами $\mathbf{U}, \mathbf{V}$ очень просто --- достаточно домножить на $\mathbf{U}^H$ или $\mathbf{V}^H$, так как $\mathbf{U}^H \mathbf{U} = \mathbf{I}$. А решать систему с диагональной матрицей --- это вообще тривиально. Давай нарисуем, как это выглядит для матрицы $3 \times 3$: \[ \begin{pmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 & \sigma_3 \end{pmatrix} \begin{pmatrix} t_{2,1} \\ t_{2,2} \\ t_{2,3} \end{pmatrix} = \begin{pmatrix} t_{1,1} \\ t_{1,2} \\ t_{1,3} \end{pmatrix} \quad \Rightarrow \quad \begin{cases} \sigma_1 t_{2,1} = t_{1,1} \\ \sigma_2 t_{2,2} = t_{1,2} \\ \sigma_3 t_{2,3} = t_{1,3} \end{cases} \quad \Rightarrow \quad t_{2,k} = \frac{t_{1,k}}{\sigma_k} \] Видите? Каждое уравнение --- это просто деление одной компоненты на соответствующее сингулярное число. И вот тут --- самое важное. Представим, что мы посчитали сингулярные числа исходной матрицы и заметили, что самое большое сингулярное число $\sigma_1$ очень-очень сильно больше самого маленького $\sigma_N$. Тогда на шаге $\vec{t}_2 = \mathbf{\Sigma}^{-1} \vec{t}_1$ происходит следующее: \begin{itemize} \item В векторе $\vec{t}_1$ ошибка измерений распределена более-менее равномерно по всем компонентам. \item После умножения на $\mathbf{\Sigma}^{-1}$ компонента, соответствующая \textit{маленькому} $\sigma_N$, возрастает в $\sigma_1/\sigma_N$ раз! \item А если $\sigma_N$ вообще равен нулю? Тогда мы делим на ноль --- и решение не существует. \end{itemize} Именно из-за этого люди решили обозначать отношение $\sigma_1 / \sigma_N$ как \textbf{число обусловленности} матрицы: \[ \cond(\mathbf{A}) = \frac{\sigma_{\max}}{\sigma_{\min}} \] Фактически, это --- характеристика того, \textit{насколько мы можем испортить решение такой матрицей}. Если $\cond(\mathbf{A}) \approx 10^k$, то при решении системы мы теряем примерно $k$ значащих цифр точности. Для double precision (16 значащих цифр) это означает, что при $\cond(\mathbf{A}) \approx 10^{16}$ ответ будет полностью состоять из шума. \begin{tipbox}[Важное замечание] Число обусловленности --- это характеристика \textit{матрицы}, а не алгоритма. Никакой, даже самый гениальный алгоритм не может решить плохо обусловленную систему точнее, чем позволяет $\cond(\mathbf{A})$. Это фундаментальное ограничение задачи, а не наших вычислительных методов. При этом обусловленность \textit{не влияет} на вычисление собственных векторов симметричных/эрмитовых матриц или на поиск сингулярных векторов --- эти задачи, как правило, гораздо устойчивее. \end{tipbox} \subsection{Где это всё встречается в химии и за её пределами?} Что-то у нас было много сухой теории, и надо бы вспомнить, где же это всё в химии применяется. \textbf{Квантовая химия: уравнение Шрёдингера.} Да, волновое уравнение Шрёдингера и все его приближения --- уравнения Хартри--Фока, Density Functional Theory (DFT) --- всё это задачи на собственные значения: поиск таких векторов $\vec{x}$ и значений $\lambda$, которые удовлетворяют уравнению $\mathbf{A}\vec{x} = \lambda \vec{x}$. Матрицы здесь --- матрицы Фока или матрицы Гамильтониана --- часто имеют размерность тысяч и десятков тысяч, и их обусловленность напрямую определяет, сможем ли мы вообще получить физически осмысленный ответ. \textbf{Масс-спектрометрия и спектроскопия.} Поиск снятого спектра по базе данных откликов в масс-спектрометре --- это алгоритмы, которые базируются на решении линейных систем. Если у вас есть смесь веществ, и вы хотите понять, из чего она состоит, вы решаете систему $\mathbf{A}\vec{x} = \vec{b}$, где $\mathbf{A}$ --- матрица эталонных спектров, $\vec{b}$ --- измеренный спектр смеси, а $\vec{x}$ --- концентрации компонентов. И если спектры компонентов похожи --- матрица плохо обусловлена, и концентрации будут ``прыгать''. \textbf{Крио-электронная микроскопия.} Реконструкция 3D-структуры белков из тысяч 2D-проекций --- это гигантская разреженная система уравнений. Обусловленность там --- ключевой фактор, определяющий разрешение итоговой структуры. \textbf{ЯМР-спектроскопия.} Обращение свёртки (деконволюция) --- это классическая плохо обусловленная задача. Именно поэтому в ЯМР используют регуляризацию Тихонова и преобразование Фурье --- они превращают плохо обусловленную задачу в диагональную. \textbf{Калибровка аналитических приборов.} PLS-регрессия (Partial Least Squares) --- это, по сути, SVD с регуляризацией. Обусловленность матрицы спектров напрямую влияет на точность предсказания концентраций. \textbf{Поиск в интернете (PageRank).} Даже очень-очень старый Google-поиск был устроен на том, что все заранее найденные документы классифицировались в огромную матрицу, и для неё считалось сингулярное разложение (или, что эквивалентно, ищется главный собственный вектор матрицы переходов), что помогало моментально найти искомый документ по совпадающей комбинации слов. В современном ИИ SVD используется очень-очень часто --- от сжатия весов нейросетей до понижения размерности в рекомендательных системах. \textbf{Обработка сигналов и изображений.} Деконволюция размытых изображений в микроскопии, подавление шума в ЭКГ/ЭЭГ, сжатие аудио (MP3) и видео --- всё это задачи, в которых обусловленность играет ключевую роль. Когда-то эти три задачи --- решение линейных систем, поиск собственных значений и SVD --- были камнем преткновения при решении больших индустриальных задач. Даже были фирмы, которые \textit{только} разрабатывали такие солверы --- и ничегошеньки больше. Правда, это было ещё в 90-ые годы прошлого века. \subsection{Сколько это стоит? Вычислительная сложность} Грубо говоря, если у нас есть квадратная матрица $N \times N$ и мы решаем либо линейную систему с ней, либо ищем её собственные или сингулярные числа и вектора, то нам надо потратить примерно $\mathcal{O}(N^3)$ арифметических операций, если мы никак специально не используем структуру или свойства этой матрицы. Это --- ``цена'' полной задачи. Но если матрица обладает специальными свойствами, цену можно сильно снизить --- иногда до $\mathcal{O}(N)$ или $\mathcal{O}(N \log N)$. И об этом --- ниже. % === СВЯЗ SVD И СОБСТВЕННЫХ ЗНАЧЕНИЙ === \subsection{Связь сингулярного разложения и собственных значений} А теперь --- один из самых красивых фактов во всей линейной алгебре. Пусть у нас есть сингулярное разложение: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] Давайте посмотрим, что получится, если мы умножим $\mathbf{A}$ на её эрмитово сопряжённую $\mathbf{A}^H$: \begin{align*} \mathbf{A} \mathbf{A}^H &= (\mathbf{U} \mathbf{\Sigma} \mathbf{V}^H)(\mathbf{V} \mathbf{\Sigma} \mathbf{U}^H) \\ &= \mathbf{U} \mathbf{\Sigma} (\mathbf{V}^H \mathbf{V}) \mathbf{\Sigma} \mathbf{U}^H \\ &= \mathbf{U} \mathbf{\Sigma}^2 \mathbf{U}^H \end{align*} Здесь мы использовали, что $\mathbf{V}^H \mathbf{V} = \mathbf{I}$ (унитарность $\mathbf{V}$), и $\mathbf{\Sigma}^T = \mathbf{\Sigma}$ (диагональная матрица). Аналогично: \begin{align*} \mathbf{A}^H \mathbf{A} &= (\mathbf{V} \mathbf{\Sigma} \mathbf{U}^H)(\mathbf{U} \mathbf{\Sigma} \mathbf{V}^H) \\ &= \mathbf{V} \mathbf{\Sigma} (\mathbf{U}^H \mathbf{U}) \mathbf{\Sigma} \mathbf{V}^H \\ &= \mathbf{V} \mathbf{\Sigma}^2 \mathbf{V}^H \end{align*} \textbf{Что это означает?} \begin{itemize} \item Сингулярные векторы матрицы $\mathbf{A}$ --- это \textit{в точности} собственные векторы матриц $\mathbf{A}\mathbf{A}^H$ и $\mathbf{A}^H\mathbf{A}$. \item Сингулярные числа $\sigma_k$ матрицы $\mathbf{A}$ --- это квадратные корни из собственных чисел $\lambda_k$ матриц $\mathbf{A}\mathbf{A}^H$ или $\mathbf{A}^H\mathbf{A}$: \[ \sigma_k = \sqrt{\lambda_k} \] \end{itemize} Это --- глубокая связь между двумя, казалось бы, разными задачами: сингулярным разложением и поиском собственных значений. \begin{warningbox}[Осторожно: квадрат обусловленности!] Но есть одна коварная деталь. Число обусловленности матриц $\mathbf{A}\mathbf{A}^H$ и $\mathbf{A}^H\mathbf{A}$ --- это \textbf{квадрат} числа обусловленности исходной матрицы: \[ \cond(\mathbf{A}\mathbf{A}^H) = \cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2 \] Если $\cond(\mathbf{A}) = 10^8$, то $\cond(\mathbf{A}^H\mathbf{A}) = 10^{16}$ --- и мы потеряем все 16 значащих цифр точности! Поэтому переходить от задачи сингулярного разложения к задаче собственных значений для $\mathbf{A}^H\mathbf{A}$ надо \textbf{крайне осторожно}. В современных библиотеках (LAPACK) SVD считается напрямую, без явного построения $\mathbf{A}^H\mathbf{A}$ --- именно чтобы избежать этой катастрофы. \end{warningbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Классификация матриц: кто есть кто} Теперь, когда понятно, что практически всё базируется на этих ``простых'' математических задачах, давай немного систематизируем наши знания о таких решениях. Ведь всегда у нас есть матрица, и от свойств этой матрицы зависит очень многое. \subsection{По форме} \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Тип} & \textbf{Размер} & \textbf{Описание} \\ \midrule Квадратная & $N \times N$ & Одинаковое число строк и столбцов. Только для таких матриц определён определитель, собственные значения и обратная матрица. \\ ``Стоячая'' (высокая) & $M \times N, \; M > N$ & Уравнений больше, чем неизвестных. Обычно не имеет точного решения --- решаем через МНК. \\ ``Лежачая'' (широкая) & $M \times N, \; M < N$ & Неизвестных больше, чем уравнений. Решений бесконечно много --- ищем минимальное по норме. \\ \bottomrule \end{tabularx} \subsection{По симметрии и структуре элементов} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Тип} & \textbf{Определение и свойства} & \textbf{Сложность} \\ \midrule Симметричная действительная & $\mathbf{A} = \mathbf{A}^T$, элементы $\in \R$. Собственные значения \textbf{действительные}, собственные векторы --- ортогональны. & $\mathcal{O}(N^3/3)$ \\ Эрмитова комплексная & $\mathbf{A} = \mathbf{A}^H$, элементы $\in \C$. Собственные значения \textbf{действительные}, собственные векторы --- ортонормированы. & $\mathcal{O}(4N^3/3)$ \\ Кососимметричная & $\mathbf{A} = -\mathbf{A}^T$. Собственные значения --- чисто мнимые или нули. & $\mathcal{O}(N^3/3)$ \\ Ортогональная / Унитарная & $\mathbf{A}^T\mathbf{A} = \mathbf{I}$ / $\mathbf{A}^H\mathbf{A} = \mathbf{I}$. Сохраняет длины и углы. Собственные значения по модулю равны 1. & $\mathcal{O}(N^2)$ для умножения \\ Единичная & $\mathbf{A} = \mathbf{I}$. Диагональ из единиц, остальное --- нули. & $\mathcal{O}(N)$ \\ Диагональная & $A_{ij} = 0$ при $i \neq j$. Хранится как вектор из $N$ чисел. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} \subsection{По положительной определённости} Это --- подкласс симметричных/эрмитовых матриц, и он критически важен для оптимизации и статистики. \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Тип} & \textbf{Определение и свойства} \\ \midrule Положительно определённая ($\mathbf{A} \succ 0$) & $\forall \vec{x} \neq 0: \; \vec{x}^H \mathbf{A} \vec{x} > 0$. Все собственные значения строго положительны. Обусловленность --- отношение максимального и минимального собственных значений. Решается методом Холецкого за $\mathcal{O}(N^3/3)$. \\ Неотрицательно определённая ($\mathbf{A} \succeq 0$) & $\forall \vec{x}: \; \vec{x}^H \mathbf{A} \vec{x} \ge 0$. Все собственные значения $\ge 0$. Может быть вырожденной. \\ \bottomrule \end{tabularx} \subsection{По структуре расположения ненулевых элементов} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Тип} & \textbf{Определение и свойства} & \textbf{Сложность} \\ \midrule Ленточная (banded) & Ненулевые элементы только вблизи главной диагонали: $A_{ij} = 0$ при $|i-j| > k$. Хранится как $N \times (2k+1)$. & $\mathcal{O}(N k^2)$ \\ Теплицева & $A_{ij}$ зависит только от $i-j$. Каждая диагональ --- константа. Возникает в задачах со стационарными процессами. & $\mathcal{O}(N^2)$, через БПФ --- $\mathcal{O}(N \log N)$ \\ Циркулянтная & Теплицева + периодичность: каждая строка --- циклический сдвиг предыдущей. Диагонализуется преобразованием Фурье. & $\mathcal{O}(N \log N)$ через БПФ \\ Блочная & Состоит из блоков, каждый из которых --- матрица. Позволяет использовать рекурсивные алгоритмы. & Зависит от структуры блоков \\ Разреженная (sparse) & Большинство элементов --- нули. Хранится в специальных форматах (CSR, CSC). Основа больших вычислений. & $\mathcal{O}(\text{nnz})$, где nnz --- число ненулевых \\ \bottomrule \end{tabularx} \subsection{Вырожденность} \textbf{Вырожденная (сингулярная) матрица} --- это матрица, определитель которой равен нулю, или, что эквивалентно, у которой хотя бы одно сингулярное значение равно нулю. Для такой матрицы: \begin{itemize} \item не существует обратной матрицы, \item система $\mathbf{A}\vec{x} = \vec{b}$ либо не имеет решений, либо имеет бесконечно много, \item $\rank(\mathbf{A}) < N$. \end{itemize} На практике матрицы почти никогда не бывают \textit{точно} вырожденными --- но они могут быть \textit{почти} вырожденными, то есть с очень маленьким $\sigma_{\min}$. И вот это --- гораздо более коварный случай, потому что формально решение существует, но оно полностью определяется шумом в данных. % === ДОПОЛНЕНИЕ К КЛАССИФИКАЦИИ МАТРИЦ === \subsection{Дополнительные важные типы матриц} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Тип} & \textbf{Определение и свойства} & \textbf{Сложность} \\ \midrule Верхнетреугольная & $A_{ij} = 0$ при $i > j$. Все ненулевые элементы выше или на главной диагонали. Решение системы --- обратная подстановка. & $\mathcal{O}(N^2)$ \\ Нижнетреугольная & $A_{ij} = 0$ при $i < j$. Все ненулевые элементы ниже или на главной диагонали. Решение системы --- прямая подстановка. & $\mathcal{O}(N^2)$ \\ Матрица перестановок & В каждой строке и каждом столбце ровно одна единица, остальное --- нули. Умножение на такую матрицу --- это перестановка строк/столбцов. $\mathbf{P}^T = \mathbf{P}^{-1}$. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} Матрицы перестановок --- это не просто абстракция. Они критически важны для численной устойчивости алгоритмов (например, в LU-разложении с выбором главного элемента), и мы ещё встретимся с ними. \subsection{Сводная таблица: что и как решать} \begin{longtable}{l c c c} \toprule \textbf{Задача} & \textbf{Общий случай} & \textbf{Симметричная} & \textbf{Разреженная} \\ \midrule Линейная система $\mathbf{A}\vec{x} = \vec{b}$ & LU, $\mathcal{O}(N^3)$ & Холецкий, $\mathcal{O}(N^3/3)$ & Итерационные, $\mathcal{O}(\text{nnz} \cdot k)$ \\ МНК $\min \|\mathbf{A}\vec{x} - \vec{b}\|_2$ & QR, $\mathcal{O}(MN^2)$ & --- & LSQR, итерационные \\ Собственные значения & $\mathcal{O}(N^3)$ & $\mathcal{O}(N^3/3)$, все действительные & Lanczos, $\mathcal{O}(\text{nnz} \cdot k)$ \\ SVD & $\mathcal{O}(MN^2)$ & Через собственные $\mathbf{A}^T\mathbf{A}$ & Усечённое SVD, итерационные \\ \bottomrule \end{longtable} Здесь $k$ --- число итераций, которое зависит от желаемой точности и обусловленности. \begin{tipbox}[Главный вывод] Прежде чем решать задачу, \textit{посмотри на матрицу}. Её свойства --- симметрия, разреженность, ленточная структура --- могут снизить сложность с кубической до линейной. А её обусловленность скажет тебе, стоит ли вообще доверять полученному ответу. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Как решают задачи линейной алгебры: от теории к практике} \subsection{Одного решения на все случаи --- нет} Так как же решают эти задачи линейной алгебры? Наверное, есть одно или может быть пара самых простых и надёжных решений? К сожалению, даже сейчас одного решения на все случаи жизни --- нет и не предвидится. Автор этого скрипта участвовал в своё время в написании более сотни различных таких алгоритмов. Да, таких алгоритмов --- очень много, и надо хотя бы коротко, пусть даже ни разу не имплементировав это, понимать, когда и что можно использовать. \subsection{Классификация методов решения} Методы решения можно классифицировать примерно так: \subsubsection{Прямые методы разложения} \textbf{LU-разложение:} $\mathbf{A} = \mathbf{L}\mathbf{U}$, где $\mathbf{L}$ --- нижнетреугольная, $\mathbf{U}$ --- верхнетреугольная матрица. С такой факторизацией решение системы $\mathbf{A}\vec{x} = \vec{b}$ сводится к двум простым подстановкам: \begin{enumerate} \item Решаем $\mathbf{L}\vec{y} = \vec{b}$ (прямая подстановка, $\mathcal{O}(N^2)$) \item Решаем $\mathbf{U}\vec{x} = \vec{y}$ (обратная подстановка, $\mathcal{O}(N^2)$) \end{enumerate} Но в общем виде $\cond(\mathbf{L}) \cdot \cond(\mathbf{U}) \ge \cond(\mathbf{A})$, а если не делать перестановок, то произведение обусловленностей может стать существенно больше $\cond(\mathbf{A})$. То есть мы можем испортить нашу исходную задачу! Поэтому на практике почти всегда используют \textbf{LUP-разложение}: $\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}$, где $\mathbf{P}$ --- матрица перестановок. Перестановки выбираются так, чтобы на диагонали $\mathbf{U}$ стояли максимально возможные элементы (выбор главного элемента --- pivoting). Это гарантирует численную устойчивость. \textbf{Для специальных матриц:} \begin{itemize} \item Если $\mathbf{A}$ --- ленточная, то $\mathbf{L}$ и $\mathbf{U}$ не выходят за пределы ленты. Сложность $\mathcal{O}(N k^2)$, где $k$ --- ширина ленты. \item Если $\mathbf{A}$ --- разреженная, то $\mathbf{L}$ и $\mathbf{U}$ могут содержать существенно больше ненулевых элементов (fill-in), и это может быть большой проблемой. \item Для симметричной матрицы --- $\mathbf{A} = \mathbf{L}\mathbf{D}\mathbf{L}^H$ или $\mathbf{P}\mathbf{A}\mathbf{P}^H = \mathbf{L}\mathbf{D}\mathbf{L}^H$. \item Для положительно определённой --- $\mathbf{A} = \mathbf{L}\mathbf{L}^H$, так называемый \textbf{метод Холецкого}. Но проблема роста числа обусловленности есть даже для положительно определённых матриц. \end{itemize} \textbf{QR-разложение:} $\mathbf{A} = \mathbf{Q}\mathbf{R}$, где $\mathbf{Q}$ --- унитарная ($\mathbf{Q}^H\mathbf{Q} = \mathbf{I}$), $\mathbf{R}$ --- верхнетреугольная. Этот метод \textbf{гарантирует не увеличение числа обусловленности} решения, так как унитарные преобразования сохраняют длины векторов. Хорош для плотных матриц без структуры. Но если исходная матрица --- разреженная, то $\mathbf{Q}$ практически всегда будет плотной (только для очень специальных случаев и методов), что сильно ограничивает применимость этого метода для решения огромных задач, когда сама матрица помещается в память, а её плотная версия $N^2$ --- не помещается. \textbf{Разложение Шура:} $\mathbf{A} = \mathbf{Q}\mathbf{G}\mathbf{Q}^H$, где $\mathbf{G}$ --- верхнетреугольная (для действительных матриц --- верхняя квазитреугольная с блоками $2 \times 2$ на диагонали). Применяется для поиска собственных значений и собственных векторов, которые находят, решая задачу на собственные значения для уже треугольной матрицы $\mathbf{G}$ (а для треугольной матрицы собственные значения --- это просто диагональные элементы!). \subsubsection{Итерационные методы} Это методы решения линейной системы или задачи на собственные значения, когда мы можем эффективно умножить на матрицу из-за её структуры, и строим решение итерационно. Тут есть много методов, но важно заметить: \textbf{сходимость итерационных методов зависит от числа обусловленности}. Чем оно больше, тем медленнее они сходятся. (Есть исключения: если в матрице есть только несколько сильно больших и сильно маленьких сингулярных чисел, тогда сходимость может быть существенно быстрее.) Все эти методы можно также разделить на подклассы: \begin{enumerate} \item \textbf{Специальные методы для положительно определённых матриц} с маленькими затратами по дополнительной памяти и разумным числом итераций. Самый известный --- \textbf{метод сопряжённых градиентов (Conjugate Gradients, CG)}. \item \textbf{Приведение несимметричной задачи к симметричной:} вместо решения $\mathbf{A}\vec{x} = \vec{b}$ решать $\mathbf{A}^H\mathbf{A}\vec{x} = \mathbf{A}^H\vec{b}$. Если число обусловленности $\cond(\mathbf{A}^H\mathbf{A})$ не настолько огромное по сравнению с $\cond(\mathbf{A})$, то это разумно. Но помните: $\cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2$, так что это палка о двух концах. \item \textbf{Предобуславливатели (preconditioners):} такие специальные матрицы $\mathbf{P} \approx \mathbf{A}^{-1}$, что вместо решения $\mathbf{A}\vec{x} = \vec{b}$ мы решаем $\mathbf{P}\mathbf{A}\vec{x} = \mathbf{P}\vec{b}$, предполагая, что $\cond(\mathbf{P}\mathbf{A}) \ll \cond(\mathbf{A})$, что сильно ускоряет решение такими итерационными методами. Кстати, часто для больших разреженных матриц предобуславливатели строят как \textbf{неполное LU-разложение}: то есть строят для исходной матрицы $\mathbf{A} \approx \mathbf{L}\mathbf{U}$, но стараются ``выбросить'' новые ненулевые элементы при построении $\mathbf{L}\mathbf{U}$. В таком случае применять эту матрицу как что-то похожее на исходную ещё можно, поэтому, решая с ней, мы предобуславливаем исходную задачу и повышаем сходимость. \end{enumerate} \subsection{Два мира: плотные и разреженные задачи} На самом деле, большинство алгоритмов решения можно разделить на два класса: \begin{enumerate} \item Когда мы готовы потратить $\mathcal{O}(N^3)$ арифметических операций и у нас есть $\mathcal{O}(N^2)$ памяти. \item Когда нам надо ужаться, а матрица настолько огромная и имеет какую-то специальную структуру, что мы на неё можем умножать, но взять её в полном виде --- очень сложно. \end{enumerate} В первом случае --- это \textbf{методы разложения} (LU, QR, Холецкий). Во втором --- это \textbf{итерационные методы}. Причём из-за специфичной структуры устройства современных процессоров и памяти, итерационные методы имеют немного худшую производительность, чем методы полного разложения. Поэтому применять итерационные методы для очень маленьких матриц практически всегда не бывает оправдано. \subsection{Примеры из квантовой химии: DFT} \textbf{Пример 1: Небольшая система.} Небольшая система на несколько электронов методом DFT, и размер гамильтониана --- около $1000 \times 1000$. Мы вписываемся в память (это всего $\sim 8$ МБ для double precision). Мы можем найти все собственные значения и набор тех собственных векторов, которые нас интересуют, и для этого очевидный выбор --- прямой метод (например, разложение Шура или QR-алгоритм). \textbf{Пример 2: Большая химическая структура.} Довольно большая химическая структура методом DFT, в которой есть около тысячи электронных пар на внешних орбиталях. Для этого у нас построился гамильтониан уравнения Шрёдингера, и нам надо найти на каждую электронную пару по своему собственному вектору. А гамильтониан у нас задан так, что мы взяли под сотню базисных функций на каждую орбиталь, то есть размерность этого гамильтониана --- около $100\,000 \times 100\,000$. Это вроде бы и не очень много, но только сама матрица такого гамильтониана уже занимает в оперативной памяти примерно \textbf{80 ГБ}, и само решение будет ещё почти гигабайт занимать. Тут гораздо очевиднее итерационное решение (например, метод Ланцоша или Davidson), которое находит только несколько нужных собственных векторов, не работая с полной матрицей. \subsection{Не изобретайте велосипед: библиотеки} Сейчас во многих языках программирования есть хорошо и годами вылизанные библиотеки по решению таких систем уравнений, поиску сингулярных чисел и векторов и собственных векторов и чисел. Достаточно поискать в документации на \textbf{NumPy} для Python, и в \textbf{LAPACK/BLAS} для C/C++ --- и всё будет сразу понятно. Более того, код очень хорошо оптимизирован под современные компьютеры и довольно сложный. Например, современная версия SVD в LAPACK содержит около \textbf{полумиллиона строк кода}, и повторить его или даже сделать лучше --- это реально почти невыполнимая задача. \begin{tipbox}[Практический совет] Никогда не пишите свои солверы для линейной алгебры, если только это не учебная задача. Используйте: \begin{itemize} \item \textbf{Python:} NumPy (\texttt{numpy.linalg}), SciPy (\texttt{scipy.linalg}) \item \textbf{C/C++:} LAPACK, BLAS, Eigen \item \textbf{Fortran:} LAPACK (это его родная стихия) \item \textbf{MATLAB:} встроенные функции (\texttt{eig}, \texttt{svd}, \texttt{lu}, \texttt{qr}) \end{itemize} Эти библиотеки прошли десятилетия оптимизации и тестирования. Они знают про кэш-процессора, про SIMD-инструкции, про многопоточность --- всё то, что вы не сможете повторить вручную. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Нелинейные операторы: когда мир перестает быть прямым} \subsection{Что такое нелинейный оператор?} До сих пор мы говорили о линейных операторах --- тех, которые можно записать как $\mathbf{A}\vec{x}$ или $\int K(x,y)g(y)\diff y$. Но реальный мир --- нелинеен. \textbf{Нелинейный оператор} --- это такое отображение $\mathcal{F}$, которое действует на функцию или вектор $\vec{x}$, но которое \textit{не} обладает свойствами аддитивности и однородности: \[ \mathcal{F}(\alpha \vec{x} + \beta \vec{y}) \neq \alpha \mathcal{F}(\vec{x}) + \beta \mathcal{F}(\vec{y}) \] Грубо говоря, дефинируем теперь, что у нас есть какая-то функция $f(\vec{x})$, которая зависит от одного или нескольких параметров. И мы хотим либо: \begin{itemize} \item Решить уравнение $f(\vec{x}) = 0$ (систему нелинейных уравнений), \item Найти минимум $\min_{\vec{x}} f(\vec{x})$ (задача оптимизации). \end{itemize} \subsection{Пример из хроматографии: модель тарелок} Классический пример --- система уравнений баланса и констант диссоциации на каждой хроматографической тарелке. Мы хотим моделировать хроматограф, вернее его колонку. В ней идет разделение, и подвижная фаза, перемещаясь вдоль колонки, на каждом шаге фактически создает такое равновесие. Вы скажете --- ой, да тут всего какие-то 4--5 уравнений, и столько же неизвестных! Да, согласен, немного. Но: \begin{enumerate} \item Тарелок-то много (пусть тысяча) --- это раз. \item Мы разбиваем время прохождения вещества по колонке на множество шагов по времени --- это два. \end{enumerate} То есть на каждом шаге по времени мы имеем тысячу раз одну и ту же систему уравнений с разными параметрами входных концентраций. И таких шагов по времени будет существенно больше, чем число тарелок, ведь вещество-то всё-таки колонкой удерживается. То есть нам надо решить эдак под \textbf{десять миллионов раз} систему уравнений из всего-то 5 уравнений с 5 неизвестными. Но стоит хоть один раз не решить --- мы не сможем получить ни одно дальнейшее решение. К сожалению, такая система --- не является линейной. Уравнения баланса --- линейны, но вот уравнения связи констант диссоциации содержат произведения и отношения неизвестных друг на друга. И, на удивление, такая система иногда может действительно очень плохо решаться. \subsection{Классификация нелинейных задач} Прежде чем говорить о методах, давайте классифицируем сами задачи: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Свойство} & \textbf{Описание} \\ \midrule \textbf{Гладкость:} & \\ \quad Гладкие & Функция имеет непрерывные производные (как минимум первую, а желательно и вторую). Пример: $f(x) = x^2 + \sin x$. \\ \quad Почти гладкие & Функция гладкая почти везде, но есть точки излома или разрывы производных. Пример: $f(x) = |x|$. \\ \quad Разрывные & Функция имеет разрывы или очень шумная. Пример: экспериментальные данные с артефактами. \\ \midrule \textbf{Число минимумов:} & \\ \quad Унимодальные & Есть только один глобальный минимум (или максимум). Пример: выпуклые функции. \\ \quad Мультимодальные & Есть несколько локальных минимумов, и задача --- найти глобальный. Пример: потенциалы сложных молекул. \\ \bottomrule \end{tabularx} \subsection{Методы решения нелинейных задач} Теперь давайте систематизируем методы решения. Каждый метод что-то требует, от чего-то зависит, и у каждого есть свои сильные и слабые стороны. \subsubsection{Методы Монте-Карло} \textbf{Идея:} Случайный поиск. Генерируем случайные точки в пространстве параметров, вычисляем в них функцию, выбираем лучшую. \textbf{Требования:} Ничего не требует. Может работать даже с разрывными функциями. \textbf{Плюсы:} \begin{itemize} \item Не застревает в локальных минимумах (при достаточном числе итераций), \item Прост в реализации, \item Параллелится идеально. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Очень медленная сходимость: ошибка убывает как $\mathcal{O}(1/\sqrt{N})$, где $N$ --- число итераций, \item Не использует информацию о структуре функции. \end{itemize} \textbf{Когда использовать:} Когда функция очень шумная, разрывная, или когда нужно просто найти ``хоть какое-то'' решение для инициализации других методов. \subsubsection{Градиентные методы (Gradient Descent)} \textbf{Идея:} Двигаться в направлении, противоположном градиенту: \[ \vec{x}_{k+1} = \vec{x}_k - \alpha_k \nabla f(\vec{x}_k) \] где $\alpha_k$ --- шаг обучения (learning rate). \textbf{Требования:} Функция должна быть дифференцируемой (гладкой). \textbf{Плюсы:} \begin{itemize} \item Прост в реализации, \item Гарантированная сходимость к локальному минимуму (при правильном выборе $\alpha$), \item Хорошо работает в высоких размерностях. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Застревает в локальных минимумах, \item Может очень медленно сходиться в ``оврагах'' (когда собственные значения гессиана сильно различаются), \item Требует выбора шага $\alpha$ (слишком большой --- расходится, слишком маленький --- медленно). \end{itemize} \subsubsection{Метод сопряженных направлений (Conjugate Gradient для оптимизации)} \textbf{Идея:} Строить направления поиска так, чтобы они были ``сопряжены'' относительно гессиана функции. Это позволяет найти минимум квадратичной функции за $N$ шагов (где $N$ --- размерность). \textbf{Важное замечание:} Метод сопряженных градиентов для решения линейных систем итерационным методом \textit{только называется} одинаково с методом сопряженных градиентов для нелинейной минимизации. На самом деле, это один и тот же алгоритм, просто примененный к разным задачам: \begin{itemize} \item Для линейных систем $\mathbf{A}\vec{x} = \vec{b}$ --- это минимизация квадратичной функции $f(\vec{x}) = \frac{1}{2}\vec{x}^T\mathbf{A}\vec{x} - \vec{b}^T\vec{x}$, \item Для нелинейной оптимизации --- это обобщение той же идеи на произвольные функции. \end{itemize} \textbf{Требования:} Гладкая функция, желательно с непрерывным градиентом. \textbf{Плюсы:} \begin{itemize} \item Сходится быстрее обычного градиентного спуска, \item Не требует хранения гессиана (в отличие от метода Ньютона), \item Память --- $\mathcal{O}(N)$. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Застревает в локальных минимумах, \item Для неквадратичных функций требуется ``перезапуск'' (reset) направлений. \end{itemize} \subsubsection{Методы Ньютона} \textbf{Идея:} Использовать не только градиент, но и вторые производные (гессиан). Разлагаем функцию в ряд Тейлора до второго порядка: \[ f(\vec{x} + \vec{h}) \approx f(\vec{x}) + \nabla f(\vec{x})^T \vec{h} + \frac{1}{2} \vec{h}^T \mathbf{H}(\vec{x}) \vec{h} \] и находим минимум этой квадратичной аппроксимации. Получаем итерацию: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{H}^{-1}(\vec{x}_k) \nabla f(\vec{x}_k) \] где $\mathbf{H}$ --- матрица вторых производных (гессиан). \textbf{Требования:} Функция должна быть дважды дифференцируемой. Гессиан должен быть вычислим. \textbf{Плюсы:} \begin{itemize} \item \textbf{Квадратичная сходимость:} если начали близко к решению, то число верных цифр удваивается на каждой итерации, \item Не чувствителен к ``оврагам'' (в отличие от градиентного спуска). \end{itemize} \textbf{Минусы:} \begin{itemize} \item Требует вычисления гессиана --- $\mathcal{O}(N^2)$ элементов, \item Требует решения линейной системы с гессианом --- $\mathcal{O}(N^3)$ операций, \item Может расходиться, если начали далеко от решения, \item Гессиан может быть вырожденным или не положительно определенным. \end{itemize} \subsubsection{Метод Ньютона--Рафсона и блочный Ньютон для квантовой механики} \textbf{Метод Ньютона--Рафсона} --- это вариант метода Ньютона для решения систем нелинейных уравнений $\vec{F}(\vec{x}) = \vec{0}$: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{J}^{-1}(\vec{x}_k) \vec{F}(\vec{x}_k) \] где $\mathbf{J}$ --- матрица Якоби (матрица первых производных $\partial F_i / \partial x_j$). \textbf{Блочный Ньютон для квантовой механики:} В задачах квантовой химии (например, в методе Хартри--Фока или DFT) часто возникает система уравнений, которую можно разбить на блоки: \begin{itemize} \item Уравнения для орбиталей (коэффициенты разложения), \item Уравнения для плотности, \item Уравнения для энергии. \end{itemize} Блочный метод Ньютона учитывает эту структуру: гессиан представляется в блочном виде, и каждый блок обрабатывается отдельно. Это позволяет: \begin{itemize} \item Эффективно использовать структуру задачи, \item Параллелить вычисления, \item Применять разные методы к разным блокам (например, Ньютон для одного блока, градиентный спуск для другого). \end{itemize} \subsubsection{Методы BFGS и L-BFGS} \textbf{BFGS (Broyden--Fletcher--Goldfarb--Shanno)} --- это квазиньютоновский метод. Идея: не вычислять гессиан явно, а приближать его на каждой итерации, используя информацию о изменении градиента. \textbf{Требования:} Гладкая функция с непрерывным градиентом. \textbf{Плюсы:} \begin{itemize} \item Сходимость почти как у Ньютона, но без вычисления гессиана, \item Гарантированная положительная определенность приближения гессиана, \item Хорошо работает на практике --- один из самых популярных методов. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Требует хранения матрицы $N \times N$ (приближение гессиана) --- $\mathcal{O}(N^2)$ памяти. \end{itemize} \textbf{L-BFGS (Limited-memory BFGS)} --- модификация BFGS для больших задач. Идея: не хранить полную матрицу, а хранить только последние $m$ пар векторов $(\vec{s}_k, \vec{y}_k)$, где: \[ \vec{s}_k = \vec{x}_{k+1} - \vec{x}_k, \quad \vec{y}_k = \nabla f_{k+1} - \nabla f_k \] \textbf{Плюсы:} \begin{itemize} \item Память --- $\mathcal{O}(mN)$, где $m \sim 5\text{--}20$ (не зависит от $N^2$!), \item Идеален для больших задач ($N > 1000$), \item Очень популярен в машинном обучении и квантовой химии. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Чуть медленнее сходится, чем полный BFGS, \item Требует настройки параметра $m$. \end{itemize} \subsubsection{Симплекс-методы (Nelder--Mead)} \textbf{Идея:} Строим симплекс (в $N$-мерном пространстве это $N+1$ точка), вычисляем функцию в вершинах, и итеративно отражаем, сжимаем или расширяем симплекс, двигаясь к минимуму. \textbf{Требования:} Никаких! Функция может быть негладкой, шумной, даже разрывной. \textbf{Плюсы:} \begin{itemize} \item Не требует градиентов, \item Прост в реализации, \item Хорошо работает для малых размерностей ($N < 10$). \end{itemize} \textbf{Минусы:} \begin{itemize} \item Очень медленный для больших $N$, \item Может застревать, \item Нет теоретических гарантий сходимости. \end{itemize} \textbf{Когда использовать:} Когда функция очень шумная или когда $N$ маленькое и лень выводить градиенты. \subsubsection{Simulated Annealing (Имитация отжига)} \textbf{Идея:} Вдохновлен физическим процессом отжига металлов. Начинаем с высокой ``температуры'' $T$, которая позволяет алгоритму ``перепрыгивать'' через локальные минимумы. Постепенно снижаем $T$, и алгоритм ``застывает'' в глобальном минимуме. На каждом шаге: \begin{enumerate} \item Предлагаем случайное изменение $\vec{x} \to \vec{x}'$, \item Вычисляем $\Delta f = f(\vec{x}') - f(\vec{x})$, \item Если $\Delta f < 0$ --- принимаем изменение, \item Если $\Delta f > 0$ --- принимаем с вероятностью $P = \exp(-\Delta f / T)$. \end{enumerate} \textbf{Требования:} Никаких. \textbf{Плюсы:} \begin{itemize} \item Может найти глобальный минимум (при правильном ``расписании отжига''), \item Не застревает в локальных минимумах на ранних этапах, \item Прост в реализации. \end{itemize} \textbf{Минусы:} \begin{itemize} \item Очень медленный, \item Требует настройки расписания $T(t)$, \item Нет гарантий сходимости за разумное время. \end{itemize} \textbf{Когда использовать:} Когда задача мультимодальная (много локальных минимумов) и нужно найти глобальный. \subsection{Вычисление градиентов} Все градиентные методы требуют вычисления $\nabla f(\vec{x})$. Как это делать? \subsubsection{Конечные разности} \textbf{Идея:} Приближаем производную: \[ \frac{\partial f}{\partial x_j} \approx \frac{f(\vec{x} + h \vec{e}_j) - f(\vec{x})}{h} \] где $\vec{e}_j$ --- $j$-й базисный вектор. \textbf{Проблемы:} \begin{enumerate} \item \textbf{Дорого:} Для вычисления градиента в $N$-мерном пространстве нужно $N+1$ вычислений функции (или $2N$ для центральной разности). \item \textbf{Неустойчиво:} Если $h$ слишком большое --- большая ошибка аппроксимации. Если $h$ слишком маленькое --- потеря точности из-за вычитания близких чисел. Оптимальное $h \sim \sqrt{\epsilon_{\text{mach}}}$, где $\epsilon_{\text{mach}}$ --- машинная точность. \item \textbf{Шум:} Если функция вычисляется с шумом (например, экспериментальные данные), то конечные разности усиливают шум. \end{enumerate} \subsubsection{Автоматическое дифференцирование (AD): Метод Бауэра--Штрассена: аналитическое вычисление градиента} А теперь --- магия! Оказывается, можно вычислить градиент функции \textit{аналитически}, используя автоматическое дифференцирование, и сделать это всего в \textbf{несколько раз} дороже, чем вычисление самой функции! \textbf{Идея:} Любая функция $f(\vec{x})$ --- это композиция элементарных операций ($+$, $-$, $\times$, $\div$, $\sin$, $\cos$, $\exp$, $\log$, и т.д.). Мы можем: \begin{enumerate} \item Представить вычисление функции как \textbf{вычислительный граф}, \item Применить \textbf{цепное правило} (chain rule) в обратном порядке (reverse mode). \end{enumerate} \textbf{Результат:} Градиент вычисляется за $\mathcal{O}(1) \times$ стоимость функции (на практике --- в 3--5 раз дороже, но не в $N$ раз!). \textbf{Как это работает на пальцах:} \begin{enumerate} \item Прямой проход: вычисляем $f(\vec{x})$, сохраняя все промежуточные результаты. \item Обратный проход: идем по графу в обратном направлении, применяя цепное правило: \[ \frac{\partial f}{\partial x_j} = \sum_{\text{пути}} \frac{\partial f}{\partial u_1} \frac{\partial u_1}{\partial u_2} \cdots \frac{\partial u_k}{\partial x_j} \] \end{enumerate} \textbf{Где это используется:} \begin{itemize} \item \textbf{PyTorch, TensorFlow:} Основа всех современных нейросетей --- это именно reverse-mode automatic differentiation. \item \textbf{Квантовая химия:} Программы типа PySCF, Psi4 используют AD для вычисления градиентов энергии по координатам ядер. \item \textbf{Оптимизация:} Все современные библиотеки (SciPy, JAX) используют AD. \end{itemize} \begin{successbox}[Главный вывод] Никогда не вычисляйте градиенты конечными разностями, если можно использовать автоматическое дифференцирование! Это быстрее, точнее и устойчивее. \end{successbox} \subsection{Живой пример: силовые поля и конформеры} Ещё один живой пример из химии --- моделирование молекул. Для моделирования небольших молекул часто используют такое представление, когда общая энергия записана как сумма простых взаимодействий: \begin{enumerate} \item \textbf{Ван-дер-Ваальсовы взаимодействия} между несвязанными атомами (потенциал Леннард-Джонса): \[ E_{\text{vdW}} = \sum_{i < j} 4\varepsilon_{ij} \left[ \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6 \right] \] \item \textbf{Гармонические связи} между связанными атомами: \[ E_{\text{bond}} = \sum_{\text{bonds}} \frac{1}{2} k_b (r - r_0)^2 \] \item \textbf{Угловые взаимодействия} между тремя последовательными атомами: \[ E_{\text{angle}} = \sum_{\text{angles}} \frac{1}{2} k_\theta (\theta - \theta_0)^2 \] \item \textbf{Торсионные (диэдральные) взаимодействия} между четырьмя последовательными атомами: \[ E_{\text{torsion}} = \sum_{\text{dihedrals}} \frac{V_n}{2} [1 + \cos(n\phi - \gamma)] \] \end{enumerate} Их будет не много, но около квадрата от числа атомов (для ван-дер-ваальса). Каждая такая функция --- это какая-то несложная формула: иногда с синусами, иногда с экспонентами, иногда со степенями. А суммарно у молекулы --- $3N$ степеней свободы, так как каждый атом имеет 3 пространственные координаты. \subsubsection{Проблема конформеров} Казалось бы, задача ясна: надо найти минимум этой функции $E(\vec{x})$, где $\vec{x} \in \R^{3N}$ --- координаты всех атомов. Но тут начинается самое интересное. Функция энергии молекулы --- это \textbf{мультимодальный ландшафт} с огромным числом локальных минимумов. Каждый локальный минимум соответствует своей \textbf{конформации} (конформеру) молекулы --- своему способу скручивания молекулы в пространстве. Например, для белка из 100 аминокислот число возможных конформеров может достигать $10^{100}$ --- это знаменитый \textbf{парадокс Левинталя}. И глобальный минимум (нативная структура) --- это лишь одна из $10^{100}$ возможностей. \subsubsection{Почему простая минимизация не работает} Если мы возьмём случайную начальную конформацию и запустим обычный градиентный спуск или BFGS, то мы очень быстро (за несколько сотен итераций) сойдёмся в \textit{ближайший} локальный минимум. Но это почти наверняка \textit{не} будет глобальный минимум --- то есть не та конформация, которую молекула принимает в реальности. \subsubsection{Решение: комбинация методов} Поэтому на практике используют комбинацию методов: \begin{enumerate} \item \textbf{Simulated Annealing (имитация отжига):} Начинаем с высокой ``температуры'', которая позволяет алгоритму перепрыгивать через энергетические барьеры и исследовать разные области конформационного пространства. Постепенно снижаем температуру, ``замораживая'' систему в низкоэнергетической области. \item \textbf{Локальная минимизация (BFGS):} После того как simulated annealing нашёл ``хорошую'' область, запускаем быстрый локальный метод (BFGS или метод сопряжённых градиентов) для точной сходимости к локальному минимуму. \item \textbf{Повторение:} Запускаем эту комбинацию много раз с разными начальными условиями и выбираем конформер с наименьшей энергией. \end{enumerate} \begin{tipbox}[Популярные программы] Этот подход реализован в популярных программах молекулярной динамики: \begin{itemize} \item \textbf{AMBER, GROMACS, NAMD:} Используют силовые поля (AMBER, CHARMM, OPLS) и методы типа simulated annealing + локальная минимизация. \item \textbf{Rosetta:} Для предсказания структуры белков использует фрагментный подход + Monte Carlo + локальная минимизация. \end{itemize} \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Быстрое преобразование Фурье: когда $O(N^2)$ превращается в $O(N \log N)$} \subsection{Задача поиска паттерна} Иногда нам надо найти какой-то определенный фрагмент в огромной последовательности. Если это --- очень короткий фрагмент, то простой перебор всей последовательности и сравнение с этим коротким фрагментом --- это довольно хорошая стратегия. Так часто поступают при поиске коротких фрагментов в белковых структурах. Но иногда размерности того, что надо сравнить, почти совпадают. Например, у нас есть какая-то периодическая кристаллическая структура, и мы знаем, что в ней есть какое-то отклонение, и даже понимаем, что отклонение только частично похоже на то, с чем мы хотим сравнить. Формально: исходные данные --- это $v_i$, $i = 1, \dots, N$, и то, с чем надо сравнить --- $w_j$, $j = 1, \dots, M$, где $M < N$. Мы можем явно в цикле посчитать: \[ \forall k = 0, \dots, N-M: \quad s_k = \sum_{j=1}^M w_j v_{j+k} \] и найти максимум по $s$. Но вычислительная сложность такого поиска будет квадратична по $N$: если $M \simeq N/2$, то надо выполнить $N^2/2$ арифметических операций. Для $N = 10^6$ это $5 \times 10^{11}$ операций --- слишком много! \subsection{Матричная формулировка} Тут люди придумали, как найти такие $s_k$ гораздо быстрее. Давайте представим вычисление этих $s$ в матричном виде. Пусть у нас есть матрица: \[ \mathbf{W} = \begin{pmatrix} w_1 & w_2 & \cdots & w_M & 0 & \cdots & 0 \\ 0 & w_1 & w_2 & \cdots & w_M & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & & \ddots & 0 \\ 0 & \cdots & 0 & w_1 & w_2 & \cdots & w_M \end{pmatrix} \] размера $(N-M+1) \times N$. Если мы умножим её на наш вектор $\vec{v}$, то получим как раз наши заветные $s_k$: \[ \vec{s} = \mathbf{W} \vec{v} \] Но как же нам это сделать быстро? \subsection{Циркулянтная матрица} Если мы дополним эту матрицу снизу ещё $M-1$ строками вида: \begin{align*} & w_M, 0, \dots, 0, w_1, \dots, w_{M-1} \\ & w_{M-1}, w_M, 0, \dots, 0, w_1, \dots, w_{M-2} \\ & \vdots \\ & w_2, \dots, w_M, 0, \dots, 0, w_1 \end{align*} то после умножения мы получим тот же вектор $\vec{s}$, но в конце к нему будет приписано ещё $M-1$ каких-то чисел, которые мы можем отбросить. Но эта расширенная матрица --- это \textbf{циркулянтная матрица} $\mathbf{C}$! Она обладает замечательным свойством: каждая следующая строка получается циклическим сдвигом предыдущей. \subsection{Диагонализация циркулянтной матрицы} У циркулянтной матрицы есть интересная запись через матрицу Фурье: \[ \mathbf{C} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \] где: \begin{itemize} \item $\vec{c}$ --- первый столбец исходной матрицы $\mathbf{C}$, \item $\mathbf{F}$ --- матрица дискретного преобразования Фурье (ДПФ), \item $\mathbf{F}^H$ --- эрмитово сопряжённая (комплексно сопряжённая и транспонированная) матрица. \end{itemize} Матрица Фурье определяется как: \[ \mathbf{F} = \{f_{jk}\}_{j,k=0}^{N-1}, \quad f_{jk} = e^{-2\pi i j k / N} \] \subsection{Быстрое умножение через FFT} Теперь самое интересное. Если мы умеем быстро умножать матрицу Фурье на вектор, то мы сможем быстрее вычислить этот $\vec{s}$: \[ \vec{s} = \mathbf{C} \vec{v} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \vec{v} \] Это вычисление состоит из трёх шагов: \begin{enumerate} \item $\vec{a} = \mathbf{F} \vec{v}$ --- прямое ДПФ вектора $\vec{v}$, \item $\vec{b} = \diag(\mathbf{F} \vec{c}) \vec{a}$ --- поэлементное умножение (это $\mathcal{O}(N)$), \item $\vec{s} = \frac{1}{N} \mathbf{F}^H \vec{b}$ --- обратное ДПФ. \end{enumerate} \subsection{Разложение матрицы Фурье} Как же быстро умножить $\mathbf{F}$ на вектор? Оказывается, $\mathbf{F}$ можно представить в виде произведения специальных матриц: \begin{enumerate} \item Одной матрицы перестановки (bit-reversal permutation), \item Нескольких пар матриц: \begin{itemize} \item Блочных матриц вида $\begin{pmatrix} \mathbf{I} & \mathbf{I} \\ \mathbf{I} & -\mathbf{I} \end{pmatrix}$, \item Диагональных матриц с элементами $e^{-2\pi i k / N}$ (так называемые ``twiddle factors''). \end{itemize} \end{enumerate} Давайте наберём эти матрицы явно для $N = 8$: \textbf{Матрица перестановки (bit-reversal):} \[ \mathbf{P} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \] (строки переставлены согласно обращению бит: $0, 4, 2, 6, 1, 5, 3, 7$) \textbf{Блочная матрица (первый уровень):} \[ \mathbf{B}_1 = \begin{pmatrix} 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & -1 \end{pmatrix} \] \textbf{Диагональная матрица (twiddle factors, первый уровень):} \[ \mathbf{D}_1 = \diag\left(1, 1, 1, 1, 1, e^{-2\pi i / 8}, e^{-4\pi i / 8}, e^{-6\pi i / 8}\right) \] И так далее для $\log_2 N$ уровней. Тогда каждое такое умножение будет требовать только $N$ и $2N$ арифметических операций, а всего таких шагов будет $\log_2 N$. \textbf{Итого:} умножить матрицу Фурье на вектор можно за $\mathcal{O}(N \log_2 N)$ арифметических операций! А значит, и посчитать все $s_k$ можно за те же $\mathcal{O}(N \log_2 N)$, а не за $\mathcal{O}(N^2)$. Для $N = 10^6$: \begin{itemize} \item Наивный алгоритм: $10^{12}$ операций, \item FFT: $2 \times 10^7$ операций (в 50\,000 раз быстрее!). \end{itemize} \subsection{Быстрое преобразование Фурье (FFT)} Это и есть знаменитый \textbf{метод Быстрого Преобразования Фурье} (Fast Fourier Transform, FFT), открытый Кули и Тьюки в 1965 году (хотя Гаусс знал его ещё в 1805-м!). Он является одним из самых востребованных алгоритмов в современной химии. \subsubsection{Применение в ЯМР} С его помощью, например, преобразуют исходные FID (Free Induction Decay) из ЯМР в спектральную область. Если посмотреть на каждую строку матрицы Фурье, то в ней можно увидеть осциллирующую функцию: \begin{itemize} \item Чем ближе строка к центру матрицы, тем больше осцилляция, \item У первой строки её вообще нет --- это просто константа, \item У самой нижней (или второй) строки период этой осцилляции равен длине строки. \end{itemize} Именно поэтому, умножив матрицу Фурье на FID, мы получаем максимумы там, где есть резонансные частоты: просто именно такая строчка совпала по частоте осцилляции с резонансной частотой. \begin{tipbox}[Почему FFT так важен для ЯМР] Современные FID очень длинные --- их часто оцифровывают по сотни тысяч и даже миллионы чисел. Если бы быстрого преобразования Фурье не было, умножение на такую матрицу длилось бы на современных компьютерах даже дольше, чем сам съём ЯМР-спектров --- то есть реально минуты и десятки минут даже на современных рабочих станциях. Благодаря FFT преобразование занимает \textbf{доли секунды}. \end{tipbox} \subsubsection{Другие применения FFT в химии} \begin{itemize} \item \textbf{Крио-электронная микроскопия:} Реконструкция 3D-структур из 2D-проекций. \item \textbf{Рентгеноструктурный анализ:} Преобразование дифракционной картины в электронную плотность. \item \textbf{Молекулярная динамика:} Вычисление электростатических взаимодействий через метод PPPM (Particle-Particle Particle-Mesh). \item \textbf{Обработка сигналов:} Фильтрация шума, сжатие данных. \item \textbf{Хемометрика:} Быстрая свёртка и корреляция спектров. \end{itemize} \begin{successbox}[Итог] FFT --- это один из тех алгоритмов, которые изменили мир. Без него не было бы современной спектроскопии, обработки сигналов, сжатия аудио и видео, и многого другого. Это must-know для любого химика, работающего с данными. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Когда БПФ не видит очевидного: утечка спектра и метод Прони} \subsection{Парадокс: глаз видит синусоиду, а БПФ --- нет} Хотя БПФ --- это must-have алгоритм практически везде в экспериментальной химии, у него тоже есть свои особенности, на которые надо однозначно обращать внимание. Рассмотрим пример. Мы сняли спектр, сигнал осциллирует, и это прямо видно глазами на графике. Но мы смогли записать только слегка больше одного периода. Когда мы применили БПФ, получили что-то очень странное: пик вроде бы есть, но и вроде бы его нет --- он получился каким-то очень размазанным. А если у нас было много другого шума, то мы получили в другой части спектра что-то очень похожее, и да, случайно шум оказался ``ярче'' сигнала, и наша программа выдала совершенно неправильный результат. Но мы же глазами видим эту синусоиду! Почему же БПФ её не видит? \subsection{Утечка спектра (Spectral Leakage)} БПФ хорош тогда, и \textit{только тогда}, когда мы применили его к сигналу, в котором укладывается \textbf{точно целое число периодов} этого сигнала. Почему? Потому что БПФ неявно предполагает, что наш конечный сигнал \textit{периодически продолжается} за пределы окна измерения. Если в окне укладывается ровно $k$ полных периодов, то периодическое продолжение будет гладким, и БПФ даст один чёткий пик. Но если сигнал содержит, скажем, $1.4$ периода, то при периодическом продолжении возникнет \textbf{разрыв} на границе окна. Этот разрыв БПФ вынужден аппроксимировать множеством гармоник --- и энергия ``утекает'' из основного пика во все остальные частоты. Это явление называется \textbf{утечкой спектра} (spectral leakage). Математически: если истинная частота сигнала попадает \textit{между} двумя соседними частотными бинами БПФ, то вместо одной дельта-функции мы получаем размазанный ``горб'' (функцию вида $\sin(x)/x$, так называемый \textit{синк}). \begin{warningbox}[Правило БПФ] БПФ видит только те частоты, которые укладываются в окно измерения целое число раз. Всё остальное --- утечка. \end{warningbox} \subsection{Zero-padding: простое, но не идеальное решение} Что же нам тогда делать? Самый простой вариант --- добавить много нулей в конец сигнала (zero-padding) и надеяться, что суммарно нам существенно больше повезёт. Ведь если мы, к примеру, добавим нули так, что исходные данные растянутся в 5 раз, то у нас точно получится 7 периодов (вместо 1.4). А если увеличим, скажем, в 7 раз, то будет 9.8 периода --- тоже очень близко к целому. Люди часто добавляют разное число нулей, а потом пытаются объединить получаемые спектры, угадывая, какие пики оказались более точными на соответствующих расширенных спектрах. Это реально помогает! \begin{tipbox}[Важный нюанс zero-padding] Zero-padding \textit{не увеличивает реальное частотное разрешение} --- он только интерполирует спектр, делая его более гладким. Разрешение определяется \textit{длиной исходного сигнала}, а не длиной дополненного нулями. Но на практике zero-padding помогает ``попасть'' в целое число периодов и уменьшить утечку. \end{tipbox} \subsection{Дрейф частоты: когда сигнал ``уплывает''} Но есть ещё одна проблема, с которой zero-padding не справится. Иногда мы измеряем долго и нудно какой-то спектр, а он ``слегка'' взял и поплыл. То есть в начале у нас была частота $\omega$, а под конец $\omega + \varepsilon$, и этой маленькой поправки достаточно, чтобы мы совершенно не попали по БПФ и получили вместо чёткого единичного пика что-то очень размытое. Но ведь глазом-то видно --- вот она синусоида, и тут, и там, и в начале, и в конце! Только вот БПФ её опять не видит. Причина в том, что БПФ предполагает \textbf{стационарность} сигнала --- то есть что частоты не меняются со временем. Если частота дрейфует, то ни одна гармоника БПФ не может хорошо аппроксимировать весь сигнал целиком. Что же нам делать? \subsection{Метод Прони: идея сдвигов} На этот вопрос нам несколько веков назад ответил Гаспар Прони (Gaspard de Prony, 1795). Вернее, он придумал, как решать свою задачу (разложение сигнала на сумму экспонент), но она как раз решает и нашу проблему --- когда спектр во времени слегка уплывает. \textbf{Идея:} Возьмём сигнал и сохраним его в вектор $\vec{s} = (s_0, s_1, \dots, s_{N-1})^T$. Далее сдвинем этот вектор вниз на один отсчёт, потом ещё раз, и снова поставим рядом. Сделаем так $L$ раз. Что же мы будем тут иметь? Если у нас был периодический сигнал с какой-то неизвестной синусоидой $\sin(\omega t)$, то после сдвига на $p$ отсчётов мы получаем $\sin(\omega t + \omega p)$. Вспоминая формулы сложения из главы 1: \[ \sin(\omega t + \omega p) = \sin(\omega t)\cos(\omega p) + \cos(\omega t)\sin(\omega p) \] Поскольку $\cos(\omega p)$ и $\sin(\omega p)$ --- это \textit{константы} для всего вектора (они не зависят от $t$), то набор таких сдвинутых векторов будет иметь только \textbf{столько ненулевых сингулярных чисел, сколько у нас удвоенное число различных осциллирующих гармоник}. Более того, если какая-то частота со временем немного ``ушла'', это \textit{никак не повлияет} на такую аппроксимацию, но совершенно разрушит результат БПФ, как мы отмечали выше. \subsection{Метод Прони, вариант 1: линейное предсказание} Пусть наш сигнал моделируется как сумма $K$ комплексных экспонент: \[ s_n = \sum_{m=1}^{K} c_m z_m^n, \quad n = 0, 1, \dots, N-1 \] где $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ --- комплексные ``частоты'' (содержащие и затухание $\alpha_m$, и осцилляцию $\omega_m$), а $c_m$ --- комплексные амплитуды. \textbf{Ключевой факт:} Если сигнал состоит из $K$ экспонент и нет шума, то он удовлетворяет \textbf{линейному рекуррентному соотношению} порядка $K$: \[ s_{n} + a_1 s_{n-1} + a_2 s_{n-2} + \cdots + a_K s_{n-K} = 0 \quad \forall n \ge K \] Это означает, что $(n+1)$-й отсчёт можно \textit{точно предсказать} через $K$ предыдущих! Коэффициенты $a_1, \dots, a_K$ --- это коэффициенты характеристического полинома, корнями которого являются $z_m$: \[ P(z) = z^K + a_1 z^{K-1} + a_2 z^{K-2} + \cdots + a_K = \prod_{m=1}^{K}(z - z_m) \] \textbf{Как найти коэффициенты?} Запишем рекуррентное соотношение для всех доступных отсчётов в матричном виде: \[ \underbrace{ \begin{pmatrix} s_{K-1} & s_{K-2} & \cdots & s_0 \\ s_{K} & s_{K-1} & \cdots & s_1 \\ \vdots & \vdots & \ddots & \vdots \\ s_{N-2} & s_{N-3} & \cdots & s_{N-K-1} \end{pmatrix} }_{\mathbf{S} \;\; (N-K) \times K} \underbrace{ \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_K \end{pmatrix} }_{\vec{a}} = - \underbrace{ \begin{pmatrix} s_K \\ s_{K+1} \\ \vdots \\ s_{N-1} \end{pmatrix} }_{\vec{s}_{\text{future}}} \] Это переопределённая система (если $N > 2K$), и мы решаем её методом наименьших квадратов: \[ \vec{a} = -(\mathbf{S}^H \mathbf{S})^{-1} \mathbf{S}^H \vec{s}_{\text{future}} \] После нахождения $\vec{a}$ мы ищем корни полинома $P(z)$ --- это и есть наши $z_m$, из которых извлекаются частоты $\omega_m$ и затухания $\alpha_m$. \subsection{Метод Прони, вариант 2: SVD матрицы сдвигов} Второй вариант --- более устойчивый к шуму и более элегантный. \textbf{Шаг 1: Строим матрицу Ганкеля из сдвигов сигнала.} \[ \mathbf{H} = \begin{pmatrix} s_0 & s_1 & \cdots & s_{L-1} \\ s_1 & s_2 & \cdots & s_L \\ \vdots & \vdots & \ddots & \vdots \\ s_{N-L} & s_{N-L+1} & \cdots & s_{N-1} \end{pmatrix} \] размера $(N-L+1) \times L$, где $L > K$ (мы выбираем $L$ заведомо больше ожидаемого числа гармоник). \textbf{Шаг 2: Делаем SVD.} \[ \mathbf{H} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] Если сигнал состоит из $K$ экспонент и нет шума, то $\rank(\mathbf{H}) = K$, и последние $L - K$ сингулярных чисел равны нулю: \[ \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_K > 0, \quad \sigma_{K+1} = \cdots = \sigma_L = 0 \] В присутствии шума малые сингулярные числа не равны нулю, но они \textit{существенно меньше} первых $K$. Мы отбрасываем их (это называется \textbf{усечённое SVD} или \textbf{регуляризация}). \textbf{Шаг 3: Извлекаем полином из нуль-пространства.} Сингулярный вектор $\vec{v}_{\min}$, соответствующий \textit{самому маленькому} сингулярному числу (последний столбец $\mathbf{V}$), лежит в (приближённом) нуль-пространстве матрицы $\mathbf{H}$. Это означает: \[ \mathbf{H} \vec{v}_{\min} \approx \vec{0} \] Если записать компоненты $\vec{v}_{\min} = (v_0, v_1, \dots, v_{L-1})^T$, то это соотношение эквивалентно тому, что: \[ \sum_{j=0}^{L-1} v_j s_{n+j} \approx 0 \quad \forall n \] То есть $\vec{v}_{\min}$ содержит коэффициенты полинома: \[ Q(z) = v_0 + v_1 z + v_2 z^2 + \cdots + v_{L-1} z^{L-1} \] \textbf{Шаг 4: Ищем корни полинома.} Корни $Q(z)$ содержат $K$ ``сигнальных'' корней $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ (лежащих внутри или на единичной окружности) и $L - K$ ``шумовых'' корней (разбросанных случайно). Из сигнальных корней извлекаем: \[ \omega_m = \frac{\arg(z_m)}{\Delta t}, \quad \alpha_m = \frac{\ln|z_m|}{\Delta t} \] \begin{tipbox}[Почему SVD-вариант лучше?] SVD-вариант метода Прони устойчивее к шуму, потому что: \begin{enumerate} \item Усечение малых сингулярных чисел автоматически фильтрует шум. \item Мы не решаем переопределённую систему напрямую (что может быть плохо обусловлено), а используем ортогональное разложение. \item Число гармоник $K$ определяется автоматически по ``скачку'' в спектре сингулярных чисел. \end{enumerate} \end{tipbox} \subsection{Ограничения метода Прони} Но Прони --- это не панацея. Стоит попасться в исходном сигнале составляющей, которая хорошо аппроксимируется как самой функцией, так и её сдвигами (например, белый шум или очень широкополосный сигнал), как мы сразу получаем на неё ``корень'', и потом оказывается, что он --- \textbf{ложный}. Другие ограничения: \begin{itemize} \item Метод предполагает, что сигнал --- это сумма \textit{экспонент} (включая синусоиды как частный случай). Если сигнал имеет другую структуру (например, прямоугольные импульсы), метод будет работать плохо. \item Число гармоник $K$ нужно знать или угадывать заранее (хотя SVD помогает с этим). \item Поиск корней полинома высокой степени ($L > 50$) сам по себе может быть численно неустойчивым. \item Метод чувствителен к сильному шуму: при низком отношении сигнал/шум ложные корни могут ``замаскироваться'' под настоящие. \end{itemize} В то же время этот метод очень активно и довольно издавна применяется в ЯМР-спектроскопии и хорошо дополняет БПФ. На практике часто используют \textbf{гибридный подход}: \begin{enumerate} \item БПФ для быстрого обзора спектра и определения числа пиков, \item Метод Прони (или его современные варианты: MUSIC, ESPRIT, Matrix Pencil) для точного определения частот, затуханий и амплитуд отдельных пиков. \end{enumerate} \begin{successbox}[Итог] БПФ --- это мощный, но ``слепой'' инструмент: он видит только то, что укладывается в его частотную сетку. Метод Прони ``видит'' частоты между бинами и устойчив к дрейфу, но требует больше вычислений и осторожности. В реальной ЯМР-спектроскопии они работают в паре. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Наименьшие квадраты, невязки и регуляризация по Тихонову} \subsection{Постановка задачи} Итак, мы часто решаем задачи, которые мы называем задачами наименьших квадратов. Давай рассмотрим какие-нибудь из них, чтобы поточнее понять, как же это можно решать. \subsection{Пример из медицины: компьютерная томография (КТ)} Возьмём пример из смежной области --- из медицины. Пусть у нас есть рентгеновский КТ-сканер. Мы поместили пациента (или то, что мы исследуем) в неподвижном виде на какую-то площадку, и вокруг него с противоположных сторон располагаем передатчик и приёмники рентгеновских лучей. Мы выполняем измерения, насколько сильно излучение поглотилось исследуемым телом, и располагаем эти приёмники и передатчики в различных ракурсах. Обычно пациент лежит на кушетке, а приёмник и передатчики расположены на поверхности кольца. Это кольцо расположено так, что ось находится примерно вдоль кушетки, и само кольцо как вращается вдоль своей оси, так и передвигается вдоль кушетки. Что же тут происходит? Мягкие ткани почти не поглощают рентгеновское излучение, в то время как кости --- поглощают довольно сильно. Мы можем в уме разбить всё пространство, где расположен пациент, на небольшие (обычно одинаковые) кубики --- \textbf{воксели} (volume elements). Если мы знаем точно расположение передатчика и приёмника, то мы можем провести линию (рентгеновский луч), которая пересекает наши кубики. Пусть в каждом таком кубике у нас есть одно неизвестное --- это интенсивность поглощения рентгеновского излучения. Мы пронумеруем все эти кубики и запишем интенсивность их поглощения как неизвестный вектор $\vec{x} = (x_1, \dots, x_N)^T$. Тогда одно взаимное расположение передатчика и приёмника даст нам довольно разреженный набор коэффициентов --- длины прохождения луча через соответствующий кубик. Пусть это всё записано в строке матрицы $\mathbf{A}$, и число столбцов в ней равно $N$ (число вокселей), а число строк --- как раз равно числу проведённых измерений (пусть равное $M$). А вот сами измеренные значения (у нас может быть, например, один источник рентгена и много приёмников; тогда каждое взаимное расположение приёмника и источника --- это одна строка) будут записаны в векторе $\vec{b}$ тоже длины $M$. Тогда мы можем сформулировать задачу наименьших квадратов как: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2 \] то есть мы хотим, чтобы каждая строка матрицы $\mathbf{A}$, умноженная на $\vec{x}$ (суммарная интенсивность поглощения), была равна тому, что мы измеряем. \subsection{Нормальные уравнения} Мы помним, что такую задачу можно решить, но давай посмотрим на неё внимательно. Функция невязки: \[ E = \|\mathbf{A}\vec{x} - \vec{b}\|_2^2 = (\mathbf{A}\vec{x} - \vec{b})^H (\mathbf{A}\vec{x} - \vec{b}) = \vec{x}^H \mathbf{A}^H \mathbf{A} \vec{x} - 2 \vec{x}^H \mathbf{A}^H \vec{b} + \|\vec{b}\|_2^2 \] Если мы найдём производную $E$ по всем $x_j$, приравняем к нулю (в минимуме), то получим систему уравнений: \[ 2 \mathbf{A}^H \mathbf{A} \vec{x} - 2 \mathbf{A}^H \vec{b} = 0 \quad \Rightarrow \quad (\mathbf{A}^H \mathbf{A}) \vec{x} = \mathbf{A}^H \vec{b} \] Это так называемые \textbf{нормальные уравнения}. Тут вроде всё хорошо: у нас получается неотрицательно определённая матрица $\mathbf{A}^H \mathbf{A}$ и какая-то правая часть, с которой мы можем решить эту задачу. \subsection{Проблема вырожденности} А что будет, если мы построили такие кубики, но у нас не было ни одного эксперимента, в котором рентгеновский луч прошёл через, например, $i$-ый кубик? Очевидно, что тогда в исходной матрице $\mathbf{A}$ соответствующий $i$-ый столбец будет содержать только нули, а матрица $\mathbf{A}^H \mathbf{A}$ будет иметь нули как в $i$-ом столбце, так и в $i$-ой строке, то есть иметь гарантированно одно нулевое сингулярное значение. Пусть этот кубик был где-то в пространстве, и в нём нет части исследуемого тела, то есть --- на наше счастье --- нам этот кубик и не нужен. Но такая постановка \textit{испортит} решение: мы просто не сможем решить эту задачу, так как матрица $\mathbf{A}^H \mathbf{A}$ будет вырожденной. Более того, мы не можем гарантировать, что даже если в $\mathbf{A}^H \mathbf{A}$ нет таких нулевых столбцов и строк, то её обусловленность будет хорошей. То есть мы, вместо решения, можем получить совершенно неверные числа. \subsection{Решение через SVD и псевдообратную матрицу} Как же в этом случае быть? Вернёмся к сингулярному разложению матрицы $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$. Мы можем заметить, что: \begin{itemize} \item Если в $\mathbf{\Sigma}$ есть нули на диагонали, то они относятся к тем кубикам, через которые не проходил рентгеновский луч. \item Если есть что-то очень маленькое, то оно относится к таким экспериментам, когда рентгеновский луч только несколько раз и, возможно, только касаясь, попадал в какой-то кубик или набор кубиков. \end{itemize} То есть нам надо бы отказаться учитывать такие данные, но их очень сложно отфильтровать ещё на стадии формирования матрицы $\mathbf{A}$. Но раз эти данные малоинформативны, давай мы их просто ``выбросим''! То есть выполним это SVD и занулим маленькие диагональные значения в $\mathbf{\Sigma}$, и реконструируем $\mathbf{A}$. Это, конечно, хорошо, но мы всё равно не знаем, как решить такую задачу, так как $\mathbf{A}^H \mathbf{A}$ тоже будет содержать нулевые сингулярные значения. Но если для квадратной невырожденной системы $\mathbf{A}\vec{x} = \vec{b}$ можно представить решение как: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \quad \Rightarrow \quad \vec{x} = \mathbf{V} \mathbf{\Sigma}^{-1} \mathbf{U}^H \vec{b} \] то и нашу задачу наименьших квадратов мы, наверное, могли бы решить по аналогии. Правда, там, где в $\mathbf{\Sigma}$ у нас будут нули, в $\mathbf{\Sigma}^{-1}$ мы их не будем инвертировать, а запишем туда тоже нуль. Тогда неправильно будет называть её обратной, и будем называть её \textbf{псевдообратной}, обозначая как $\mathbf{\Sigma}^+$ (или $\mathbf{A}^+$ для всей матрицы). \subsection{Проблема размера: когда SVD не помещается в память} Все бы ничего, только вот в реальных задачах матрица $\mathbf{A}$ бывает огромной. Часто её элементы вычисляются аналитически, и нет необходимости их хранить. А плотная матрица $\mathbf{U}$, у которой размерность равна $\mathbf{A}$, обычно не помещается в память. Например, для обычного КТ-скана используют кубики размером около 1 мм. За секунду происходит около 100--500 измерений на 100--500 приёмниках, и всё измерение длится около минуты. То есть $N$ и $M$ легко могут быть порядка 10--100 миллионов, и полная сингулярная матрица потребует 100 петабайт памяти, что неразумно много. \subsection{Регуляризация по Тихонову} Можно заметить, что если к вырожденной матрице $\mathbf{A}^H \mathbf{A}$ добавить единичную $\lambda \mathbf{I}$, так что $\lambda$ будет в диапазоне тех самых маленьких сингулярных чисел, которые мы занулили, то такое добавление можно записать так: \begin{align*} \mathbf{A}^H \mathbf{A} + \lambda \mathbf{I} &= \mathbf{V} \mathbf{\Sigma} \mathbf{U}^H \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H + \lambda \mathbf{V} \mathbf{V}^H \\ &= \mathbf{V} (\mathbf{\Sigma}^2 + \lambda \mathbf{I}) \mathbf{V}^H \end{align*} И мы видим, что мы добавляем эту лямбду к каждому сингулярному числу, ``превращая'' вырожденную или плохо обусловленную матрицу в хорошо обусловленную. Более того, при малых $\lambda$ (порядка малых сингулярных чисел) мы почти не изменяем большие сингулярные числа, а вот маленькие просто заменяем на $\lambda$. Тогда, решая систему с такой матрицей, мы просто не ``портим'' большие сингулярные числа и существенно уменьшаем отклик у малых сингулярных чисел. Итак, обусловленность матрицы стала существенно меньшей, а мы не ``портим'' матрицей решение. Зная, что матрица разреженная, мы можем применить какие-то итерационные методы (даже метод сопряжённых градиентов) и очень быстро сойтись. Фактически, такая $\lambda$ --- это так называемый \textbf{регуляризатор по Тихонову}, ведь мы можем вместо исходной задачи решать такую: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2^2 + \lambda \|\vec{x}\|_2^2 \] фактически потребовав, чтобы одновременно минимизировались как невязка, так и норма решения, а их взаимная пропорция регулируется как раз значением этой $\lambda$. Тут есть много красивой теории, в своё время выведенной в 70-х годах прошлого века Тихоновым и его последователями, но суть остаётся простой: такой дополнительный ``регуляризатор'' помогает решать \textbf{плохо поставленные} (ill-posed) задачи. \subsection{Итеративное уменьшение $\lambda$} Более того, часто вначале ставят $\lambda$ довольно большим, тогда сопряжённый градиент сходится очень быстро (обусловленность у матрицы будет очень маленькая), а потом уменьшают $\lambda$, каждый раз используя предыдущее решение как начальное значение. Это позволяет: \begin{itemize} \item Быстро получить ``черновое'' решение с большим $\lambda$, \item Постепенно уточнять его, уменьшая $\lambda$, \item Избежать проблем с обусловленностью на ранних этапах. \end{itemize} \begin{successbox}[Итог] Регуляризация по Тихонову --- это мощный инструмент для решения плохо обусловленных задач наименьших квадратов. Она добавляет штраф за норму решения, что стабилизирует решение и позволяет использовать итерационные методы даже для вырожденных систем. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Базисные функции: от конечных разностей до сплайнов} \subsection{Идея: аппроксимировать неизвестное через известное} До сих пор мы работали с дискретными векторами и матрицами. Но реальные задачи --- дифференциальные уравнения, интегралы, волновые функции --- живут в непрерывном пространстве. Как же нам свести бесконечномерную задачу к конечномерной? Ответ: \textbf{базисные функции}. Мы представляем неизвестную функцию $f(x)$ как линейную комбинацию известных базисных функций $\phi_k(x)$: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] и ищем коэффициенты $c_k$. Это превращает задачу поиска функции в задачу поиска вектора $\vec{c} \in \R^N$. \subsection{Метод Ритца (вариационный метод)} Один из самых общих подходов --- \textbf{метод Ритца}. Пусть у нас есть функционал $E[f]$, который мы хотим минимизировать (например, энергия в квантовой механике). Мы подставляем разложение: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] и получаем функцию от коэффициентов: \[ E(\vec{c}) = E\left[\sum_{k=1}^N c_k \phi_k\right] \] Теперь мы минимизируем $E(\vec{c})$ по $\vec{c}$ --- это уже обычная задача оптимизации в $\R^N$. \subsection{Конечные разности (Finite Difference Method, FDM)} Самый простой выбор базисных функций --- это \textbf{локальные константы} или \textbf{линейные функции} на равномерной сетке. \subsubsection{Пример из химии: диффузия} Пусть у нас есть уравнение диффузии концентрации вещества $C(x,t)$: \[ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} \] Мы дискретизируем пространство с шагом $h$: $x_j = jh$, и время с шагом $\Delta t$: $t_n = n\Delta t$. Вторую производную аппроксимируем через центральную разность: \[ \frac{\partial^2 C}{\partial x^2}\bigg|_{x_j} \approx \frac{C_{j+1} - 2C_j + C_{j-1}}{h^2} \] где $C_j^n \approx C(x_j, t_n)$. Тогда схема Эйлера во времени: \[ \frac{C_j^{n+1} - C_j^n}{\Delta t} = D \frac{C_{j+1}^n - 2C_j^n + C_{j-1}^n}{h^2} \] Это даёт явную рекуррентную формулу: \[ C_j^{n+1} = C_j^n + \frac{D\Delta t}{h^2} (C_{j+1}^n - 2C_j^n + C_{j-1}^n) \] \begin{warningbox}[Устойчивость] Явная схема устойчива только при $\frac{D\Delta t}{h^2} \le \frac{1}{2}$. Если шаг по времени слишком большой --- решение ``развалится''. \end{warningbox} \subsection{Конечные элементы (Finite Element Method, FEM)} Более сложный, но более гибкий подход --- \textbf{конечные элементы}. Мы разбиваем область на маленькие элементы (треугольники, тетраэдры) и на каждом элементе аппроксимируем функцию \textbf{локальными полиномами} (обычно линейными или квадратичными). Базисные функции --- это \textbf{кусочно-линейные ``шляпы''} (hat functions), каждая из которых равна 1 в одном узле и 0 во всех остальных. Преимущества FEM: \begin{itemize} \item Можно работать со сложными геометриями, \item Можно адаптировать сетку (сгущать там, где решение быстро меняется), \item Хорошо теоретически обоснован. \end{itemize} \subsection{Базисные функции в квантовой химии: гауссовы орбитали} А теперь --- самое интересное для химиков. В квантовой химии (метод Хартри--Фока, DFT) мы решаем уравнение Шрёдингера: \[ \hat{H} \Psi = E \Psi \] где $\hat{H}$ --- оператор Гамильтона, $\Psi$ --- волновая функция. Мы представляем молекулярные орбитали $\psi_i(\vec{r})$ как линейную комбинацию \textbf{атомных орбиталей} (LCAO --- Linear Combination of Atomic Orbitals): \[ \psi_i(\vec{r}) = \sum_{\mu=1}^K c_{\mu i} \phi_\mu(\vec{r}) \] где $\phi_\mu(\vec{r})$ --- базисные функции, центрированные на атомах. \subsubsection{Гауссовы функции (Gaussian Type Orbitals, GTO)} В большинстве современных программ (Gaussian, ORCA, Q-Chem) в качестве базисных функций используют \textbf{гауссовы орбитали}: \[ \phi_\mu(\vec{r}) = (x - X_A)^l (y - Y_A)^m (z - Z_A)^n \exp(-\alpha |\vec{r} - \vec{R}_A|^2) \] где: \begin{itemize} \item $\vec{R}_A = (X_A, Y_A, Z_A)$ --- координаты атома $A$, \item $l, m, n$ --- угловые квантовые числа ($s, p, d, f$-орбитали), \item $\alpha$ --- показатель экспоненты (контролирует ``размер'' орбитали). \end{itemize} Почему именно гауссовы функции? \begin{itemize} \item \textbf{Произведение гауссов --- снова гаусс:} Это критически важно для вычисления многоцентровых интегралов (электрон-электронное отталкивание), \item \textbf{Аналитические интегралы:} Все интегралы вычисляются аналитически, \item \textbf{Контракты:} Реальные атомные орбитали аппроксимируются суммой нескольких гауссов (контракты), что повышает точность. \end{itemize} \subsubsection{Метод Ритца в квантовой химии} Мы подставляем разложение LCAO в функционал энергии Хартри--Фока или DFT и минимизируем по коэффициентам $c_{\mu i}$. Это приводит к задаче на собственные значения: \[ \mathbf{F} \vec{c}_i = \varepsilon_i \mathbf{S} \vec{c}_i \] где: \begin{itemize} \item $\mathbf{F}$ --- матрица Фока (или матрица Кона--Шэма в DFT), \item $\mathbf{S}$ --- матрица перекрывания ($S_{\mu\nu} = \langle \phi_\mu | \phi_\nu \rangle$), \item $\varepsilon_i$ --- орбитальные энергии, \item $\vec{c}_i$ --- коэффициенты разложения $i$-й молекулярной орбитали. \end{itemize} Это --- обобщённая задача на собственные значения, и она решается итерационно (методом самосогласованного поля, SCF). \begin{tipbox}[Базисные наборы] Популярные базисные наборы в квантовой химии: \begin{itemize} \item \textbf{STO-3G:} Минимальный базис (3 гаусса на каждую атомную орбиталь), \item \textbf{6-31G*:} Двойной дзета-качество с поляризационными функциями, \item \textbf{cc-pVTZ:} Корреляционно-согласованный тройной дзета (высокая точность). \end{itemize} Чем больше базис --- тем точнее результат, но тем дороже вычисления ($\mathcal{O}(N^4)$ для Хартри--Фока, где $N$ --- число базисных функций). \end{tipbox} \subsection{Сплайны: гладкие кусочно-полиномиальные функции} А теперь --- про сплайны. Это --- базисные функции, которые широко используются в обработке сигналов, компьютерной графике и численных методах. \subsubsection{Что такое сплайн?} \textbf{Сплайн} --- это кусочно-полиномиальная функция, которая является \textbf{гладкой} (имеет непрерывные производные до некоторого порядка) в узлах сшивки. Самый популярный --- \textbf{кубический сплайн} (степень 3, гладкость $C^2$ --- непрерывны функция, первая и вторая производные). \subsubsection{B-сплайны (Basis Splines)} \textbf{B-сплайны} --- это специальный набор базисных сплайнов с \textbf{конечным носителем} (compact support). Каждый B-сплайн отличен от нуля только на небольшом интервале. Преимущества: \begin{itemize} \item \textbf{Локальность:} Изменение одного коэффициента влияет только на небольшую область, \item \textbf{Численная устойчивость:} Матрицы хорошо обусловлены, \item \textbf{Рекурсивное определение:} B-сплайны строятся через рекурсию (формула де Бура). \end{itemize} \subsubsection{Двумерные сплайны} B-сплайны легко обобщаются на двумерный случай через \textbf{тензорное произведение}: \[ B_{ij}(x,y) = B_i(x) \cdot B_j(y) \] Это используется в: \begin{itemize} \item Обработке изображений (сглаживание, интерполяция), \item Конечных элементах (высшие порядки), \item Компьютерной графике (поверхности). \end{itemize} \subsubsection{Сплайны Безье (Bézier Splines)} \textbf{Сплайны Безье} --- это параметрические кривые, определяемые \textbf{контрольными точками}. Кривая не проходит через контрольные точки, но ``притягивается'' к ним. Формула кривой Безье степени $n$: \[ \vec{B}(t) = \sum_{i=0}^n \binom{n}{i} (1-t)^{n-i} t^i \vec{P}_i, \quad t \in [0,1] \] где $\vec{P}_i$ --- контрольные точки, а $\binom{n}{i} (1-t)^{n-i} t^i$ --- \textbf{полиномы Бернштейна}. Применение: \begin{itemize} \item Компьютерная графика (векторная графика, шрифты TrueType), \item CAD-системы (AutoCAD, SolidWorks), \item Анимация и траектории. \end{itemize} \subsubsection{Связь с конечными элементами} Интересно, что B-сплайны с конечным носителем \textbf{плавно перетекают} в базисные конечные элементы. Если взять B-сплайны степени 0 --- это кусочно-постоянные функции (как в простейших конечных элементах). Степени 1 --- кусочно-линейные ``шляпы''. Степени 2 и выше --- более гладкие базисы. Это позволяет строить \textbf{изопараметрические конечные элементы} высших порядков, которые дают более точную аппроксимацию при меньшем числе узлов. \begin{successbox}[Итог] Базисные функции --- это мост между непрерывным миром дифференциальных уравнений и дискретным миром компьютеров. Конечные разности --- самые простые, конечные элементы --- самые гибкие, гауссовы орбитали --- специализированы для квантовой химии, а сплайны --- для гладкой интерполяции и графики. Понимание их свойств --- ключ к выбору правильного метода для вашей задачи. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Интегрирование: от аналитики до Монте-Карло} \subsection{Зачем нам интегрирование?} Интеграл --- это одна из самых фундаментальных операций в математике и физике. Мы постоянно вычисляем: \begin{itemize} \item Площади и объёмы, \item Математические ожидания и дисперсии, \item Нормировочные константы в квантовой механике, \item Многомерные интегралы в статистической физике, \item Свёртки в обработке сигналов. \end{itemize} И если в школе мы учились брать интегралы аналитически, то в реальной жизни --- особенно в химии и физике --- аналитика часто заканчивается очень быстро. \subsection{Аналитическое интегрирование: когда оно работает} Вспомним азы. Если у нас есть функция $f(x)$, заданная аналитически, и мы знаем её первообразную $F(x)$, то: \[ \int_a^b f(x)\,dx = F(b) - F(a) \] Например: \[ \int_0^1 x^2\,dx = \frac{x^3}{3}\bigg|_0^1 = \frac{1}{3} \] Красиво, точно, быстро. Но: \begin{itemize} \item Не все функции имеют элементарную первообразную (например, $e^{-x^2}$ --- интеграл вероятностей), \item Не все функции заданы аналитически --- часто у нас есть только набор точек (экспериментальные данные), \item В многомерных случаях аналитика почти всегда бессильна. \end{itemize} \subsection{Численное интегрирование в 1D: метод Симпсона} Что же делать, если аналитически взять интеграл нельзя? Ответ: аппроксимировать функцию чем-то простым и проинтегрировать это простое. Мы в прошлой главе научились строить сплайны третьей степени. А тут --- мы строим тоже полиномы, и интегрируем такими кусками. \subsubsection{Метод прямоугольников} Самый простой подход: разбиваем отрезок $[a,b]$ на $N$ равных частей шириной $h = (b-a)/N$, и на каждом отрезке аппроксимируем функцию константой (значением в середине): \[ \int_a^b f(x)\,dx \approx h \sum_{k=0}^{N-1} f\left(a + \left(k + \frac{1}{2}\right)h\right) \] Погрешность --- $\mathcal{O}(h^2)$. \subsubsection{Метод трапеций} Чуть лучше: аппроксимируем функцию линейным полиномом на каждом отрезке: \[ \int_a^b f(x)\,dx \approx \frac{h}{2} \left[ f(a) + 2\sum_{k=1}^{N-1} f(a+kh) + f(b) \right] \] Погрешность --- $\mathcal{O}(h^2)$. \subsubsection{Метод Симпсона} Ещё лучше: аппроксимируем функцию \textbf{квадратичным полиномом} на каждой паре отрезков: \[ \int_a^b f(x)\,dx \approx \frac{h}{3} \left[ f(a) + 4\sum_{k=1,3,5,\dots}^{N-1} f(a+kh) + 2\sum_{k=2,4,6,\dots}^{N-2} f(a+kh) + f(b) \right] \] где $N$ --- чётное число. Погрешность --- $\mathcal{O}(h^4)$! Это уже очень хорошо для гладких функций. \begin{tipbox}[Почему Симпсон так хорош?] Метод Симпсона --- это, по сути, интегрирование кусочно-квадратичных сплайнов. Если функция гладкая (имеет непрерывные производные до 4-го порядка), то погрешность убывает как $h^4$, что намного быстрее, чем у трапеций. Для очень гладких функций есть ещё более точные методы --- квадратуры Гаусса, которые используют оптимально выбранные узлы и веса. \end{tipbox} \subsection{Проклятие размерности} Но что будет, если мы будем интегрировать 2-, 3-, 4-, 10-мерную функцию? Нам надо тогда $N^{10}$ точек интегрирования, где $N$ --- хотя бы 100... Это как-то много! Давайте посчитаем. Для 10-мерного интеграла с $N = 100$ точками по каждой оси: \[ 100^{10} = 10^{20} \text{ точек} \] Даже если каждая точка вычисляется за 1 наносекунду, это займёт: \[ 10^{20} \text{ нс} = 10^{11} \text{ с} \approx 3000 \text{ лет} \] Это --- \textbf{проклятие размерности} (curse of dimensionality). \subsubsection{Где это встречается в химии?} \begin{itemize} \item \textbf{Квантовая химия:} Вычисление многоэлектронных интегралов (электрон-электронное отталкивание) --- это 6-мерные интегралы (по 3 координаты на каждый электрон). \item \textbf{Статистическая механика:} Статистическая сумма --- это интеграл по фазовому пространству всех частиц. Для $N$ частиц --- это $6N$-мерный интеграл (3 координаты + 3 импульса на частицу). \item \textbf{Молекулярная динамика:} Усреднение по конфигурациям --- многомерное интегрирование. \item \textbf{Машинное обучение:} Байесовский вывод --- интегрирование по пространству параметров модели. \end{itemize} \subsection{Метод Монте-Карло для интегрирования} И тут на сцену выходит метод Монте-Карло --- один из самых элегантных и удивительных алгоритмов в вычислительной математике. \subsubsection{Идея на пальцах} Представим, что мы хотим вычислить интеграл: \[ I = \int_{[a,b]^d} f(\vec{x})\,d\vec{x} \] где $d$ --- размерность (может быть очень большой). Метод Монте-Карло говорит: ``Давай мы просто накидаем $N$ случайных точек $\vec{x}_1, \vec{x}_2, \dots, \vec{x}_N$ равномерно в область интегрирования, вычислим в них функцию и усредним''. Формула: \[ I \approx \frac{V}{N} \sum_{k=1}^N f(\vec{x}_k) \] где $V = (b-a)^d$ --- объём области интегрирования, а $\vec{x}_k$ --- случайные точки, равномерно распределённые в $[a,b]^d$. \subsubsection{Почему это работает?} Это --- просто закон больших чисел. Математическое ожидание случайной величины $f(\vec{X})$, где $\vec{X}$ равномерно распределена в $[a,b]^d$, равно: \[ \E[f(\vec{X})] = \frac{1}{V} \int_{[a,b]^d} f(\vec{x})\,d\vec{x} = \frac{I}{V} \] По закону больших чисел, среднее значение $\frac{1}{N}\sum f(\vec{x}_k)$ сходится к $\E[f(\vec{X})]$ при $N \to \infty$. \subsubsection{Скорость сходимости} А вот тут --- магия! Погрешность метода Монте-Карло убывает как: \[ \text{Погрешность} \sim \frac{\sigma}{\sqrt{N}} \] где $\sigma$ --- стандартное отклонение функции $f(\vec{x})$ в области интегрирования. \textbf{Важно:} эта скорость сходимости \textit{не зависит от размерности} $d$! Для сравнения: \begin{itemize} \item Метод Симпсона в 1D: погрешность $\sim h^4 \sim N^{-4}$, \item Метод Симпсона в $d$-мерном случае: погрешность $\sim N^{-4/d}$ (катастрофически медленно при больших $d$), \item Метод Монте-Карло в любой размерности: погрешность $\sim N^{-1/2}$. \end{itemize} Да, Монте-Карло сходится медленнее, чем Симпсон в 1D. Но в 10-мерном случае: \begin{itemize} \item Симпсон: $N^{-4/10} = N^{-0.4}$, \item Монте-Карло: $N^{-0.5}$. \end{itemize} Монте-Карло \textit{быстрее}! \begin{warningbox}[Цена Монте-Карло] Метод Монте-Карло сходится как $N^{-1/2}$ --- это значит, чтобы увеличить точность в 10 раз, нужно в 100 раз больше точек. Это медленно по абсолютным меркам, но это \textit{единственный} метод, который работает в высоких размерностях. \end{warningbox} \subsection{Пример из квантовой химии: вычисление орбиталей} Давайте посмотрим, как метод Монте-Карло применяется в реальной квантовой химии. \subsubsection{Задача: нормировка волновой функции} В квантовой механике волновая функция $\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N)$ описывает состояние $N$ электронов. Она должна быть нормирована: \[ \int |\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N)|^2 \,d\vec{r}_1\,d\vec{r}_2\cdots d\vec{r}_N = 1 \] Это --- $3N$-мерный интеграл! Для молекулы воды ($N = 10$ электронов) --- это 30-мерный интеграл. \subsubsection{Применение Монте-Карло} Мы генерируем $M$ случайных конфигураций электронов: \[ \{\vec{r}_1^{(k)}, \vec{r}_2^{(k)}, \dots, \vec{r}_{10}^{(k)}\}, \quad k = 1, \dots, M \] где каждая координата $\vec{r}_i^{(k)} = (x_i^{(k)}, y_i^{(k)}, z_i^{(k)})$ выбирается из некоторого распределения (например, гауссова). Вычисляем $|\Psi|^2$ в каждой конфигурации и усредняем: \[ \int |\Psi|^2 \approx \frac{V^{30}}{M} \sum_{k=1}^M |\Psi(\vec{r}_1^{(k)}, \dots, \vec{r}_{10}^{(k)})|^2 \] \subsubsection{Вариационный метод Монте-Карло (Variational Monte Carlo, VMC)} Но это ещё не всё! В квантовой химии есть метод \textbf{Variational Monte Carlo (VMC)}, который использует Монте-Карло для \textit{минимизации} энергии. Мы выбираем пробную волновую функцию $\Psi_T(\vec{r}_1, \dots, \vec{r}_N; \vec{\alpha})$ с параметрами $\vec{\alpha}$ и вычисляем энергию: \[ E(\vec{\alpha}) = \frac{\int \Psi_T^* \hat{H} \Psi_T \,d\vec{r}_1\cdots d\vec{r}_N}{\int |\Psi_T|^2 \,d\vec{r}_1\cdots d\vec{r}_N} \] Оба интеграла --- $3N$-мерные, и мы вычисляем их методом Монте-Карло. Затем минимизируем $E(\vec{\alpha})$ по параметрам $\vec{\alpha}$ --- и получаем приближённую волновую функцию. \subsubsection{Quantum Monte Carlo (QMC)} Есть ещё более продвинутые методы --- \textbf{Diffusion Monte Carlo (DMC)}, которые решают уравнение Шрёдингера напрямую, моделируя ``диффузию'' электронов в воображаемом времени. Эти методы дают \textit{почти точные} решения для небольших молекул и используются как эталон для проверки других методов. \begin{tipbox}[Где используется Монте-Карло в химии?] \begin{itemize} \item \textbf{Quantum Monte Carlo (QMC):} Точное решение уравнения Шрёдингера для небольших систем, \item \textbf{Молекулярная динамика Монте-Карло:} Сэмплирование конфигураций в статистической механике, \item \textbf{Интегрирование в DFT:} Численное интегрирование обменно-корреляционного функционала, \item \textbf{Байесовская оптимизация:} Поиск глобального минимума энергии молекулы. \end{itemize} \end{tipbox} \subsection{Улучшения метода Монте-Карло} Базовый метод Монте-Карло --- это просто равномерное случайное блуждание. Но есть много улучшений: \subsubsection{Важностное сэмплирование (Importance Sampling)} Вместо равномерного распределения используем распределение, пропорциональное $|f(\vec{x})|$. Тогда точки будут чаще попадать в области, где функция большая, и погрешность уменьшится. \subsubsection{Метод Метрополиса (Metropolis-Hastings)} Для сложных многомерных интегралов используем марковские цепи: каждая следующая точка зависит от предыдущей. Это позволяет эффективно исследовать пространство, даже если оно очень большое. \subsubsection{Квази-Монте-Карло (Quasi-Monte Carlo)} Вместо случайных точек используем \textbf{детерминированные} последовательности с низким отклонением (последовательности Соболля, Халтона). Они ``равномернее'' покрывают пространство, чем случайные точки, и сходятся быстрее: погрешность $\sim N^{-1} (\log N)^d$. \begin{successbox}[Итог] Метод Монте-Карло --- это единственный практичный способ вычисления многомерных интегралов. Его скорость сходимости не зависит от размерности, что делает его незаменимым в квантовой химии, статистической физике и машинном обучении. Хотя он сходится медленнее, чем детерминированные методы в низких размерностях, в высоких размерностях у него просто нет конкуренции. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Статистика: как отделить информацию от шума} До сих пор мы в основном рассматривали математические задачи, в которых числа считаются известными. Но экспериментальная химия устроена иначе. Если мы измеряем концентрацию вещества, положение спектральной линии или интенсивность сигнала, прибор никогда не сообщает нам ``истинное'' значение с бесконечной точностью. Каждое измерение немного отличается от другого. Поэтому вместо одного числа нам приходится иметь дело с \textbf{случайной величиной}. Пусть мы несколько раз измеряем одну и ту же величину. Получаем: $(x_1,x_2,\ldots,x_N)=\vec{x}$, причём $\vec{x}=\mu+\vec{\varepsilon},$ где \begin{itemize} \item $\mu$ --- неизвестное истинное или среднее значение; \item $\vec{\varepsilon}$ --- случайная ошибка измерения. \end{itemize} Это простое уравнение является одной из основных идей статистики. Статистика начинается там, где мы понимаем: \textit{мы наблюдаем данные, но интересуемся скрытыми величинами.} \subsection{Среднее} Самая простая оценка $\mu$ --- среднее арифметическое: % $$ \mu = \frac{1}{N} \sum_{i=1}^N x_i. $$ % Почему именно оно? Если ошибки имеют среднее около нуля, то среднее арифметическое минимизирует сумму квадратов отклонений: $\displaystyle \min_{\mu} \sum_{i=1}^N |x_i - \mu|^2$, и решение даётся формулой выше. Кстати, это уже объясняет, почему повторение эксперимента обычно повышает точность результата. \subsection{Дисперсия и стандартное отклонение} Среднее говорит нам, где находится центр данных, но ничего не говорит о том, насколько сильно результаты разбросаны. Для этого вводят дисперсию --- меру того, насколько сильно наши измерения ``разбросаны'' вокруг истинного значения (к которому стремится среднее): $$ \sigma^2 = \frac{1}{N} \sum_{i=1}^N |x_i - \mu|^2 $$ Да, мы же только что её минимизировали? Верно, но значение минимума будет равно нулю тогда и только тогда, когда все $x_i$ одинаковы; иначе это значение и есть дисперсия. А чтобы измерять ее в тех же единицах, мы просто вспомним, что вторая норма --- содержит в себе корень, и запишем $$ \sigma = \frac{1}{\sqrt{N}} \sqrt{\sum_{i=1}^N |x_i - \mu|^2} $$ То есть дисперсия --- это, по существу, \textbf{квадрат длины вектора отклонений от среднего}. Статистика снова превращается в линейную алгебру. \subsection{Нормальное распределение} Во многих физических и химических измерениях случайные ошибки приблизительно описываются нормальным распределением: $$ p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right). $$ Не обязательно запоминать эту формулу. Гораздо важнее понять смысл двух параметров: $\mu$ и $\sigma$. Первый определяет положение центра распределения, второй --- его ширину, то есть чем больше $\sigma$, тем сильнее разброс измерений. \begin{tipbox}[Интуиция] \begin{itemize} \item Среднее отвечает на вопрос: \textit{``Где находится результат?''} \item Стандартное отклонение отвечает на вопрос: \textit{``Насколько сильно результаты разбросаны?''} \end{itemize} \end{tipbox} \subsection{Почему среднее становится точнее при повторении эксперимента?} Пусть мы сделали $N$ независимых измерений. Для независимых ошибок дисперсия среднего равна $\displaystyle \frac{\sigma^2}{N}$, следовательно, \textbf{ошибка среднего уменьшается как} $\frac{\sigma}{\sqrt N}$. Поэтому увеличение числа измерений повышает точность, но не линейно (как в Монте-Карло!). \begin{tipbox} Чтобы уменьшить случайную ошибку в $10$ раз, в идеализированном случае понадобится примерно в $100$ раз больше независимых измерений. \end{tipbox} \subsection{Случайная ошибка и систематическая ошибка} Здесь необходимо сделать очень важное различие. Если прибор дает $x=\mu+\varepsilon$, где $\varepsilon$ случайно колеблется около нуля, многократные измерения могут уменьшить влияние этой ошибки. Но представим, что прибор систематически завышает результат: $x=\mu+b+\varepsilon$, где $b$ --- постоянное смещение. Тогда усреднение дает $\bar{x}\approx\mu+b$. И сколько бы раз мы ни повторяли эксперимент, смещение $b$ никуда не исчезнет. \begin{warningbox}[Важно] Повторение измерений уменьшает случайную ошибку, но не устраняет систематическую ошибку. Статистика не может исправить неправильно откалиброванный прибор. \end{warningbox} Это одна из причин, почему в химии так важны калибровка, контрольные образцы и независимые методы измерения. \subsection{Две величины могут быть связаны} Пусть одновременно измеряются две величины: $x_i$ и $y_i$. Например: % \begin{itemize} \item концентрация вещества и интенсивность сигнала; \item температура и скорость реакции; \item давление и объем; \item две спектральные характеристики. \end{itemize} Нас может интересовать вопрос: \textit{изменяются ли $x$ и $y$ согласованно?} Для этого используют ковариацию: $\displaystyle \operatorname{Cov}(x,y) = \frac{1}{N} (x-\mu_x)^T (y-\mu_y)$ Если большие значения $x$ обычно соответствуют большим значениям $y$, ковариация положительна. Если большие значения $x$ соответствуют маленьким $y$, она отрицательна. Если линейной связи нет, ковариация может быть близка к нулю. \subsection{Корреляция} Ковариация зависит от единиц измерения. Поэтому часто используют нормированную величину: $\displaystyle \rho_{xy} = \frac{\operatorname{Cov}(x,y)} {\sigma_x\sigma_y}.$ Для нее $\displaystyle -1\leq\rho_{xy}\leq1$. Значение около $1$ означает сильную положительную линейную связь, значение около $-1$ --- сильную отрицательную, а значение около $0$ означает отсутствие выраженной линейной корреляции. Но здесь необходимо запомнить одну очень важную вещь: \begin{warningbox}[Корреляция не означает причинность] Если две величины коррелируют, это еще не означает, что одна вызывает другую. Корреляция говорит о статистической связи между данными. Чтобы установить причинную связь, необходимы дополнительные физические, химические или экспериментальные аргументы. \end{warningbox} \subsection{Ковариационная матрица} Если у нас не две, а $p$ измеряемых величин, $x_1,x_2,\ldots,x_p$, то ковариации всех пар можно собрать в одну матрицу: $$ \Sigma= \begin{pmatrix} \operatorname{Cov}(x_1,x_1) & \operatorname{Cov}(x_1,x_2) & \cdots\\ \operatorname{Cov}(x_2,x_1) & \operatorname{Cov}(x_2,x_2) & \cdots\\ \vdots & \vdots & \ddots \end{pmatrix}. $$ Эта матрица симметрична: $\Sigma=\Sigma^T$, и она не просто хранит статистическую информацию. Она является математическим объектом линейной алгебры. Именно поэтому собственные значения и собственные векторы, с которыми мы уже встречались раньше, снова возникают в статистике. \subsection{PCA: статистика встречается с линейной алгеброй} Представим, что у нас есть $N$ химических образцов, для каждого из которых измерено $p$ характеристик. Получаем матрицу данных: $X\in\mathbb{R}^{N\times p}$. Некоторые характеристики могут быть сильно связаны между собой. Например, если две измеряемые величины почти всегда меняются вместе, информация о них частично избыточна. Хотелось бы найти новые координаты, в которых: \begin{itemize} \item первая координата содержит максимально возможную вариацию данных; \item вторая --- максимально возможную оставшуюся вариацию; \item третья --- следующую и так далее. \end{itemize} Это и есть основная идея \textbf{метода главных компонент} (Principal Component Analysis, PCA). Математически PCA тесно связан с собственными векторами ковариационной матрицы: $\displaystyle \Sigma v_i=\lambda_i v_i$. Собственные векторы $v_i$ задают новые направления, а собственные значения $\lambda_i$ показывают, сколько вариации данных приходится на соответствующее направление. Если несколько первых собственных значений значительно больше остальных, $$ \lambda_1,\lambda_2,\ldots,\lambda_k \gg \lambda_{k+1},\ldots,\lambda_p, $$ то большая часть информации может быть представлена всего несколькими координатами. Например, $$ 5000\text{ измеряемых параметров} \quad\longrightarrow\quad 20\text{ главных компонент}. $$ Это не означает, что мы ``выбросили'' химическую информацию. Мы нашли более компактное математическое представление данных. И здесь снова появляется SVD, с которым мы уже встречались: $$ X=U\Sigma V^T. $$ PCA и SVD оказываются двумя сторонами одной и той же линейно-алгебраической идеи. \subsection{Регрессия: когда мы хотим построить модель} Предположим, что мы измерили концентрацию $c_i$ и соответствующий сигнал $y_i$. Мы предполагаем линейную модель: $$ y_i=ac_i+b+\varepsilon_i. $$ Поскольку измерения содержат ошибки, точки обычно не лежат точно на одной прямой. Поэтому ищем такие $a$ и $b$, которые минимизируют сумму квадратов ошибок: $$ \min_{a,b} \sum_{i=1}^N \left(y_i-ac_i-b\right)^2. $$ Это и есть \textbf{метод наименьших квадратов}. Таким образом, регрессия не является чем-то совершенно новым. Мы уже знаем эту математическую конструкцию: $$ \boxed{ \text{данные} \rightarrow \text{модель} \rightarrow \text{невязка} \rightarrow \text{минимизация} } $$ Именно поэтому линейная алгебра и оптимизация так важны для статистики. \subsection{Переобучение} Теперь возникает более сложная проблема. Если данных немного, мы можем подобрать очень сложную функцию, которая пройдет практически через каждую экспериментальную точку. Например, вместо прямой можно использовать полином высокой степени: $$ y=a_0+a_1x+a_2x^2+\cdots+a_kx^k. $$ При достаточно большом $k$ можно почти идеально описать имеющиеся измерения. Но это не обязательно означает, что мы нашли правильный физический закон. Мы могли просто подстроиться под шум. Это называется \textbf{переобучением} (overfitting). Именно поэтому в статистике и машинном обучении нас интересует не только вопрос: \textit{``Насколько хорошо модель описывает известные данные?''}, но и: \textit{``Насколько хорошо она работает на новых данных?''} \subsection{Обучение, проверка и тестирование} Если мы строим модель по имеющимся данным, полезно разделить данные на части: $$ \boxed{ \text{training} \quad+\quad \text{validation} \quad+\quad \text{test} } $$ На training set модель обучается. Validation set используется для выбора параметров модели. Test set должен оставаться независимым и используется для окончательной оценки способности модели работать на новых данных. Эта идея будет особенно важна, когда мы дойдем до машинного обучения. \subsection{Что статистика действительно пытается сделать} Можно сказать, что статистика занимается не столько ``подсчетом средних'', сколько более общей задачей: \textbf{извлечь надежную информацию из несовершенных данных.} Мы наблюдаем: $$ \text{данные} = \text{структура} + \text{шум}. $$ Наша задача состоит в том, чтобы по данным восстановить интересующую нас структуру и одновременно оценить, насколько мы можем ей доверять. В самом простом случае это выглядит так: $$ x_1,x_2,\ldots,x_N \longrightarrow \bar{x}\pm\text{uncertainty}. $$ В более сложном случае: $$ X \longrightarrow \text{PCA} \longrightarrow \text{низкоразмерное представление}. $$ А еще сложнее: $$ X \longrightarrow \text{модель} \longrightarrow \text{предсказание} \longrightarrow \text{проверка на новых данных}. $$ Именно отсюда естественным образом вырастает современное машинное обучение. \begin{successbox}[Что важно запомнить] \begin{enumerate} \item Реальное измерение содержит случайную и, возможно, систематическую ошибку. \item Среднее описывает центральное значение данных. \item Стандартное отклонение описывает их разброс. \item Случайная ошибка среднего уменьшается как $1/\sqrt N$. \item Систематическая ошибка усреднением не устраняется. \item Ковариация описывает совместное изменение величин. \item PCA ищет наиболее информативные направления в многомерных данных. \item Регрессия строит математическую модель по измерениям. \item Хорошая модель должна работать не только на известных данных, но и на новых. \end{enumerate} \end{successbox} \subsection{И снова линейная алгебра} И, пожалуй, самое приятное для нас наблюдение состоит в том, что статистика оказалась не таким уж чужим предметом. Мы начали с измерений: $\displaystyle x_1,x_2,\ldots,x_N, $ перешли к векторам: $\displaystyle \vec x, $ затем к нормам: $\displaystyle \|\vec x\|_2, $ к матрицам ковариаций: $\displaystyle \Sigma, $ к собственным значениям: $\displaystyle \Sigma v=\lambda v, $ к SVD: $\displaystyle X=U\Sigma V^T, $ и наконец к оптимизации: $\displaystyle \min \|Ax-b\|_2^2. $ То есть статистика не разрушает нашу математическую картину. Она показывает, как та же самая математика начинает работать с \textbf{реальными, несовершенными и зашумленными данными}. И это именно та ситуация, с которой современный химик сталкивается почти каждый день. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{От SVD к сжатым измерениям: когда данных меньше, чем нужно} \subsection{Практическая задача: разделение спектров} Давайте рассмотрим практическую задачу. Нам в лабораторию принесли массу цветастой органики, под рукой был только VIS/IR спектрометр для флуоресцентных спектров. Сказали, что в семплах есть несколько продуктов реакции с предположительно отличающимися спектрами, и попросили определить как концентрации, так и спектры чистых веществ. Мы подобрали не сильно УФ-источник, чтобы флюоресценция была очень яркой и давала разные спектры для разных веществ. Кажется, это задача для SVD! Ведь если мы поставим каждый такой спектр в свой столбец матрицы $\mathbf{A}$, положим, что у нас гипотетически имеется матрица чистых спектров этих веществ $\mathbf{P}$, и для каждого семпла есть коэффициенты концентрации самих веществ $\mathbf{Q}$, тогда: \[ \mathbf{A} = \mathbf{P}\mathbf{Q} \] и мы можем получить $\mathbf{P}$ как левые сингулярные вектора матрицы $\mathbf{A}$, а $\mathbf{Q}$ --- как правые сингулярные вектора, умноженные слева на диагональ. Мы всегда можем также записать сингулярное разложение как: \[ \forall i, j: \quad a_{ij} = \sum_{r=1}^R d_r u_{ir} v_{rj} \] Более того, если у нас было бы одно единственное вещество ($R=1$), именно так бы это и было. \subsection{Но что-то пошло не так...} Вместо красивых чистых спектров с положительными значениями, мы получили какие-то странные графики в $\mathbf{P}$ --- с отрицательными значениями, осцилляциями, ``призрачными'' пиками. На самом деле --- всё так и должно было быть, спектры так получить \textit{нельзя}. Сами спектры --- это положительно заданные графики, значит скалярное произведение любой пары хотя и может быть очень близким к нулю, не обязательно им является. В сингулярном разложении мы требуем, чтобы все столбцы $\mathbf{U}$ были \textbf{ортогональны} друг другу. А реальные спектры веществ --- не ортогональны! Они просто ``похожи'' или ``не похожи''. Нам нужно что-то ещё, чтобы растащить эти спектры так, чтобы наша математика увидела их по-отдельности. \subsection{Магия третьего измерения: теорема Крускала} Если мы так поставим эксперимент, что получим 3D-данные, вместо 2D как в примере выше, то можно математически доказать, что существует и \textbf{не ортогональное} разложение. Например, если мы сможем снять спектры флюоресценции с несколькими разными длинами волн возбуждения. Тогда у нас уже будут трёхмерные данные: \[ \forall i, j, k: \quad a_{ijk} = \sum_{r=1}^R \alpha_r \, b_{ir} \, c_{jr} \, d_{kr} \] и, по \textbf{теореме Крускала} (Kruskal, 1977), у нас нет требований на ортогональность каждой из таких матриц! \subsubsection{Что говорит теорема Крускала} Теорема Крускала --- это один из краеугольных камней многомерного анализа данных. Она даёт условия \textbf{единственности} разложения тензора (так называемого CANDECOMP/PARAFAC разложения). Для тензора 3-го порядка $\mathcal{A} = \sum_{r=1}^R \vec{a}_r \circ \vec{b}_r \circ \vec{c}_r$ (где $\circ$ --- внешнее произведение векторов), разложение \textbf{единственно} с точностью до перестановки и масштабирования компонентов, если: \[ k_A + k_B + k_C \geq 2R + 2 \] где $k_A$ --- это \textbf{k-ранг} (Kruskal rank) матрицы $\mathbf{A}$, то есть максимальное число $k$ такое, что любое подмножество из $k$ столбцов матрицы $\mathbf{A}$ линейно независимо. \begin{tipbox}[Почему это так важно?] В обычном SVD мы имеем \textit{бесконечно много} разложений $\mathbf{A} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^H$ --- любое ортогональное преобразование внутри даёт новое валидное разложение. Поэтому SVD не может выделить ``физические'' компоненты --- только математически удобные (ортогональные). Теорема Крускала говорит: для тензоров 3-го и более высокого порядка при выполнении условия на k-ранги разложение \textit{единственно}! Это значит, что мы можем восстановить именно те физические компоненты, которые были в данных --- без требования ортогональности. \end{tipbox} \subsubsection{Методы решения: PARAFAC и ALS} На практике разложение тензора ищут итерационными методами. Самый популярный --- \textbf{ALS} (Alternating Least Squares): \begin{enumerate} \item Фиксируем $\mathbf{B}$ и $\mathbf{C}$, решаем задачу наименьших квадратов для $\mathbf{A}$, \item Фиксируем $\mathbf{A}$ и $\mathbf{C}$, решаем для $\mathbf{B}$, \item Фиксируем $\mathbf{A}$ и $\mathbf{B}$, решаем для $\mathbf{C}$, \item Повторяем до сходимости. \end{enumerate} Это --- классическая задача нелинейной минимизации (мы её обсуждали в главе про нелинейные операторы), и она может застревать в локальных минимумах. Поэтому на практике запускают ALS много раз с разными начальными приближениями. \begin{orangebox}[Применимость в химии] PARAFAC/CANDECOMP широко используется в: \begin{itemize} \item Флуоресцентной спектроскопии (EEM --- Excitation-Emission Matrix), \item Хромато-масс-спектрометрии, \item ЯМР-спектроскопии, \item Анализе метаболических данных (метаболомика). \end{itemize} \end{orangebox} \subsection{Многомерные спектры ЯМР} А теперь --- самое интересное для химиков-структурщиков. В ЯМР-спектроскопии белков можно построить очень многомерные спектры: 2D, 3D, 4D и даже 5D. Каждое измерение --- это своя частота (химический сдвиг определённого ядра: $^1$H, $^{13}$C, $^{15}$N). Кросс-пики в таких спектрах показывают, какие ядра находятся близко друг к другу в пространстве, что позволяет восстановить 3D-структуру белка. Но вот проблема: если мы хотим снять 4D-спектр с разрешением $1024 \times 256 \times 256 \times 64$ точек, то общее число точек --- это $\sim 4 \times 10^9$. А каждое измерение требует времени (чтобы дождаться релаксации), и в сумме это --- \textbf{недели} работы спектрометра! Когда-то были времена, когда, чтобы определить все кросс-пики у убиквитина (пептид из 76 аминокислот), приходилось снимать спектр около \textbf{трёх недель}. Белок за это время мог деградировать! \subsection{Sparse sampling: меньше --- значит больше} Но можно снять только случайно разреженные одномерные спектры --- причём даже не 10\%, и не 1\%, а реально очень-очень немного, и применить такое же разложение, только для разреженных данных (sparse data), и получить такие же надёжные результаты! Применяя такие методы ещё 20 лет назад, автор с его коллегами ускорил съёмку такого многомерного спектра без потери информации с трёх недель до \textbf{15 минут}. Идея проста: если мы знаем, что спектр \textit{разрежен} в некотором базисе (то есть состоит из небольшого числа пиков), то нам не нужно измерять все точки. Достаточно измерить случайное подмножество, и затем \textit{восстановить} пропущенные точки через оптимизацию. \subsection{Compressed Sensing: революция в измерениях} Это подводит нас к одной из самых красивых идей в современной математике и обработке сигналов --- \textbf{Compressed Sensing} (сжатые измерения, или compressive sampling). \subsubsection{Классическая теория: Найквист--Шеннон} Традиционная теория дискретизации (Найквиста--Шеннона) говорит: чтобы восстановить сигнал с максимальной частотой $f_{\max}$, нужно измерять его с частотой не менее $2f_{\max}$. Для 4D-ЯМР спектра с разрешением $1024 \times 256 \times 256 \times 64$ это означает, что мы \textit{обязаны} измерить все $4 \times 10^9$ точек. Меньше --- нельзя, иначе будет aliasing (наложение частот). \subsubsection{Революция: сжатые измерения} Но в 2004--2006 годах Донехо (Donoho), Кантор (Candès), Тао (Tao) и другие математики показали: если сигнал \textbf{разрежен} (sparse) в некотором базисе, то его можно восстановить из \textbf{существенно меньшего} числа измерений! Формально: пусть сигнал $\vec{x} \in \R^N$ имеет разреженное представление $\vec{x} = \mathbf{\Psi}\vec{s}$, где $\vec{s}$ имеет только $K \ll N$ ненулевых элементов. Тогда мы можем измерить $\vec{y} = \mathbf{\Phi}\vec{x}$, где $\mathbf{\Phi} \in \R^{M \times N}$ --- матрица измерений с $M \sim K \log(N/K) \ll N$, и восстановить $\vec{x}$ через решение задачи: \[ \min_{\vec{s}} \|\vec{s}\|_1 \quad \text{при условии} \quad \vec{y} = \mathbf{\Phi}\mathbf{\Psi}\vec{s} \] \begin{tipbox} \textbf{Почему L1-норма, а не L0?} Идея в том, что мы хотим найти \textit{самое разреженное} решение (минимизировать $\|\vec{s}\|_0$ --- число ненулевых элементов). Но минимизация L0-нормы --- это NP-трудная задача (перебор всех комбинаций). Магия compressed sensing в том, что при определённых условиях на матрицу $\mathbf{\Phi}$ (так называемая Restricted Isometry Property, RIP), минимизация L1-нормы даёт \textit{точно тот же результат}, что и минимизация L0! А L1-минимизация --- это выпуклая задача, которая решается эффективно (это линейное программирование). \end{tipbox} \subsubsection{Условия применимости} Для успешного восстановления нужно два условия: \begin{enumerate} \item \textbf{Разреженность:} Сигнал должен быть разрежен в некотором базисе (например, спектр ЯМР --- это набор дельта-функций, то есть очень разрежен в частотной области). \item \textbf{Некоррелированность:} Матрица измерений $\mathbf{\Phi}$ должна быть ``некоррелирована'' с базисом $\mathbf{\Psi}$. На практике это означает, что точки измерений должны быть выбраны \textbf{случайно} (или псевдослучайно). \end{enumerate} \subsection{Примеры Compressed Sensing в химии и за её пределами} \subsubsection{Быстрая МРТ (Magnetic Resonance Imaging)} Одно из самых известных применений --- ускорение МРТ в медицине. Классическая МРТ требует длительного сканирования (k-space заполняется построчно). С compressed sensing можно измерить только \textbf{случайные} линии k-space и восстановить полное изображение через L1-минимизацию. Это сокращает время сканирования в 5--10 раз --- критично для пациентов, которые не могут долго лежать неподвижно. \subsubsection{Разреженная ЯМР-спектроскопия} В многомерной ЯМР белков compressed sensing позволяет: \begin{itemize} \item Измерять только случайное подмножество точек в многомерном k-space, \item Восстанавливать полный спектр через L1-минимизацию, \item Сокращать время эксперимента с недель до часов или минут. \end{itemize} Это --- именно то, что автор скрипта делал 20 лет назад, и что сейчас стало стандартом в современных ЯМР-спектрометрах (методы sparse sampling, Poisson-gap sampling, и т.д.). \subsubsection{Single-Pixel Camera (камера Райса)} Удивительный пример: камера, у которой \textbf{нет матрицы пикселей}, а есть только один детектор. Изображение восстанавливается через сжатые измерения с использованием DMD (Digital Micromirror Device). Это работает, потому что реальные изображения разрежены в вейвлет-базисе. \subsubsection{Масс-спектрометрия} В масс-спектрометрии compressed sensing используется для ускорения FT-ICR (Fourier Transform Ion Cyclotron Resonance) и Orbitrap --- можно измерять FID короче и всё равно получить высокое разрешение. \subsubsection{Сжатие данных} JPEG использует дискретное косинусное преобразование (близкий родственник Фурье), и изображения разрежены в этом базисе --- поэтому JPEG так хорошо сжимает. Compressed sensing --- это следующая ступень: мы не просто сжимаем уже измеренные данные, а сразу измеряем меньше. \subsection{Связь с предыдущими главами} Давайте посмотрим, как compressed sensing связывает всё, что мы прошли: \begin{itemize} \item \textbf{Линейная алгебра:} Мы решаем переопределённую систему $\mathbf{\Phi}\mathbf{\Psi}\vec{s} = \vec{y}$ через минимизацию L1-нормы. \item \textbf{Обусловленность:} Матрица $\mathbf{\Phi}$ должна быть хорошо обусловлена и удовлетворять Restricted Isometry Property (RIP). \item \textbf{SVD и тензоры:} Для многомерных данных используем тензорные разложения (PARAFAC), которые дают единственность без ортогональности. \item \textbf{Нелинейная оптимизация:} L1-минимизация --- это выпуклая задача, но не гладкая (в нуле нет производной). Используем методы типа ISTA (Iterative Shrinkage-Thresholding Algorithm) или ADMM. \item \textbf{FFT:} Оператор $\mathbf{\Psi}$ часто --- это преобразование Фурье, и мы используем FFT для быстрого умножения. \end{itemize} \begin{successbox}[Главный вывод] Compressed sensing --- это не просто ещё один алгоритм. Это \textit{философия измерений}: если мы знаем, что сигнал простой (разреженный), то мы можем измерять его существенно меньше, чем требует классическая теория. Это революционизировало МРТ, ЯМР-спектроскопию, масс-спектрометрию и многие другие области. И всё это базируется на красивой математике: выпуклой оптимизации, теории вероятностей (случайные матрицы) и линейной алгебре. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Сжатие информации: от архиваторов до JPEG} \subsection{Два мира сжатия} Сжатие информации --- это одна из самых практичных областей прикладной математики. Мы постоянно сжимаем и разжимаем данные: архивы, картинки, видео, музыку. Но не все сжатия одинаковы. Есть два принципиально разных подхода: \begin{itemize} \item \textbf{Сжатие без потерь (lossless):} Мы можем восстановить исходные данные \textit{бит в бит}. Примеры: ZIP, PNG, FLAC. \item \textbf{Сжатие с потерями (lossy):} Мы жертвуем частью информации ради большего сжатия. Примеры: JPEG, MP3, MP4. \end{itemize} \subsection{Информационная энтропия Шеннона} Прежде чем говорить о методах, давайте поймём, \textit{сколько} вообще можно сжать данные. Клод Шеннон в 1948 году ввёл понятие \textbf{информационной энтропии}. Если у нас есть алфавит из $n$ символов, и символ $i$ встречается с вероятностью $p_i$, то энтропия (среднее количество информации на символ) равна: \[ H = -\sum_{i=1}^n p_i \log_2 p_i \quad \text{[бит]} \] Это --- \textbf{теоретический предел} сжатия без потерь. Мы не можем сжать данные сильнее, чем $H$ бит на символ. \begin{tipbox}[Пример] В английском тексте буква ``e'' встречается чаще всего ($p \approx 0.13$), а ``z'' --- очень редко ($p \approx 0.001$). Если бы мы кодировали все буквы одинаково (8 бит на символ), мы бы тратили слишком много бит на редкие буквы. Энтропия английского текста --- около 4.2 бит на символ, значит, теоретически мы можем сжать его примерно в 2 раза. \end{tipbox} \subsection{Сжатие без потерь: кодирование Хаффмана} Идея проста: замечаем, что какие-то символы и последовательности по 2--3 символа встречаются чаще, заводим таблицу таких символов, и чем чаще встречается какой-то символ или последовательность, тем меньше бит используем для его кодирования. \subsubsection{Алгоритм Хаффмана} \begin{enumerate} \item Подсчитываем частоты всех символов в тексте. \item Строим бинарное дерево: два самых редких символа объединяем в узел с суммарной частотой, повторяем, пока не получим одно дерево. \item Код символа --- это путь от корня до листа (0 --- влево, 1 --- вправо). \end{enumerate} Частые символы оказываются ближе к корню --- их коды короче. Редкие --- дальше, коды длиннее. \begin{tipbox}[Пример] Пусть у нас есть символы A (частота 0.5), B (0.25), C (0.125), D (0.125). Дерево Хаффмана даст коды: \begin{itemize} \item A: 0 (1 бит) \item B: 10 (2 бита) \item C: 110 (3 бита) \item D: 111 (3 бита) \end{itemize} Средняя длина: $0.5 \times 1 + 0.25 \times 2 + 0.125 \times 3 + 0.125 \times 3 = 1.75$ бит на символ --- близко к энтропии! \end{tipbox} \subsubsection{Словарные методы: LZ77, LZ78, LZW} Алгоритмы семейства LZ (Lempel--Ziv) идут дальше: они ищут не только частые символы, но и \textbf{повторяющиеся последовательности}. Идея: если мы уже видели строку ``абракадабра'', то при повторном появлении мы не пишем её заново, а ссылаемся на предыдущее вхождение: ``(назад на 11, длина 11)''. Это основа форматов ZIP, GZIP, PNG. \subsection{Химический пример: поиск белковых последовательностей} А теперь --- живой пример из биологии. У нас есть база данных белковых последовательностей (миллионы аминокислот), и нам нужно найти, есть ли в ней гомологи (похожие белки) к нашему новому белку. Прямое сравнение ``наш белок vs каждый белок в базе'' --- это слишком медленно. Нужен умный алгоритм сжатия и поиска. \subsubsection{BLAST (Basic Local Alignment Search Tool)} BLAST --- это один из самых цитируемых алгоритмов в биологии. Идея: \begin{enumerate} \item Разбиваем наш белок на короткие ``слова'' (k-меры, обычно $k = 3$ для белков). \item Для каждого слова ищем похожие слова в базе (с небольшими мутациями). \item Расширяем совпадения в обе стороны, пока сходство не упадёт ниже порога. \end{enumerate} Это --- по сути, словарный метод сжатия + быстрый поиск по хеш-таблице. BLAST позволяет найти гомологи за секунды, хотя полное выравнивание заняло бы часы. \begin{orangebox}[Почему это важно?] BLAST и его варианты (PSI-BLAST, BLASTP, BLASTN) --- это основа современной биоинформатики. Без них не было бы ни расшифровки генома, ни разработки лекарств, ни эволюционных исследований. \end{orangebox} \subsection{Сжатие с потерями: JPEG и малоранговая аппроксимация} Теперь --- сжатие с потерями. Идея: мы жертвуем частью информации, которую человеческий глаз (или ухо) всё равно не заметит. \subsubsection{JPEG через дискретное косинусное преобразование (DCT)} JPEG работает так: \begin{enumerate} \item Разбиваем картинку на блоки $8 \times 8$ пикселей. \item Для каждого блока применяем \textbf{дискретное косинусное преобразование} (DCT) --- это близкий родственник Фурье. \item DCT раскладывает блок на 64 частотные компоненты (от низких частот --- общий фон, до высоких --- мелкие детали). \item Квантуем коэффициенты: делим на матрицу квантования и округляем. Высокие частоты (мелкие детали) квантуются грубее --- мы их ``теряем''. \item Оставшиеся ненулевые коэффициенты сжимаем без потерь (Хаффман). \end{enumerate} \subsubsection{Малоранговая аппроксимация через SVD} Альтернативный взгляд на сжатие картинок --- через SVD. Пусть у нас есть матрица изображения $\mathbf{A} \in \R^{M \times N}$. SVD: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H = \sum_{k=1}^{\min(M,N)} \sigma_k \vec{u}_k \vec{v}_k^H \] Мы можем аппроксимировать $\mathbf{A}$ через первые $r$ сингулярных значений: \[ \mathbf{A}_r = \sum_{k=1}^r \sigma_k \vec{u}_k \vec{v}_k^H \] Это --- \textbf{малоранговая аппроксимация} ранга $r$. По теореме Эккарта--Янга, это \textit{лучшая} аппроксимация ранга $r$ в смысле L2-нормы. \begin{tipbox}[Степень сжатия] Исходная матрица: $M \times N$ чисел. Аппроксимация ранга $r$: $r(M + N + 1)$ чисел. Если $r \ll \min(M,N)$, то сжатие огромно! Например, для картинки $1000 \times 1000$ и $r = 50$: сжатие в $1000 \times 1000 / (50 \times 2001) \approx 10$ раз. \end{tipbox} \subsection{Вейвлеты: локальные частоты} И DCT, и SVD --- это глобальные преобразования: они раскладывают всю картинку по частотам. Но у картинок есть \textit{локальные} особенности: края объектов, текстуры, точки. \textbf{Вейвлеты} --- это базисные функции, которые локализованы и в пространстве, и в частоте. Они позволяют анализировать сигнал на разных масштабах. \subsubsection{Вейвлет-преобразование} Вместо синусоид (как в Фурье) используем ``маленькие волны'' (wavelets) --- функции, которые быстро затухают. Самый популярный --- \textbf{вейвлет Хаара}: \[ \psi(t) = \begin{cases} 1, & 0 \leq t < 0.5 \\ -1, & 0.5 \leq t < 1 \\ 0, & \text{иначе} \end{cases} \] Вейвлет-преобразование раскладывает сигнал по масштабам (частотам) и позициям. Это позволяет: \begin{itemize} \item Сжимать картинки лучше, чем JPEG (формат JPEG2000 использует вейвлеты), \item Detect edges (края объектов), \item Удалять шум, сохраняя важные детали. \end{itemize} \begin{successbox}[Итог] Сжатие информации --- это баланс между математической теорией (энтропия, SVD, вейвлеты) и психофизикой (что видит глаз, что слышит ухо). От архиваторов до JPEG --- везде одна идея: найти структуру в данных и использовать её для компактного представления. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Сравнение изображений: от свёртки до нейросетей} \subsection{Задача: найти объект на картинке} Как сравнить две картинки и понять, что на них есть один и тот же объект? Мы уже сравнивали две последовательности, строя циркулянтную матрицу (глава про FFT). Наверное, картинки можно также сравнить? Пусть они 2D. Но не всё тут так хорошо, как кажется. Если одна картинка относительно другой просто \textbf{сдвинута} по одной или обеим осям --- да, всё будет работать. Мы можем использовать 2D-свёртку (через FFT) и найти максимум --- это будет позиция объекта. Но если есть \textbf{поворот}? Или \textbf{растяжение}? Или \textbf{изменение освещения}? Простая свёртка уже не поможет. \subsection{Edge detection: выделение границ} Первый шаг к пониманию содержимого картинки --- это выделение \textbf{границ} (edges). Границы --- это места, где интенсивность пикселей резко меняется. \subsubsection{Оператор Собеля} Самый простой способ --- это вычислить градиент интенсивности. Оператор Собеля использует две маски $3 \times 3$ для вычисления производных по $x$ и $y$: \[ G_x = \begin{pmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{pmatrix} * I, \quad G_y = \begin{pmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{pmatrix} * I \] где $*$ --- свёртка, $I$ --- изображение. Модуль градиента: $G = \sqrt{G_x^2 + G_y^2}$. Там, где $G$ большой --- там граница. \subsubsection{Детектор Кэнни (Canny Edge Detector)} Более продвинутый алгоритм: \begin{enumerate} \item Сглаживаем картинку гауссовым фильтром (убираем шум). \item Вычисляем градиент (как у Собеля). \item Подавляем немаксимумы (оставляем только самые сильные границы). \item Применяем гистерезис: если граница выше верхнего порога --- оставляем, ниже нижнего --- удаляем, между --- оставляем, если связана с сильной границей. \end{enumerate} Результат --- тонкие, связные линии границ. \subsection{Особые точки: ключевые точки изображения} Границы --- это хорошо, но нам нужно что-то более ``точечное'', чтобы сравнивать картинки. Для этого ищут \textbf{особые точки} (keypoints, interest points) --- места, которые легко найти и которые устойчивы к преобразованиям. \subsubsection{Детектор Харриса (Harris Corner Detector)} Идея: ищем углы --- места, где граница меняет направление. Для каждого пикселя вычисляем матрицу вторых моментов градиента: \[ \mathbf{M} = \sum_{x,y} w(x,y) \begin{pmatrix} G_x^2 & G_x G_y \\ G_x G_y & G_y^2 \end{pmatrix} \] где $w(x,y)$ --- окно (например, гауссово). Если оба собственных значения $\mathbf{M}$ большие --- это угол. Если одно большое, другое маленькое --- это граница. Если оба маленькие --- это плоская область. \subsubsection{SIFT (Scale-Invariant Feature Transform)} SIFT --- это один из самых популярных алгоритмов (Lowe, 1999). Он находит особые точки, которые устойчивы к: \begin{itemize} \item Масштабу (объект может быть ближе или дальше), \item Повороту, \item Изменению освещения, \item Небольшим изменениям ракурса. \end{itemize} Алгоритм: \begin{enumerate} \item Строим \textbf{пирамиду масштабов}: размываем картинку гауссовым фильтром с разным $\sigma$ и уменьшаем размер. \item Ищем \textbf{экстремумы} в разности гауссиан (DoG --- Difference of Gaussians) --- это приближение лапласиана. \item Для каждой особой точки вычисляем \textbf{дескриптор}: гистограмму градиентов в окрестности $16 \times 16$ пикселей. Это вектор из 128 чисел. \item Сравниваем дескрипторы между картинками (евклидово расстояние). \end{enumerate} \subsubsection{SURF, ORB, AKAZE} Есть много улучшений SIFT: \begin{itemize} \item \textbf{SURF} (Speeded-Up Robust Features) --- быстрее, использует интегральные изображения. \item \textbf{ORB} (Oriented FAST and Rotated BRIEF) --- очень быстрый, свободен от патентов. \item \textbf{AKAZE} --- использует нелинейные диффузионные масштабы. \end{itemize} \subsection{Нахождение трансформации: RANSAC} Мы нашли особые точки на двух картинках и сопоставили их (feature matching). Но не все сопоставления правильные --- есть выбросы (outliers). Как найти трансформацию (поворот, масштаб, сдвиг), которая переводит одну картинку в другую, если у нас есть выбросы? \subsubsection{RANSAC (Random Sample Consensus)} RANSAC --- это элегантный алгоритм для робастной оценки параметров: \begin{enumerate} \item Случайно выбираем минимальное число точек (для аффинной трансформации --- 3 пары точек). \item По этим точкам вычисляем параметры трансформации. \item Считаем, сколько \textit{всех} точек согласуются с этой трансформацией (inliers) --- расстояние меньше порога. \item Повторяем $N$ раз, выбираем трансформацию с максимальным числом inliers. \item Финально уточняем параметры по всем inliers (метод наименьших квадратов). \end{enumerate} \begin{tipbox}[Почему RANSAC работает?] Если доля выбросов не слишком велика (скажем, $< 50\%$), то вероятность случайно выбрать только inliers --- ненулевая. Повторяя достаточно раз, мы почти гарантированно найдём ``правильную'' трансформацию. \end{tipbox} \subsection{Методы роя частиц и оптимизация} Иногда у нас нет явных особых точек, но мы знаем, что одна картинка --- это трансформированная версия другой. Как найти параметры трансформации? Это --- задача оптимизации: мы максимизируем \textbf{метрику сходства} (например, взаимную информацию --- mutual information) по параметрам трансформации. \subsubsection{Particle Swarm Optimization (PSO)} Метод роя частиц --- это эвристический метод оптимизации, вдохновлённый поведением стаи птиц или косяка рыб: \begin{enumerate} \item Инициализируем ``рой'' из $N$ частиц --- каждая частица --- это набор параметров трансформации (например, угол поворота, масштаб, сдвиги по $x$ и $y$). \item Каждая частица имеет ``скорость'' и ``позицию''. \item На каждом шаге частица ``запоминает'' свою лучшую позицию (personal best) и знает лучшую позицию всего роя (global best). \item Скорость частицы обновляется как взвешенная сумма: инерция (продолжать двигаться в том же направлении), когнитивная составляющая (стремление к personal best) и социальная составляющая (стремление к global best). \item Повторяем до сходимости. \end{enumerate} PSO особенно полезен, когда: \begin{itemize} \item Функция сходства негладкая или имеет много локальных максимумов, \item Пространство параметров большое (например, 3D-трансформация с 12 параметрами), \item Мы не можем вычислить градиент (например, взаимная информация не дифференцируема). \end{itemize} \subsection{Нейросети: от ручных признаков к обучаемым} Все методы, которые мы обсудили выше (Собель, Харрис, SIFT) --- это \textbf{ручное конструирование признаков} (hand-crafted features). Мы сами придумываем, что такое ``граница'', ``угол'', ``интересная точка'', и как из них построить дескриптор. Но что, если мы позволим компьютеру \textit{самому} научиться находить признаки? \subsubsection{Свёрточные нейронные сети (CNN)} \textbf{Convolutional Neural Networks} --- это специальный тип нейросетей, созданный для работы с изображениями. Идея: \begin{enumerate} \item \textbf{Свёрточные слои:} Применяем набор обучаемых фильтров (как маски Собеля, но параметры учатся из данных). Каждый фильтр выделяет свой признак --- от простых границ на первых слоях до сложных текстур и частей объектов на глубоких слоях. \item \textbf{Pooling слои:} Уменьшаем размер карты признаков (max pooling --- берём максимум в окне $2 \times 2$). Это даёт инвариантность к малым сдвигам. \item \textbf{Полносвязные слои:} В конце --- обычная нейросеть, которая классифицирует или регрессирует. \end{enumerate} \subsubsection{Иерархия признаков} Что удивительно, CNN сами учатся той же иерархии, которую мы конструировали вручную: \begin{itemize} \item \textbf{Первые слои:} Простые границы, градиенты (как Собель), \item \textbf{Средние слои:} Текстуры, углы, простые формы (как SIFT-дескрипторы), \item \textbf{Глубокие слои:} Части объектов (глаза, колёса, листья), \item \textbf{Последние слои:} Целые объекты (лицо, машина, дерево). \end{itemize} Это --- \textbf{обучение признаков} (feature learning) вместо ручного конструирования. \subsection{Предобученные модели и transfer learning} Обучение CNN с нуля требует миллионов размеченных изображений и дней вычислений на GPU. Но есть хитрость --- \textbf{transfer learning} (перенос обучения). Идея: берём сеть, предобученную на огромной базе данных (например, ImageNet --- 1.4 миллиона изображений, 1000 классов), и используем её как \textbf{экстрактор признаков}. \subsubsection{Популярные архитектуры} \begin{itemize} \item \textbf{VGG (2014):} Простая, глубокая (16--19 слоёв), хорошие признаки. \item \textbf{ResNet (2015):} Очень глубокая (до 152 слоёв) с ``skip connections'' --- решает проблему затухающих градиентов. \item \textbf{EfficientNet (2019):} Оптимизирована по всем параметрам (точность, скорость, размер). \item \textbf{Vision Transformer (ViT, 2020):} Использует архитектуру трансформера (как в NLP) для изображений. \end{itemize} \subsubsection{Как использовать предобученную модель} \begin{enumerate} \item \textbf{Feature extraction:} Берём предобученную сеть, отрезаем последний классификационный слой, и используем предпоследний слой как экстрактор признаков. Получаем вектор из, скажем, 2048 чисел для каждого изображения. Сравниваем векторы (косинусное расстояние). \item \textbf{Fine-tuning:} Берём предобученную сеть и дообучаем её на своих данных (например, на микроскопических снимках клеток). Первые слои ``замораживаем'' (они уже хорошо работают), последние --- обучаем. \item \textbf{Siamese networks:} Две одинаковые сети с общими весами, обучаемые так, чтобы похожие изображения имели близкие векторы признаков, а разные --- далёкие. \end{enumerate} \subsection{Применение в химии и биологии} \subsubsection{Крио-электронная микроскопия (cryo-EM)} В cryo-EM мы получаем тысячи 2D-проекций белковых молекул в случайных ориентациях. Задача: \begin{enumerate} \item Найти все молекулы на микрофотографиях (particle picking) --- это задача детекции объектов, \item Определить ориентацию каждой проекции --- это задача сравнения изображений, \item Реконструировать 3D-структуру --- это обратная задача томографии. \end{enumerate} Современные программы (RELION, cryoSPARC) используют CNN для particle picking и классификации. Это позволило достичь ``resolution revolution'' --- определять структуры белков с атомарным разрешением. \subsubsection{Микроскопия и анализ клеток} В биологических исследованиях нужно: \begin{itemize} \item Считать клетки на микроскопических снимках, \item Классифицировать типы клеток, \item Отслеживать движение клеток во времени, \item Находить аномалии (раковые клетки). \end{itemize} CNN (особенно U-Net для сегментации) --- это стандарт в области. \subsubsection{Хемотомика и drug discovery} В поиске лекарств нужно сравнивать молекулярные структуры. Но молекулы --- это не картинки! Однако их можно представить как 2D-изображения (молекулярные графы, отрисованные на сетке) и использовать CNN для предсказания активности. \subsection{Связь с предыдущими главами} Давайте посмотрим, как сравнение изображений связывает всё, что мы прошли: \begin{itemize} \item \textbf{Линейная алгебра:} Свёртка --- это умножение матрицы на вектор (или тензор). CNN --- это последовательность линейных операций с нелинейностями. \item \textbf{SVD и малоранговая аппроксимация:} CNN можно сжимать через SVD (low-rank factorization свёрточных ядер). \item \textbf{FFT:} Свёртка через FFT --- это $\mathcal{O}(N \log N)$ вместо $\mathcal{O}(N^2)$. В CNN это критично для скорости. \item \textbf{Нелинейная оптимизация:} Обучение CNN --- это минимизация функции потерь (cross-entropy, MSE) через backpropagation + SGD/Adam (варианты градиентного спуска). \item \textbf{Compressed sensing:} В MRI и cryo-EM мы используем сжатые измерения + CNN для реконструкции. \end{itemize} \begin{successbox}[Главный вывод] Сравнение изображений --- это от простого к сложному: \begin{enumerate} \item Свёртка (для сдвигов) --- через FFT, \item Особые точки (SIFT) + RANSAC (для поворотов и масштаба), \item Оптимизация (PSO) для сложных трансформаций, \item Нейросети (CNN) --- для семантического сходства. \end{enumerate} Каждый метод --- это компромисс между универсальностью, скоростью и точностью. И все они базируются на математике, которую мы прошли: линейная алгебра, оптимизация, Фурье-анализ, теория вероятностей. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Машинное обучение: от перцептрона до AlphaFold} \subsection{Что такое машинное обучение?} Машинное обучение (Machine Learning, ML) --- это раздел искусственного интеллекта, где компьютер \textit{сам учится} решать задачи на основе данных, а не по жёстко прописанным правилам. Вместо того чтобы писать программу ``если температура > 100, то вода кипит'', мы даём компьютеру тысячи примеров ``температура --- состояние воды'' и просим найти закономерность. \subsubsection{Три типа обучения} \begin{itemize} \item \textbf{Обучение с учителем (supervised):} Есть размеченные данные (вход $\to$ выход). Задача --- научиться предсказывать выход для новых входов. Примеры: классификация изображений, предсказание свойств молекул. \item \textbf{Обучение без учителя (unsupervised):} Данные без разметки. Задача --- найти структуру, кластеры, скрытые закономерности. Примеры: кластеризация спектров, PCA. \item \textbf{Обучение с подкреплением (reinforcement):} Агент взаимодействует со средой и получает награды/штрафы. Задача --- выучить стратегию, максимизирующую награду. Примеры: игра в шахматы, управление роботом. \end{itemize} \subsection{От перцептрона к нейросетям} \subsubsection{Перцептрон (1958, Розенблатт)} Самая простая ``нейросеть'' --- это перцептрон. Он вычисляет: \[ y = \sigma\left(\sum_{i=1}^n w_i x_i + b\right) \] где $x_i$ --- входы, $w_i$ --- веса, $b$ --- смещение, $\sigma$ --- функция активации (например, ступенька или сигмоида). Перцептрон может решать только \textbf{линейно разделимые} задачи (например, логическое ИЛИ). Для XOR (исключающее ИЛИ) одного перцептрона недостаточно --- нужна многослойная сеть. \subsubsection{Многослойный перцептрон (MLP)} Сеть из нескольких слоёв: входной слой, один или несколько скрытых слоёв, выходной слой. Каждый нейрон связан со всеми нейронами следующего слоя. Формула для $l$-го слоя: \[ \vec{h}^{(l)} = \sigma\left(\mathbf{W}^{(l)} \vec{h}^{(l-1)} + \vec{b}^{(l)}\right) \] где $\mathbf{W}^{(l)}$ --- матрица весов, $\vec{b}^{(l)}$ --- вектор смещений, $\sigma$ --- нелинейная функция активации (ReLU, tanh, sigmoid). \subsubsection{Обучение: backpropagation} Как учить нейросеть? Минимизировать функцию потерь $L$ (например, MSE для регрессии или cross-entropy для классификации) через градиентный спуск. \textbf{Backpropagation} --- это эффективный алгоритм вычисления градиентов через цепное правило. Мы идём от выхода к входу, вычисляя $\partial L / \partial w_{ij}$ для каждого веса, и обновляем веса: \[ w_{ij} \leftarrow w_{ij} - \eta \frac{\partial L}{\partial w_{ij}} \] где $\eta$ --- learning rate. \subsection{Свёрточные нейросети (CNN)} Для изображений MLP неэффективен: слишком много параметров (каждый пиксель связан со всеми нейронами). CNN используют \textbf{свёрточные слои}: маленький фильтр (например, $3 \times 3$) скользит по изображению и вычисляет свёртку. Это даёт: \begin{itemize} \item \textbf{Локальность связей:} Каждый нейрон смотрит только на маленькую область, \item \textbf{Разделение весов:} Один и тот же фильтр применяется ко всему изображению, \item \textbf{Инвариантность к сдвигам:} Объект найдётся где угодно. \end{itemize} Популярные архитектуры: LeNet (1998), AlexNet (2012), VGG (2014), ResNet (2015), EfficientNet (2019). \subsection{Рекуррентные нейросети (RNN)} Для последовательностей (текст, временные ряды) нужны сети с \textbf{памятью}. RNN передают скрытое состояние от одного шага к другому: \[ \vec{h}_t = \sigma\left(\mathbf{W}_h \vec{h}_{t-1} + \mathbf{W}_x \vec{x}_t + \vec{b}\right) \] Проблема обычных RNN --- \textbf{затухающие градиенты} (vanishing gradients): сеть плохо учится на длинных последовательностях. Решения: \begin{itemize} \item \textbf{LSTM} (Long Short-Term Memory, 1997): Добавляет ``ворота'' (gates), которые контролируют поток информации. \item \textbf{GRU} (Gated Recurrent Unit, 2014): Упрощённая версия LSTM. \end{itemize} \subsection{Революция: трансформеры и attention} В 2017 году вышла статья ``Attention Is All You Need'' (Vaswani et al.), которая перевернула NLP и не только. \subsubsection{Механизм внимания (Attention)} Идея: вместо того чтобы сжимать всю последовательность в один вектор (как в RNN), мы позволяем каждому элементу ``смотреть'' на все остальные элементы и взвешивать их по важности. Для входной последовательности $\vec{x}_1, \dots, \vec{x}_n$ вычисляем: \[ \text{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}}\right)\mathbf{V} \] где $\mathbf{Q}$ (query), $\mathbf{K}$ (key), $\mathbf{V}$ (value) --- линейные проекции входов. \subsubsection{Transformer} Архитектура Transformer полностью отказывается от рекуррентности и свёрток. Вместо этого используется только механизм внимания (self-attention) и полносвязные слои. Ключевые компоненты: \begin{itemize} \item \textbf{Multi-Head Attention:} Несколько параллельных механизмов внимания, каждый из которых учится фокусироваться на разных аспектах данных. \item \textbf{Positional Encoding:} Поскольку Transformer не имеет рекуррентности, он не знает порядок элементов. Добавляем позиционные кодировки (синусы/косинусы или обучаемые векторы). \item \textbf{Layer Normalization:} Нормализация активаций для стабильности обучения. \item \textbf{Residual Connections:} ``Skip connections'' для борьбы с затухающими градиентами (как в ResNet). \end{itemize} Transformer стал основой для: \begin{itemize} \item \textbf{BERT} (2018): Bidirectional Encoder Representations from Transformers --- предобучение на маскированных словах. \item \textbf{GPT} (2018--2023): Generative Pre-trained Transformer --- авторегрессивная генерация текста. GPT-3 (175 млрд параметров), GPT-4 (мультимодальный). \item \textbf{T5, BART}: Encoder-decoder архитектуры для translation, summarization. \end{itemize} \subsection{Машинное обучение в химии: эволюция} \subsubsection{Эра QSAR (1990-е -- 2010-е)} QSAR (Quantitative Structure-Activity Relationship) --- это классический подход: вычисляем \textbf{молекулярные дескрипторы} (числовые характеристики молекулы) и строим регрессию/классификацию. Дескрипторы: \begin{itemize} \item Физико-химические: молярная масса, logP (липофильность), полярная поверхность, \item Топологические: индексы связности, формы графа молекулы, \item Электронные: заряды атомов, энергии орбиталей, \item 3D-дескрипторы: моменты инерции, радиусы gyration. \end{itemize} Методы: PLS (Partial Least Squares), Random Forest, SVM (Support Vector Machines). \begin{tipbox}[Пример] Предсказание токсичности молекулы: вычисляем 200 дескрипторов, собираем базу данных из 10\,000 молекул с известной токсичностью, обучаем Random Forest. Точность --- 80--85\%. \end{tipbox} \subsubsection{Графовые нейросети (2015 -- настоящее время)} Молекула --- это \textbf{граф}: атомы --- узлы, связи --- рёбра. Графовые нейросети (Graph Neural Networks, GNN) работают напрямую с графами. \textbf{Message Passing Neural Networks (MPNN):} \begin{enumerate} \item Каждый атом имеет начальное представление (вектор признаков: тип атома, заряд, гибридизация). \item На каждом шаге атомы ``обмениваются сообщениями'' с соседями через связи. \item После $K$ шагов каждый атом ``знает'' о своей окрестности радиуса $K$ связей. \item Финальное представление молекулы --- агрегация всех атомных представлений. \end{enumerate} Популярные архитектуры: \begin{itemize} \item \textbf{GCN} (Graph Convolutional Network): Аналог свёртки для графов. \item \textbf{GAT} (Graph Attention Network): Attention между атомами. \item \textbf{SchNet, DimeNet, SphereNet}: Учитывают 3D-координаты атомов и углы. \end{itemize} \begin{successbox}[Почему GNN так хороши для химии?] \begin{itemize} \item \textbf{Инвариантность к перестановкам:} Порядок атомов не важен, \item \textbf{Учёт структуры:} GNN видят связи и геометрию, \item \textbf{Интерпретируемость:} Можно посмотреть, какие атомы/связи важны для предсказания. \end{itemize} \end{successbox} \subsubsection{Предсказание свойств молекул} Современные must-have модели для химика: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Задача} & \textbf{Модель} & \textbf{Точность} \\ \midrule Энергия молекулы & SchNet, DimeNet++ & MAE $\sim$1 kcal/mol \\ Сольватация & SolTranX & MAE $\sim$0.5 kcal/mol \\ pKa & ChemProp, GNN & MAE $\sim$0.3 \\ LogP & MolCLR, GNN & MAE $\sim$0.2 \\ Toxicity & GraphCL, GNN & AUC $\sim$0.9 \\ Drug-likeness & MolBERT & AUC $\sim$0.85 \\ \bottomrule \end{tabularx} \subsection{Революция AlphaFold: предсказание структуры белков} \subsubsection{Проблема} Предсказание 3D-структуры белка по аминокислотной последовательности --- это одна из grand challenges биологии. Экспериментальные методы (X-ray, cryo-EM, NMR) --- дорогие и медленные. \subsubsection{AlphaFold (2020, DeepMind)} AlphaFold 2 --- это прорыв, который решил проблему предсказания структуры белков с атомарной точностью. Архитектура: \begin{enumerate} \item \textbf{Evoformer:} Transformer, который обрабатывает: \begin{itemize} \item Множественное выравнивание последовательностей (MSA) --- эволюционная информация, \item Pair representation --- попарные расстояния между остатками. \end{itemize} \item \textbf{Structure Module:} Итеративно строит 3D-координаты атомов, минимизируя функцию потерь. \item \textbf{Confidence Estimation:} Предсказывает pLDDT (per-residue confidence) --- насколько модель уверена в каждом остатке. \end{enumerate} Результаты на CASP14 (Critical Assessment of Structure Prediction): \begin{itemize} \item Средний GDT\_TS (Global Distance Test) --- 92.4 (экспериментальная точность --- $\sim$90), \item Для 2/3 белков --- точность $< 1$ Å (атомарная точность). \end{itemize} \subsubsection{AlphaFold 3 (2024)} Расширение на: \begin{itemize} \item Белок--лиганд комплексы, \item ДНК/РНК структуры, \item Пост-трансляционные модификации, \item Ионные взаимодействия. \end{itemize} \subsection{Фундаментальная модель конформеров} \subsubsection{Проблема конформационного пространства} Молекула --- это не статичная структура. Она постоянно ``дышит'', вращается вокруг связей, меняет конформацию. Для drug design критически важно знать не одну структуру, а \textbf{ансамбль} низкоэнергетических конформеров. Классические методы: \begin{itemize} \item \textbf{Molecular Dynamics (MD):} Моделируем движение атомов во времени, но это медленно (наносекунды --- микросекунды). \item \textbf{Monte Carlo:} Случайное блуждание по конформационному пространству, но неэффективно для больших молекул. \item \textbf{Systematic/Rotamer search:} Перебор всех комбинаций торсионных углов, но экспоненциальный рост. \end{itemize} \subsubsection{ML-подход: предсказание конформеров} Современные модели учатся предсказывать распределение конформеров напрямую из данных. \textbf{GeoLDM (Geometric Latent Diffusion Models, 2023):} \begin{enumerate} \item \textbf{Encoder:} Кодирует 3D-конформацию в латентное пространство. \item \textbf{Diffusion Process:} Добавляет шум к латентным векторам (как в DALL-E, Stable Diffusion). \item \textbf{Decoder:} Генерирует новые конформации из шума через обратный диффузионный процесс. \item \textbf{Energy Model:} Фильтрует сгенерированные конформации по энергии. \end{enumerate} \textbf{ConfGF (Conformation Graph Factor, 2022):} \begin{itemize} \item Использует GNN для предсказания 3D-координат, \item Генерирует ансамбль конформеров через sampling, \item Учитывает распределение Больцмана (низкоэнергетические конформеры вероятнее). \end{itemize} \subsubsection{Фундаментальная модель: Uni-Mol (2023)} Uni-Mol --- это ``foundation model'' для молекул, предобученная на миллионах конформаций. Архитектура: \begin{enumerate} \item \textbf{3D Transformer:} Работает с атомными координатами и признаками. \item \textbf{Pre-training tasks:} \begin{itemize} \item Предсказание расстояний между атомами, \item Предсказание углов и диэдральных углов, \item Маскированное предсказание атомов (как BERT), \item Контрастивное обучение (похожие конформеры --- близко). \end{itemize} \item \textbf{Fine-tuning:} Дообучение на конкретных задачах (энергия, свойства, конформеры). \end{enumerate} Результаты: \begin{itemize} \item Предсказание энергии: MAE $\sim$0.5 kcal/mol (лучше DFT для некоторых классов), \item Генерация конформеров: RMSD $< 0.5$ Å для 90\% молекул, \item Предсказание свойств: state-of-the-art на большинстве бенчмарков. \end{itemize} \subsection{Must-have инструменты для современного химика} \subsubsection{Программное обеспечение} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Инструмент} & \textbf{Назначение} \\ \midrule RDKit & Хемоинформатика: дескрипторы, fingerprints, similarity \\ DeepChem & ML для drug discovery, GNN \\ PyTorch Geometric & Графовые нейросети \\ OpenMM & Молекулярная динамика на GPU \\ ASE (Atomic Simulation Environment) & Квантовая химия + ML \\ SchNetPack & SchNet и другие GNN для химии \\ AlphaFold (ColabFold) & Предсказание структуры белков \\ Uni-Mol & Foundation model для молекул \\ \bottomrule \end{tabularx} \subsubsection{Типичный workflow} \begin{enumerate} \item \textbf{Сбор данных:} PubChem, ChEMBL, PDB (Protein Data Bank). \item \textbf{Предобработка:} RDKit для дескрипторов, генерации конформеров. \item \textbf{Моделирование:} GNN для предсказания свойств, AlphaFold для структуры. \item \textbf{Валидация:} Cross-validation, external test set, experimental verification. \item \textbf{Интерпретация:} Attention weights, SHAP values для понимания предсказаний. \end{enumerate} \subsection{Связь с предыдущими главами} Давайте посмотрим, как ML связывает всё, что мы прошли: \begin{itemize} \item \textbf{Линейная алгебра:} Нейросети --- это последовательность матричных умножений $\mathbf{W}\vec{x} + \vec{b}$. Backpropagation --- это цепное правило для матриц. \item \textbf{SVD и тензоры:} Сжатие моделей через low-rank factorization. Tensor decompositions для эффективных вычислений. \item \textbf{FFT:} Свёрточные слои используют FFT для ускорения. Attention через FFT (Linear Attention). \item \textbf{Нелинейная оптимизация:} Обучение нейросетей --- это минимизация loss через SGD, Adam (адаптивные методы), L-BFGS (для fine-tuning). \item \textbf{Compressed sensing:} Sparse training, pruning нейросетей. \item \textbf{Графы и тензоры:} GNN работают с графами молекул. 3D-модели используют тензоры координат. \item \textbf{Монте-Карло:} Diffusion models --- это стохастические процессы. Variational inference --- байесовский подход. \end{itemize} \begin{successbox}[Главный вывод] Машинное обучение --- это не замена классической химии, а \textbf{мощный инструмент}, который: \begin{itemize} \item Ускоряет вычисления в тысячи раз (предсказание свойств vs DFT), \item Открывает новые возможности (предсказание структуры белков, генерация молекул), \item Требует понимания математики (линейная алгебра, оптимизация, теория вероятностей). \end{itemize} Современный химик должен знать: \begin{itemize} \item Основы ML (нейросети, GNN, transformers), \item Инструменты (RDKit, PyTorch, DeepChem), \item Когда ML работает, а когда --- нет (физические ограничения, интерпретируемость). \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Современная вычислительная техника: от транзистора до суперкомпьютера} \subsection{Цифры, которые стоит знать} Современный компьютер --- рабочая станция --- это: \begin{itemize} \item 128--1024 ГБ оперативной памяти, \item $\sim$500 Гфлоп/с пиковой производительности (миллиардов операций с плавающей точкой в секунду), \item Несколько процессоров (от 4 до 128 ядер). \end{itemize} Но вот парадокс: скорость доступа к памяти в \textbf{сотни раз} медленнее, чем скорость арифметических операций. А если доступ случайный (не последовательная запись-чтение), то ещё в \textbf{сотню раз} медленнее. Почему так? Давайте разбираться. \subsection{Как устроен процессор} Все говорят, что в процессоре ``много транзисторов''. Но что они делают? На самом деле, процессор --- это такая сущность, которая обращается к \textbf{инструкциям} и обращается к \textbf{данным}. И то, и другое --- занимает место в памяти. \subsubsection{Команды процессора} Процессор зачитывает команды, которые обычно позволяют: \begin{itemize} \item Обратиться к данным в памяти и поместить их в \textbf{регистры} --- временную очень-очень быструю память, \item Выполнить операцию: сложение, вычитание, умножение, сравнение, \item Совершить \textbf{переход}: по умолчанию процессор исполняет команду за командой, но если он попадает на инструкцию перехода, то он может скакнуть вперёд или назад. \end{itemize} Переходы бывают: \begin{itemize} \item \textbf{Безусловные:} Переход происходит всегда, \item \textbf{Условные:} Если мы до этого сравнивали числа и имеем ответ ``да'' или ``нет'', то переход происходит только при выполнении условия. \end{itemize} \subsubsection{SIMD: одна команда --- много данных} Современная процессорная команда может работать сразу с несколькими числами. Например, одна инструкция AVX-512 может выполнить 16 пар сложений одновременно! Это называется \textbf{SIMD} (Single Instruction, Multiple Data) --- одна инструкция, множество данных. Если данные уже есть на самом процессоре, это довольно легко сделать. \begin{tipbox}[Откуда взялся SIMD?] Люди заметили, что часто надо выполнять одинаковые операции: например, одинаково трансформировать пиксели для отрисовки картинки, одинаково обрабатывать элементы массива. Зачем выполнять одну и ту же инструкцию 16 раз, если можно один раз --- но сразу для 16 чисел? \end{tipbox} Но вот чтобы каждой командой выполнять такие операции с \textit{разными} данными, нам нужна была бы очень мощная система по перекачке данных из памяти в процессор и обратно. \subsubsection{Проблема memory wall} К сожалению, такая система заняла бы в сотни и тысячи раз больше транзисторов, чем нужно для арифметического блока. А люди заметили, что в программах мы часто работаем с каким-то небольшим набором чисел, а потом переходим на другой набор, и так далее. Сделать небольшой блок памяти с быстрым доступом --- не так накладно по числу транзисторов. Вот только заставить программиста таскать туда-сюда данные в этот блок из оперативной памяти и обратно --- будет сложно. Поэтому люди автоматизировали это --- и появился \textbf{процессорный кэш}. Это быстрая память, её немного, но процессор сам ``запоминает'', что ты используешь сейчас, и некоторое время держит эти данные у себя в кэше, с расчётом, что далее они, может быть, будут тебе снова нужны. \subsection{Иерархия памяти} Современный компьютер --- это многослойная структура памяти: \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Уровень} & \textbf{Размер} & \textbf{Скорость} & \textbf{Описание} \\ \midrule Регистры & 32--128 $\times$ 64 бит & $\sim$0.3 нс & Ячейки внутри процессора, с которыми он работает напрямую \\ L1 кэш & $\sim$32--64 КБ & $\sim$1 нс & Отдельно для инструкций и данных, свой для каждого ядра \\ L2 кэш & $\sim$256 КБ -- 1 МБ & $\sim$3--5 нс & Свой для каждого ядра или общий на пару ядер \\ L3 кэш & $\sim$8--64 МБ & $\sim$10--20 нс & Общий для всех ядер процессора \\ ОЗУ (DRAM) & 128--1024 ГБ & $\sim$50--100 нс & Основная память \\ Диск (SSD) & 1--10 ТБ & $\sim$10--100 мкс & Долговременное хранение \\ \bottomrule \end{tabularx} Давайте сравним, сколько за одну наносекунду в среднем может сделать современная рабочая станция: \begin{itemize} \item \textbf{Арифметика:} Все её процессоры, с учётом длинных SIMD-инструкций, могут выполнить около \textbf{1000 операций} с плавающей точкой. \item \textbf{L1 кэш:} Можно прочитать и записать по 2--3 числа на каждом процессоре, суммарно $\sim$200 чисел. \item \textbf{L2/L3 кэш:} Эта величина падает в десятки раз. \item \textbf{ОЗУ:} Только единицы чтений и записей в наносекунду, и то --- если происходит большими блоками по 512 байт примерно последовательно. \item \textbf{Случайный доступ:} Если каждый раз обращаемся к разным блокам, скорость падает ещё в десятки раз --- может потребоваться несколько наносекунд на одно число. \end{itemize} \begin{warningbox}[Memory Wall] Случайный доступ к памяти в \textbf{десятки тысяч раз} медленнее реальных вычислительных мощностей современного процессора! Это --- главный bottleneck (``бутылочное горлышко'') современных вычислений. Именно поэтому: \begin{itemize} \item Матричные алгоритмы работают быстрее, если данные расположены ``плотно'' в памяти, \item Разреженные матрицы --- медленные (много случайных обращений), \item Кэш-ориентированные алгоритмы --- это отдельная наука. \end{itemize} \end{warningbox} \subsection{Графические процессоры (GPU)} Но людям этого было мало. Когда они рисовали на экране 3D-объекты компьютерной графики, они заметили, что всё это можно сделать почти полностью \textbf{параллельно}. Условно: если бы у нас было 1024 процессора, мы бы разбили весь экран на сетку квадратиков размера $32 \times 32$, и каждый процессор рисовал бы в своём квадратике. Так появились графические карты --- в них есть уйма небольших специализированных процессоров, которые работают полностью параллельно. \subsubsection{От графики к вычислениям} В конце 90-х такие графические карты использовались для ускорения отображения графики. Примерно в 2005 году люди стали выпускать такие карты так, чтобы на них можно было исполнять небольшие куски кода на тысячах таких небольших процессоров --- так появились \textbf{CUDA} (NVIDIA) и \textbf{OpenCL} (открытый стандарт). К 2010 году все численные методы начали перетекать на графические карты. Постепенно люди к single и double арифметике с плавающей точкой добавляли так называемые \textbf{половинчатые} (FP16) и \textbf{четвертинные} (FP8, INT8) плавающие точки --- так как с ними было быстрее выполнять вычисления и они меньше занимали памяти. В современных NVIDIA графических картах (H100, B200) имеется возможность выполнить около \textbf{миллиона операций} с такими маленькими объектами за ту же одну наносекунду --- и это позволило использовать такие графические карты для современных систем искусственного интеллекта. \subsubsection{Архитектура GPU} Современный GPU (например, NVIDIA H100) --- это: \begin{itemize} \item \textbf{Streaming Multiprocessors (SM):} $\sim$132 мультипроцессора, \item \textbf{CUDA cores:} $\sim$16\,896 ядер для FP32 арифметики, \item \textbf{Tensor cores:} $\sim$528 специализированных ядер для матричных операций (FP16, FP8, INT8), \item \textbf{Память:} 80 ГБ HBM3 (High Bandwidth Memory) с пропускной способностью $\sim$3 ТБ/с. \end{itemize} \subsubsection{Треды и warps} GPU работает с \textbf{тредами} (потоками). Тысячи тредов выполняются параллельно, но организованы они в группы: \begin{itemize} \item \textbf{Warp:} Группа из 32 тредов, которые выполняются \textbf{синхронно} --- все 32 треда выполняют одну и ту же инструкцию в один и тот же момент времени. \item \textbf{Block:} Группа варпов (до 1024 тредов), которые могут обмениваться данными через общую \textbf{shared memory}. \item \textbf{Grid:} Все блоки, запущенные на GPU. \end{itemize} \begin{orangebox}[Правило эффективности GPU] Чтобы GPU работал эффективно, программы на каждом мультипроцессоре для каждого треда должны \textbf{не расползаться}. Это значит: \begin{itemize} \item Все 32 треда в варпе должны выполнять \textit{одну и ту же} инструкцию (без условных переходов, которые разделяют треды --- ``warp divergence''), \item Все треды должны обращаться к \textit{последовательным} адресам памяти (coalesced memory access), \item Должно быть достаточно тредов, чтобы ``спрятать'' задержки памяти. \end{itemize} Если эти условия не выполнены --- GPU работает в разы медленнее. \end{orangebox} \subsection{Параллельные вычисления: таксономия Флинна} В 1966 году Майкл Флинн предложил классификацию параллельных архитектур по потокам инструкций и данных: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Тип} & \textbf{Описание} & \textbf{Пример} \\ \midrule \textbf{SISD} & Single Instruction, Single Data. Один поток инструкций, один поток данных. Классический последовательный процессор. & Старые CPU \\ \textbf{SIMD} & Single Instruction, Multiple Data. Одна инструкция применяется к множеству данных. & Векторные инструкции (AVX), GPU \\ \textbf{MISD} & Multiple Instruction, Single Data. Разные инструкции к одним данным. Редко встречается. & Некоторые специализированные архитектуры \\ \textbf{MIMD} & Multiple Instruction, Multiple Data. Множество процессоров, каждый выполняет свои инструкции над своими данными. & Современные мультипроцессоры, кластеры \\ \bottomrule \end{tabularx} \subsubsection{Общая и распределённая память} В MIMD-системах есть два подхода: \begin{itemize} \item \textbf{Общая память (shared memory):} Все процессоры имеют доступ к единой памяти. Пример: многоядерный CPU. Легко программировать (все видят одни данные), но сложно масштабировать (конфликты доступа к памяти). \item \textbf{Распределённая память (distributed memory):} Каждый процессор имеет свою локальную память, обмен данными --- через сеть. Пример: кластеры, суперкомпьютеры. Сложнее программировать (нужно явно передавать данные --- MPI), но масштабируется до миллионов ядер. \end{itemize} \subsection{TOP500: самые мощные суперкомпьютеры мира} Список TOP500 --- это рейтинг 500 самых мощных суперкомпьютеров мира, обновляемый дважды в год (июнь и ноябрь). Производительность измеряется в LINPACK benchmark --- решении большой системы линейных уравнений. \subsubsection{Современные лидеры (2024--2025)} \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Имя} & \textbf{Место} & \textbf{Производительность} & \textbf{Архитектура} \\ \midrule Frontier & \#1 & $\sim$1.2 Эфлоп/с & AMD EPYC + AMD Instinct MI250X \\ Aurora & \#2 & $\sim$1 Эфлоп/с & Intel Xeon + Intel Max PVC \\ Eagle & \#3 & $\sim$561 Пфлоп/с & AMD EPYC + NVIDIA H100 \\ Fugaku & \#4 & $\sim$442 Пфлоп/с & Fujitsu A64FX (ARM) \\ LUMI & \#5 & $\sim$231 Пфлоп/с & AMD EPYC + AMD MI250X \\ \bottomrule \end{tabularx} \begin{tipbox}[Масштаб] 1 Эфлоп/с = $10^{18}$ операций в секунду. Frontier --- это $\sim$1.2 миллиона миллиардов операций в секунду. Для сравнения: ваш ноутбук --- $\sim$0.1 Тфлоп/с = $10^{11}$ операций в секунду. Разница --- в 10 миллионов раз! \end{tipbox} \subsubsection{Как устроен современный суперкомпьютер} Типичная архитектура: \begin{enumerate} \item \textbf{Узлы (nodes):} $\sim$10\,000--100\,000 серверов, каждый с 2--4 CPU + 4--8 GPU, \item \textbf{Сеть:} Высокоскоростная interconnect (InfiniBand, Slingshot) с пропускной способностью 200--400 Гбит/с на узел, \item \textbf{Хранилище:} Параллельная файловая система (Lustre, GPFS) с пропускной способностью $\sim$1--10 ТБ/с, \item \textbf{Охлаждение:} Жидкостное охлаждение (вода или даже двухфазное), так как мощность --- 20--50 МВт. \end{enumerate} \subsection{Оцифровщики: мост между аналоговым и цифровым миром} Но процессор --- он ведь где-то должен брать для вычисления данные. И тут на сцену выходят \textbf{аналого-цифровые преобразователи} (ADC --- Analog-to-Digital Converters). Они определяются двумя параметрами: \begin{itemize} \item \textbf{Точность (разрядность):} Сколько бит используется для кодирования одного отсчёта, \item \textbf{Скорость (частота дискретизации):} Сколько отсчётов в секунду. \end{itemize} \subsubsection{Компромисс: точность vs скорость} Есть фундаментальный компромисс: чем выше точность, тем ниже скорость. \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Разрядность} & \textbf{Скорость} & \textbf{Минимальный бит} & \textbf{Применение} \\ \midrule 12 бит & $\sim$10 Гвыб/с & $\sim$500 мкВ & Осциллографы, быстрая съёмка \\ 16 бит & $\sim$1 Гвыб/с & $\sim$10--50 мкВ & Аудио, ЯМР, точные измерения \\ 24 бит & $\sim$4 Мвыб/с & $\sim$100--500 нВ & Аудио высокого разрешения, сейсмография \\ \bottomrule \end{tabularx} \begin{warningbox}[Почему не измеряют нановольты напрямую?] 24 бита дают диапазон десятков и сотен нановольт. Но нановольты никто не измеряет напрямую --- радиолокационные помехи от всего вокруг составляют около десятков микровольт! То есть измеряют там обычно \textit{токи}, а не вольты --- или используют специальные экранированные камеры. \end{warningbox} \subsection{Физика переключений: от транзистора до мегавольт} В современном процессоре всё работает, включая-выключая транзисторы на примерно 0.7 В. Скорость такого переключения --- десятки гигагерц. Но вот токи там --- очень маленькие (микроамперы на транзистор). \subsubsection{Компромисс: напряжение vs скорость} Если увеличивать напряжение, то скорость переключения начинает падать: \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Напряжение} & \textbf{Время переключения} & \textbf{Технология} \\ \midrule 0.7 В & $\sim$0.05 нс (20 ГГц) & Современные CMOS транзисторы \\ 40 В & $\sim$3--5 нс & Силовые MOSFET (обычные номиналы) \\ 1000 В & $\sim$3--5 нс & Силовые MOSF (топовые номиналы) \\ 4.5--5 кВ & $\sim$20--30 нс & IGBT, высоковольтные MOSFET \\ $>$10 кВ & $\sim$100 нс -- 1 мкс & Один полупроводник уже не справляется \\ 1 МВ & $\sim$100--500 мкс & Лампы, сборки полупроводников, умножители \\ \bottomrule \end{tabularx} \textit{Хотя современные CMOS транзисторы переключаются на частотах около 20 ГГц, для выполнения одного процессорного такта надо иметь запас в несколько таких переключений, поэтому тактовая частота процессоров находится в диапазоне около 3-4 ГГц} \subsubsection{Предел нарастания напряжения} Выше $\sim$5 кВ один полупроводник уже не справляется --- надо переходить на: \begin{itemize} \item \textbf{Старые добрые лампы} (!) --- вакуумные триоды, тиратроны, \item \textbf{Сборки полупроводников} --- последовательное включение нескольких транзисторов, \item \textbf{Умножители напряжения} --- каскадные схемы с лампами или разрядниками. \end{itemize} Создать даже 1 мегавольт --- это не сильно сложно (разрядник Маркса, каскады Кокрофта--Уолтона). Другое дело --- включить-выключить этот мегавольт за наносекунду. На это нужны реально сотни микросекунд. \begin{successbox}[Фундаментальный предел] Фактически $\sim 10^{12}$ В/с --- это современный предел нарастания напряжения практически во всём диапазоне возможных напряжений. Это --- физическое ограничение, связанное с: \begin{itemize} \item Скоростью движения носителей заряда в полупроводниках (предел $\sim 10^7$ см/с), \item Паразитными ёмкостями и индуктивностями, \item Скоростью распространения электромагнитных волн в среде. \end{itemize} \end{successbox} \subsection{Связь с предыдущими главами} Давайте посмотрим, как всё, что мы прошли, связано с ``железом'': \begin{itemize} \item \textbf{Линейная алгебра:} Матричные операции --- это основа всех вычислений. GPU оптимизированы именно под них (Tensor Cores). \item \textbf{FFT:} Быстрое преобразование Фурье --- это $\mathcal{O}(N \log N)$ операций, и оно критично для ЯМР, обработки сигналов, сжатия данных. GPU ускоряют FFT в десятки раз. \item \textbf{Итерационные методы:} Метод сопряжённых градиентов, GMRES --- это основа решения больших разреженных систем. Они работают на суперкомпьютерах с распределённой памятью (MPI + OpenMP + CUDA). \item \textbf{Машинное обучение:} Обучение нейросетей --- это миллиарды матричных умножений. Без GPU и Tensor Cores современные модели (GPT-4, AlphaFold) были бы невозможны. \item \textbf{Compressed sensing:} L1-минимизация --- это итерационный метод, который требует тысяч итераций. GPU ускоряют его в 100--1000 раз. \item \textbf{Оцифровщики:} ADC --- это мост между аналоговым миром (спектры, сигналы) и цифровым (процессор). Точность ADC определяет, какую математику мы можем применить. \end{itemize} \begin{successbox}[Главный вывод] Современная вычислительная техника --- это: \begin{itemize} \item \textbf{Иерархия памяти:} Регистры $\to$ L1 $\to$ L2 $\to$ L3 $\to$ ОЗУ $\to$ диск. Каждая ступень --- в 10--100 раз медленнее, но в 10--1000 раз больше. \item \textbf{Параллелизм:} SIMD (векторные инструкции), многоядерность (CPU), массовый параллелизм (GPU), кластеры (суперкомпьютеры). \item \textbf{Специализация:} Tensor Cores для AI, FPGA для специфичных задач, ASIC для майнинга. \item \textbf{Физические ограничения:} Memory wall, power wall, speed of light --- всё это ограничивает рост производительности. \end{itemize} Для химика это значит: \begin{itemize} \item Понимать, где ``узкое горлышко'' (память или вычисления), \item Уметь писать код, который эффективно использует кэш и GPU, \item Знать, когда задача решается на рабочей станции, а когда --- нужен суперкомпьютер. \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Послесловие: Как пользоваться этими знаниями} Ранее, после прочтения аналогичной книги, когда вы сталкиваетесь с реальной задачей, возникала необходимость имплементировать решение. Для этого надо было довольно хорошо владеть одним или несколькими языками программирования и понимать, как пользоваться сторонними программными комплексами и пакетами. В современный век развития искусственного интеллекта всё это можно делать с поддержкой систем ИИ --- и это будет значительно эффективнее. Но для этого надо немного понимать ключевые применимости языков программирования и так называемый \textbf{workflow} --- как происходит разработка программ. \subsection{Два мира языков программирования} Хотя языков программирования существует огромное количество, их можно разделить на два фундаментально разных класса: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Компилируемые} & \textbf{Интерпретируемые / байт-код} \\ \midrule C, C++, CUDA, Fortran, Rust & Python, JavaScript, Java, C\#, R \\ \midrule Код компилируется в машинные инструкции под конкретный процессор & Код интерпретируется на лету или компилируется в байт-код \\ \midrule Максимальная производительность, прямой доступ к ``железу'' & Гибкость, кроссплатформенность, быстрая разработка \\ \midrule \textbf{Когда использовать:} & \textbf{Когда использовать:} \\ \quad --- Тяжёлые вычисления (DFT, MD, ML) & \quad --- Прототипирование, анализ данных \\ \quad --- GPU-ускорение (CUDA) & \quad --- Визуализация, отчётность \\ \quad --- Встраиваемые системы & \quad --- Веб-интерфейсы, скрипты \\ \bottomrule \end{tabularx} \subsubsection{Почему интерпретируемые языки --- это удобно} Интерпретируемые языки позволяют достигать большей гибкости при переходе от одной платформы к другой. Например, мы открываем одну и ту же веб-страницу на десктопе, на смартфоне и на планшете --- и видим одинаковый контент. Этот контент обычно генерируется с использованием интерпретируемого языка JavaScript, и если бы мы должны были посылать каждый раз свой код под каждый из процессоров, то сервер должен был бы иметь все варианты такого кода заранее. Это иногда даже невозможно предусмотреть --- так как даже разные Android-смартфоны могут иметь различную процессорную архитектуру. \subsubsection{Типичный workflow химика} Для того чтобы быстро отобразить какие-то результаты, часто проще воспользоваться именно интерпретируемыми языками программирования (Python --- абсолютный стандарт в науке). А вот когда задача бывает сложной и требует много вычислительных ресурсов, возникает необходимость использовать именно компилируемые языки программирования. Типичный workflow выглядит так: \begin{enumerate} \item \textbf{Прототип на Python:} Быстро проверить идею, визуализировать данные, понять, работает ли алгоритм. \item \textbf{Оптимизация:} Если код работает, но медленно --- переписать ``узкие места'' на C++/CUDA или использовать готовые библиотеки (NumPy, SciPy, PyTorch). \item \textbf{Продакшн:} Если задача решается регулярно --- оформить в пакет, добавить тесты, документацию. \end{enumerate} \subsection{Работа с ИИ-ассистентами} Сейчас практически любой алгоритм можно имплементировать, написав правильный промпт с помощью чатов. Но стоит заметить: каждый чат --- это почти как человек. И если ему подать задачу не до конца продуманную, он может и намудрить, или, не поняв, запрограммировать что-то не совсем то. \subsubsection{Золотые правила работы с ИИ} \begin{enumerate} \item \textbf{Разбивайте задачу на небольшие блоки.} Не просите ``напиши программу для решения уравнения Шрёдингера''. Просите: ``напиши функцию, которая вычисляет матрицу Фока для заданного базиса''. \item \textbf{Просите тесты.} Для каждого блока просите чат писать тестовые проверочные программы, прогонять такие тесты и убеждаться, что всё в порядке. \item \textbf{Объединяйте постепенно.} Только после того как каждый блок работает, объединяйте их для получения окончательного результата. \item \textbf{Требуйте объяснений.} Если чат написал что-то, что вы не понимаете --- попросите объяснить. Если объяснение непонятно --- попросите проще. Это ваш скрипт, вы должны понимать каждую строчку. \item \textbf{Проверяйте на известных примерах.} Если вы решаете систему уравнений --- проверьте на системе, для которой знаете ответ. Если делаете Фурье-преобразование --- проверьте на синусоиде с известной частотой. \end{enumerate} \begin{warningbox}[Типичные ошибки] \begin{itemize} \item \textbf{Слепое доверие:} ИИ может написать код, который ``выглядит правильно'', но содержит тонкие ошибки (например, неправильная нормировка, неверные границы массивов, утечка памяти). \item \textbf{Игнорирование контекста:} ИИ не знает вашу конкретную задачу --- вы должны объяснить её подробно, с примерами, с ограничениями. \item \textbf{Отсутствие тестов:} Код без тестов --- это код, который сломается в самый неподходящий момент. \end{itemize} \end{warningbox} \subsection{Must-have инструменты для химика} Вот минимальный набор инструментов, которые вы должны знать: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Язык} & \textbf{Применение} \\ \midrule \textbf{Python} & Абсолютный стандарт: анализ данных, ML, визуализация, скрипты \\ \textbf{R} & Статистический анализ, биоинформатика \\ \textbf{MATLAB} & Инженерные расчёты, обработка сигналов (альтернатива Python) \\ \textbf{C++/CUDA} & Высокопроизводительные вычисления, GPU \\ \textbf{Bash/Shell} & Автоматизация, работа с кластерами \\ \bottomrule \end{tabularx} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Библиотека} & \textbf{Назначение} \\ \midrule \textbf{NumPy, SciPy} & Линейная алгебра, оптимизация, интеграция \\ \textbf{Pandas} & Работа с табличными данными \\ \textbf{Matplotlib, Seaborn} & Визуализация \\ \textbf{PyTorch, TensorFlow} & Машинное обучение, нейросети \\ \textbf{RDKit} & Хемоинформатика, молекулы \\ \textbf{ASE, PySCF} & Квантовая химия \\ \textbf{OpenMM, GROMACS} & Молекулярная динамика \\ \bottomrule \end{tabularx} \subsection{Финальный совет} Математика --- это не набор формул, которые надо зазубрить. Это --- \textbf{способ мышления}. Это умение видеть структуру в хаосе, находить закономерности, строить модели, проверять гипотезы. Когда вы столкнётесь с реальной задачей --- не важно, будет ли это предсказание структуры белка, анализ ЯМР-спектра, оптимизация синтеза или что-то совсем новое --- вспомните этот скрипт. Вспомните, что: \begin{itemize} \item Любая задача --- это либо система уравнений, либо задача оптимизации, либо задача аппроксимации, \item Для любой задачи есть проверенные математические методы, \item Современные инструменты (Python, GPU, ИИ) позволяют решать задачи, которые 20 лет назад были невозможны, \item Главное --- понять суть задачи, а остальное --- техника. \end{itemize} Мы верим в вас. У вас всё получится. И если когда-нибудь этот скрипт поможет вам решить задачу, которая изменит мир --- мы будем невероятно горды. \begin{successbox}[Последнее] Не бойтесь ошибаться. Не бойтесь спрашивать. Не бойтесь экспериментировать. Математика --- это не экзамен, это приключение. И вы --- в начале самого интересного пути. Удачи! \end{successbox} \vfill \begin{center} \textit{С любовью, \\ Папа и Мама} \\[0.5cm] \textit{Сентябрь 2026} \end{center} \end{document}