\documentclass[11pt,english]{article} \usepackage[english]{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} % --- Geometry and layout --- \usepackage{titlesec} \usepackage{fancyhdr} \pagestyle{fancy} \fancyhf{} \fancyhead[L]{\leftmark} \fancyfoot[C]{\thepage} % --- Colors and graphics --- \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} % --- Hyperlinks (loaded after geometry/titlesec/fancyhdr) --- \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=tumblue, urlcolor=tumblue, citecolor=tumblue } % --- User commands --- \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}{} % --- Boxes for important notes --- \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} Mathematics, Computational Methods \\ and Signal Processing}\\[0.5cm] \large A script for students of the TUM Department of Chemistry} \author{Ilgis Ibragimov \& Elena Ibragimova \& AI\footnote{Editing and Translation Assistant}} \date{September 2026 \\ DOI: 10.5281/zenodo.23002397} \maketitle \thispagestyle{empty} \vfill \begin{center} \textit{Dedicated to our daughters Maria and Alevtina. \\ May mathematics become for you not a collection of formulas, \\ but the language in which the Universe is written.} \end{center} \newpage \tableofcontents \newpage \section{Preface} Dear daughters, We so wanted everything in your life's path to be clear and simple. We know that mathematics is often taught by deliberately confusing the learner with terminology and complicated formulas --- so that the forest cannot be seen for the trees. This script is our attempt to give an overview of the mathematical and computational ideas that are especially useful for a modern chemist working with experimental data, modelling and computations, so that you have a \textit{clear picture} of the modern level of mathematical knowledge. Not all knowledge --- that is impossible in one book --- but precisely that which is used in chemistry, biochemistry, spectroscopy, molecular modelling and data processing. We have travelled a long way together: \begin{itemize} \item From the very foundations --- sets of numbers, complex numbers, vectors and matrices, \item Through linear algebra --- SVD, conditioning, methods for solving systems, \item To nonlinear problems --- optimisation, gradients, Newton's method, \item Through signal processing --- Fourier, Prony, compressed sensing, \item To numerical methods --- integration, finite elements, basis functions, \item And finally --- to machine learning, neural networks and modern computing hardware. \end{itemize} All of this is not just a collection of disconnected topics. It is a \textbf{single language} in which modern science speaks. And if you master this language, doors will open before you that you do not even suspect exist today. \subsection{How to read this book} Although the text is written to be as easy to read as possible, without detailed mathematical proofs, some places may be formulated in too abstruse a way. That is normal --- mathematics does not come all at once. Therefore the source of this book in \LaTeX{} is attached. If something is not quite clear: \begin{enumerate} \item Find the corresponding text in the source, \item Copy it into some chat (Qwen, DeepSeek, Kimi, ChatGPT, GROK, Gemini), \item Ask it to explain this more clearly, or ask the question you do not understand. \end{enumerate} You can also use these chats to translate the whole text into English or German --- simply ask them to translate while preserving the \LaTeX{} markup, and compile the result with the command \texttt{xelatex NumMathForChemists.tex} into a PDF file. \begin{tipbox}[Tip] Do not try to read everything in one sitting. Mathematics is like music: you have to ``play'' it, that is, solve problems, try code, experiment. If some chapter seems difficult --- come back to it in a week, and you will be surprised how much easier it has become. \end{tipbox} % \subsection{Why does a chemist need mathematics?} % % Before we start examining numbers, matrices, operators and Fourier transforms, it is useful to answer the most natural question: % \textit{Why is all this needed by a chemist at all?} % The answer is very simple: because modern chemistry increasingly works not directly with substances but with \textbf{data, models and computations}. % % Let us imagine a completely ordinary experiment. We place a sample into an instrument and obtain some signal. For example, a spectrometer may return several thousand numbers. What do we want to do with these numbers? % % We want to remove noise, find peaks, determine their position and intensity, compare the obtained spectrum with known spectra, determine the concentration of a substance or reconstruct the parameters of a molecule. That is, the experiment can be represented in a very simplified form: % % $$ % \boxed{ % \text{substance} % \longrightarrow % \text{measurement} % \longrightarrow % \text{data} % \longrightarrow % \text{mathematical model} % \longrightarrow % \text{chemical conclusion} % } % $$ % % And almost every transition in this chain requires mathematics. Therefore the whole subsequent text can be perceived as a journey along one chain: % % $$ % \boxed{ % \begin{array}{c} % \text{numbers}\\ % \downarrow\\ % \text{vectors and matrices}\\ % \downarrow\\ % \text{operators}\\ % \downarrow\\ % \text{linear and nonlinear problems}\\ % \downarrow\\ % \text{statistics and uncertainty}\\ % \downarrow\\ % \text{signals and Fourier}\\ % \downarrow\\ % \text{numerical methods}\\ % \downarrow\\ % \text{compression and information extraction}\\ % \downarrow\\ % \text{machine learning}\\ % \downarrow\\ % \text{modern computational chemistry} % \end{array} % } % $$ % It is not necessary to memorise all the formulas. It is far more important to gradually learn to recognise the mathematical structure of a problem. % % % \begin{itemize} % \item If a set of numbers appears before you, ask: \textit{What kind of mathematical object is this?} % \item If many measurements appear: \textit{Can they be represented as a vector or a matrix?} % \item If there are unknown parameters: \textit{How do we build a model and find these parameters?} % \item If there is noise: \textit{How do we separate information from random error?} % \item If there is too much data: \textit{Can we find hidden structure or reduce the dimension?} % \end{itemize} % % % It is precisely this way of thinking, and not the memorisation of a large number of formulas, that is the main goal of this script. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction: Notation and the Nature of Numbers} \subsection{The hierarchy of number sets} Before we start talking about complicated things, let us agree on the language. In mathematics we group numbers into sets that have certain properties. We shall use the following standard notation: \begin{itemize} \item $\N = \{1, 2, 3, \dots\}$ --- \textbf{natural numbers}. Used for counting. \item $\Z = \{\dots, -2, -1, 0, 1, 2, \dots\}$ --- \textbf{integers}. They add zero and negative numbers, allowing us to describe debt and assets, temperatures below zero, and so on. \item $\Q$ --- \textbf{rational numbers}. Any number that can be represented as a fraction $\frac{p}{q}$, where $p \in \Z$, $q \in \Z \setminus \{0\}$. \item $\R$ --- \textbf{real numbers}. They include all rational numbers, as well as irrational numbers (for example, $\sqrt{2}, \pi, e$), which cannot be represented as an ordinary fraction. They fill the entire number line continuously. \item $\C$ --- \textbf{complex numbers}. We shall discuss them in detail below. \end{itemize} \textit{Nesting of sets:} $\N \subset \Z \subset \Q \subset \R \subset \C$. \subsection{Numbers in a computer} We are used to numbers in mathematics being infinitely precise. But when we move to \textbf{computational mathematics} and programming, the computer is forced to store numbers in memory cells of fixed size (in bits and bytes). \begin{itemize} \item \textbf{Integer types} (int8, int16, int32, int64) store exact values from $\Z$, but are limited in range (for example, int32 can store numbers from $-2^{31}$ to $2^{31}-1$). \item \textbf{Real types} (float32, float64 / double) store numbers from $\R$ in exponential format (mantissa and exponent). Because of this, rounding errors arise: the computer's $0.1 + 0.2$ is not always exactly equal to $0.3$. Understanding this is critically important for numerical methods! \end{itemize} \section{When one number is not enough: Complex numbers} Our world is so complex that not everything can be described by a single number. In physics and chemistry we often have to work with entities that require \textit{two} numbers for their description (for example, amplitude and phase of a wave, or the real and imaginary parts of a wave function). \subsection{Where do they come from?} The simplest and most historical way to get acquainted with ``number-pairs'' is to try to solve the equation: \[ x^2 + 1 = 0 \quad \Rightarrow \quad x^2 = -1 \] In the set of real numbers $\R$, the square of any number is non-negative. There is no root of $-1$. But mathematicians (and physicists after them) said: ``What if we simply \textit{invent} a number which, when multiplied by itself, gives $-1$?''. Thus the \textbf{imaginary unit} $\i$ appeared. We do not ``extract the root'' of $-1$, we \textit{define} a new entity: \[ \i^2 = -1 \] \textit{Important remark: the rule $\sqrt{a}\sqrt{b} = \sqrt{ab}$ works only for $a,b \ge 0$. Therefore we do not write $\sqrt{-1}\sqrt{-1} = \sqrt{(-1)(-1)} = 1$. We simply accept $\i^2 = -1$ as a fact.} \subsection{Algebraic form and the complex plane} Any complex number $z \in \C$ can be written as: \[ z = x + \i y \] where $x = \re(z)$ is the \textbf{real part}, and $y = \im(z)$ is the \textbf{imaginary part}. Both parts $x, y \in \R$. Geometrically, a complex number is a point on the \textbf{complex plane}. \begin{itemize} \item Along the horizontal axis (the $\re$ axis) the real part $x$ is plotted. \item Along the vertical axis (the $\im$ axis) the imaginary part $y$ is plotted. \end{itemize} A complex number can also be viewed as a \textbf{vector} going from the origin $(0,0)$ to the point $(x,y)$. \subsection{Arithmetic of complex numbers} \textbf{Addition} is intuitively clear. If we have two numbers $z_1 = x_1 + \i y_1$ and $z_2 = x_2 + \i y_2$, we simply add their real and imaginary parts separately: \[ z_1 + z_2 = (x_1 + x_2) + \i(y_1 + y_2) \] Geometrically this works as the \textbf{triangle rule} for vectors: we draw two vectors from the origin, move the start of the second vector to the end of the first, and the sum is the vector from the start of the first to the end of the second. \textbf{Multiplication} in algebraic form looks a little more complicated. We expand the brackets as in ordinary algebra, remembering that $\i^2 = -1$: \begin{align*} z_1 \cdot z_2 &= (x_1 + \i y_1)(x_2 + \i y_2) \\ &= x_1 x_2 + \i x_1 y_2 + \i y_1 x_2 + \i^2 y_1 y_2 \\ &= (x_1 x_2 - y_1 y_2) + \i(x_1 y_2 + y_1 x_2) \end{align*} But to understand the \textit{true magic} of complex multiplication, we need to move to another form of notation. \subsection{Trigonometric and exponential forms} Instead of the coordinates $(x,y)$, a point on the plane can be specified using polar coordinates: the distance from the origin (the modulus) $r$ and the angle $\varphi$ relative to the positive direction of the $\re$ axis. The connection with the algebraic form is obvious from trigonometry: \[ x = r \cos \varphi, \quad y = r \sin \varphi \] Then the complex number takes the \textbf{trigonometric form}: \[ z = r(\cos \varphi + \i \sin \varphi) \] \subsubsection{Euler's formula and the exponential form} Here the greatest formula of mathematics enters the stage --- \textbf{Euler's formula}: \[ e^{\i\varphi} = \cos \varphi + \i \sin \varphi \] Thanks to it, any complex number can be written in the incredibly convenient \textbf{exponential form}: \[ z = r e^{\i\varphi} \] \subsubsection{The geometric meaning of multiplication} Let us see what happens when two numbers are multiplied in exponential form: \[ z_1 = r_1 e^{\i\varphi_1}, \quad z_2 = r_2 e^{\i\varphi_2} \] \[ z_1 \cdot z_2 = (r_1 e^{\i\varphi_1}) \cdot (r_2 e^{\i\varphi_2}) = (r_1 r_2) e^{\i(\varphi_1 + \varphi_2)} \] \textbf{The conclusion to remember forever:} Multiplication of complex numbers is the \textbf{multiplication of their lengths (moduli)} and the \textbf{addition of their angles (arguments)}! Geometrically: multiplication by a complex number is a rotation of the vector by the angle $\varphi$ and its stretching/compression by a factor of $r$. \subsection{Complex exponentials and logarithms} \subsubsection{Properties of the exponential} Since $e^{\i\varphi}$ behaves like an ordinary exponential, all the standard rules work for it: \[ 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} \] For a complex number $z = x + \i y$ the exponential splits into real and imaginary parts: \[ e^z = e^{x + \i y} = e^x \cdot e^{\i y} = e^x (\cos y + \i \sin y) \] The modulus of this number is $e^x$, and the argument is $y$. \subsubsection{The complex logarithm} The logarithm is the operation inverse to the exponential. If $e^w = z$, then $w = \ln z$. Let $z = r e^{\i\varphi}$ and $w = u + \i v$. Then: \[ e^{u + \i v} = r e^{\i\varphi} \quad \Rightarrow \quad e^u e^{\i v} = r e^{\i\varphi} \] Equating moduli and arguments, we get: \[ e^u = r \quad \Rightarrow \quad u = \ln r \] \[ v = \varphi + 2\pi k, \quad k \in \Z \] Thus the \textbf{complex logarithm} is multivalued: \[ \ln z = \ln|z| + \i(\argg z + 2\pi k) \] The principal value of the logarithm (at $k=0$ and $\argg z \in (-\pi, \pi]$) is denoted $\operatorname{Ln} z$. The multivaluedness arises because a rotation by $2\pi$ returns us to the same point on the complex plane. \subsection{Cheat sheet: Trigonometry and hyperbolic functions} Hyperbolic functions often frighten chemistry students, but in fact they are close relatives of ordinary sines and cosines, simply ``living'' in the complex plane. \subsubsection{Definitions via the exponential} \begin{align*} \cos z &= \frac{e^{\i z} + e^{-\i z}}{2}, & \sin z &= \frac{e^{\i z} - e^{-\i z}}{2\i} \\ \cosh z &= \frac{e^{z} + e^{-z}}{2}, & \sinh z &= \frac{e^{z} - e^{-z}}{2} \end{align*} \subsubsection{Connection between trigonometric and hyperbolic functions} If we substitute $\i z$ for $z$ in the definitions, we see an amazing symmetry (Osborn's rules): \begin{align*} \cos(\i z) &= \cosh z, & \cosh(\i z) &= \cos z \\ \sin(\i z) &= \i\sinh z, & \sinh(\i z) &= \i\sin z \end{align*} \textit{Mnemonic:} When passing from trigonometry to hyperbolics the argument is multiplied by $\i$, and imaginary units may appear in front of the sine/cosine. Hyperbolic functions describe not oscillations (like the sine) but exponential growth/decay (like $e^x$), which is critically important for damped processes. \subsection{Why is this needed in real chemistry and physics?} It may seem that imaginary numbers are just a beautiful mathematical abstraction. But our world is arranged such that \textit{without} the imaginary space we could not describe real things. \subsubsection{Quantum mechanics and orbitals} The Schrödinger equation, which describes the behaviour of electrons in atoms, contains the imaginary unit $\i$ explicitly: \[ \i\hbar \frac{\partial \Psi}{\partial t} = \hat{H}\Psi \] The wave function $\Psi$ is complex-valued. For the ground state of the hydrogen atom (the 1s orbital) the wave function is real, and it can be drawn in real space. But as soon as the electron goes into an excited state (for example, a 2p orbital with magnetic quantum number $m \neq 0$), its wave function contains the phase factor $e^{\i m\varphi}$. Without complex numbers it is simply impossible to describe the angular momentum of the electron and the fine structure of spectra. \subsubsection{Nuclear Magnetic Resonance (NMR) and Fourier analysis} You will work a lot with NMR spectroscopy. When you place a sample in a magnetic field and apply a radio-frequency pulse, the nuclei begin to precess. The detector registers a signal decaying in time --- this is called the \textbf{FID} (Free Induction Decay). The signal from one type of nucleus looks like a damped oscillation: \[ S(t) = A e^{-t/T_2} \cos(\omega t) \] If the sample contains many different nuclei, the signal turns into a mess of many superimposed cosines. How can we find out which frequencies $\omega$ are hidden there? Here the \textbf{Fourier transform} comes to the rescue. But working with cosines is difficult. It is much easier to pass to complex exponentials using Euler's formula: \[ \cos(\omega t) = \frac{e^{\i\omega t} + e^{-\i\omega t}}{2} \] In NMR spectrometers quadrature detection is used, which allows measuring a complex signal directly: \[ S_{complex}(t) = A e^{-t/T_2} e^{\i\omega t} \] When we apply the fast Fourier transform (FFT) to this signal, time $t$ passes into frequency $\omega$. A long oscillating function in the time domain turns into a \textbf{narrow peak (a Lorentzian)} in the frequency domain! Each type of nucleus corresponds to its own frequency $\omega$, and in the resulting NMR spectrum we see individual peaks. All modern signal processing (from MRI in medicine to audio MP3) works precisely thanks to the magic of complex numbers and Euler exponentials. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{When one measurement is not enough: Vectors and norms} \subsection{From the plane to N-dimensional space} Well, if we have complex numbers (essentially pairs of numbers), then surely there is something more complicated? Quaternions? Octonions? Yes, indeed, mathematicians have generalised numbers from the plane to higher dimensions. But beyond complex numbers and quaternions (which, by the way, have taken root excellently in robotics and 3D graphics for describing rotations in our three-dimensional space), special ``types of numbers'' are practically not used. Instead, people simply group numbers into sets and call them \textbf{N-dimensional vectors}. This is simply N-dimensional space. $N$ can be two (then it is a plane), three (our ordinary physical space) or very, very large. We shall denote such vectors in bold, for example $\vec{v} \in \R^N$. This is a vector: \[ \vec{v} = (v_1, v_2, \dots, v_N)^T \] \textit{Here $T$ means transposition, so as to write a column in one line of text.} We are free to invent what to do with this space, but the first thing we need is to understand how to measure the ``size'' or ``length'' of such a vector. \subsection{Why do we need vectors? Examples from chemistry} At first glance, a vector is just a set of numbers. But where did the numbers come from? Suppose this is data from a chromatograph or an NMR spectrometer. Imagine we want to find out which is the ``brightest'' peak here. That is simply the maximum among these numbers, right? And now let us take a chromatogram of some complex dirt and a chromatogram of a single pure substance. We inject both into the chromatograph in the same molar volume (let us agree that the detector senses all molecules equally). \begin{itemize} \item In the case of the pure substance we get one beautiful, very sharp and tall peak. \item In the case of the dirt we get a mountain of peaks, perhaps even merged into one continuous gentle hump. \end{itemize} But here is what is interesting: the \textbf{area} (integral intensity) of both chromatograms will be the same! The peak of the pure substance will be very tall compared to the mountain, but the sum of all responses is equal to the amount of substance. From this problem three different ways naturally arise to estimate our data vector: \begin{enumerate} \item \textbf{Maximum.} We care about the maximum concentration (or maximum absorption, if the peak ``looks'' downwards). We take the largest value. \item \textbf{Sum.} We care about the total amount of substance. We add up all the values (usually in absolute value, since the baseline may not be ideal). \item \textbf{Root of the sum of squares.} And why is this needed? Imagine that we describe a molecule not by a chromatogram but by three parameters: \textit{polarity, molar mass, boiling point}. This is a vector in $\R^3$. To understand how \textit{similar} two molecules are to each other (for example, for drug design), we cannot simply add the differences of the parameters or take the maximum difference. We need precisely the \textbf{Euclidean distance} --- the shortest path in this ``chemical space''. Or another example: in physics the quadratic norm of a vector (a spectrum) is proportional to the \textbf{total energy} of that signal. \end{enumerate} \subsection{The mathematical language of norms} All these three intuitive notions are called \textbf{norms of vectors} in mathematics and are denoted by double vertical bars $\| \vec{v} \|$. \begin{itemize} \item \textbf{Maximum norm ($L_\infty$):} \[ \| \vec{v} \|_\infty = \max_{1 \le k \le N} |v_k| \] \item \textbf{Manhattan distance / sum of absolute values ($L_1$):} \[ \| \vec{v} \|_1 = \sum_{k=1}^N |v_k| \] \item \textbf{Euclidean norm / length of the vector ($L_2$):} \[ \| \vec{v} \|_2 = \sqrt{\sum_{k=1}^N |v_k|^2} \] \end{itemize} But mathematicians like order, so they combined all this into one general formula for the \textbf{$L_p$-norm}: \[ \| \vec{v} \|_p = \left( \sum_{k=1}^N |v_k|^p \right)^{1/p} \] In this formula $L_1$ is the case $p=1$. $L_2$ is $p=2$. And what will $p$ be for the maximum? Yes, exactly, $p \to \infty$. If we take the limit of this formula as $p \to \infty$, we rigorously obtain precisely the maximum element of the vector. \subsection{The continuous case: Functions as vectors} Moreover, our vector can be entirely continuous. Then it is in fact a \textbf{function} $f(x)$. For now let us just remember that a function can be viewed as an ``infinite-dimensional vector'' whose values are given at each point $x$. And for functions the norms work in exactly the same way, only the sum is replaced by an integral: \[ \| f \|_p = \left( \int |f(x)|^p \, dx \right)^{1/p} \] We shall sometimes operate with an array (a vector), sometimes with a function, and later, in the course on signal processing, we shall see that these are two sides of the same coin. \subsection{Minkowski's inequality} There is one critically important point to remember. We shall often compare $\| \vec{x} + \vec{y} \|$ with $\| \vec{x} \|$ and $\| \vec{y} \|$. It is intuitively clear that if we go from point A to point B and then from B to C, the total path will be no shorter than if we went directly from A to C. This geometric property (the length of a side of a triangle is not greater than the sum of the other two) for any $L_p$-norms is formalised in \textbf{Minkowski's inequality}: \[ \| \vec{x} + \vec{y} \|_p \le \| \vec{x} \|_p + \| \vec{y} \|_p \] \textit{How to use this?} It allows us to estimate errors. If $\vec{x}$ is the true signal and $\vec{y}$ is noise (measurement error), then Minkowski's inequality guarantees that the norm (size) of our measured signal $(\vec{x}+\vec{y})$ will not exceed the sum of the norm of the true signal and the norm of the noise. \textit{(A rigorous and beautiful proof of this inequality can be found in \href{https://en.wikipedia.org/wiki/Minkowski_inequality}{Wikipedia}, so as not to overload this script).} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Two dimensions and more: Matrices and tensors} \subsection{From vectors to matrices} Since we have a vector $\vec{a} \in \R^N$ (a one-dimensional entity, the discrete analogue of the function $f(x)$), it is logical to ask: what if we have a function of two variables $f(x,y)$? What is its discrete analogue? And for $f(x,y,z)$? Yes, two-dimensional objects in discrete space are \textbf{matrices}. And as with complex numbers, everything here is very, very well studied. And everything that has three or more dimensions is usually called \textbf{tensors}. By the way, if you ever come across old scientific literature on chemometrics or data analysis, you may be surprised by the zoo of names for tensors. They are called \textit{multilinear, multidimensional, procrustes, multimodal, three-way, multi-way arrays}. Do not be frightened: under all these beautiful words there is simply a multidimensional array of numbers. For now let us stop at two-dimensional objects: the function $f(x,y)$ and the matrix $\mathbf{A} \in \R^{N \times M}$. \subsection{Complex vectors and matrices} Why have I so far written only $\R$? Both a vector and a matrix can, of course, be complex! We can consider the spaces $\C^N$ and $\C^{N \times M}$. Here everything is simple: instead of one real number in each cell of the vector or matrix there is a complex pair. But how do we compute the norm of a complex number? In exactly the same way! The modulus of a complex number $z = x + \i y$ is computed via the conjugate number $\bar{z} = x - \i y$: \[ |z| = \sqrt{x^2 + y^2} = \sqrt{z \cdot \bar{z}} \] Accordingly, the $L_2$-norm of a complex vector $\vec{z} \in \C^N$ looks like this: \[ \| \vec{z} \|_2 = \sqrt{\sum_{k=1}^N |z_k|^2} = \sqrt{\sum_{k=1}^N z_k \bar{z}_k} \] In matrix notation this is written via the Hermitian conjugate (transpose + complex conjugate) $\vec{z}^H$: \[ \| \vec{z} \|_2 = \sqrt{\vec{z}^H \vec{z}} \] \subsection{Norms of matrices} And for matrices we can also define norms! But there is an important nuance here. Most simple norms for matrices are computed as follows: we mentally ``cut'' the matrix into rows, arrange all its $N \times M$ cells into one long column vector, and take the norm of that vector. The most popular of these norms is the \textbf{Frobenius norm} (or Euclidean norm for matrices), which is the analogue of the $L_2$-norm: \[ \| \mathbf{A} \|_F = \sqrt{\sum_{i=1}^N \sum_{j=1}^M |a_{ij}|^2} \] For it, as for the ordinary $L_2$, there are many important and beautiful properties which we shall soon consider. However, besides such ``element-wise'' norms, in linear algebra there are special \textbf{operator (induced) norms}. They answer the question: ``How strongly can a matrix stretch a vector if we multiply it by it?''. But that is a story for the next chapter, where we shall deal closely with linear operators. %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% \section{Operators: when mathematics begins to act} \subsection{What is an operator?} We have come to a very important concept which is encountered everywhere in mathematics and in chemistry. This is the \textbf{operator}. An operator is some \textit{action} on an object. But here it immediately becomes somehow confusing: an action on what? On a number? On a vector? On a function? Let us consider several examples together, and everything will fall into place. \subsection{A matrix as an operator} Suppose we have a square matrix $\mathbf{A} \in \C^{N \times N}$ (I shall henceforth almost always use $\C$, since $\R$ is merely a subset of $\C$, and everything that works for real numbers works for complex ones too). Consider the action $\mathbf{A}\vec{b}$, that is, multiplication of a matrix by a vector. The result is again a vector from $\C^N$. The matrix \textit{transforms} one vector into another. This is the simplest example of a linear operator. \subsection{Analogy with a functional operator} And now --- the most interesting thing. The analogue of a matrix operator in the world of functions is the \textbf{integral operator}: \begin{equation} (\hat{K}g)(x) = \int_a^b K(x,y)\, g(y)\, \diff y \label{IntegralOperator} \end{equation} Here $K(x,y)$ is the \textit{kernel} of the operator (the analogue of the matrix $\mathbf{A}$), $g(y)$ is the input function (the analogue of the vector $\vec{b}$), and the result is a new function of $x$. Where, interestingly, does such an abstraction occur in real life? It turns out, everywhere! For example, in quantum mechanics the Hamilton operator $\hat{H}$ acts on the wave function $\Psi$, and this is precisely an integro-differential operator. In spectroscopy the response of an instrument to an input signal is often described by a convolution integral --- that too is an operator. \subsection{Back to school: a system of equations} But let us start with the simplest thing. At school we already solved systems of equations, for example: \[ \begin{cases} 2x - 3y = -4 \\ 3x + 2y = 9 \end{cases} \] You are probably thinking: ``So what? We could solve that by substitution or by elimination''. Yes, we could. But let us rewrite it in matrix form $\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}} \] You may reasonably ask me: ``Why complicate it like that?''. And I answer: no, we are not complicating --- we are \textit{putting in order} our knowledge. We are isolating the \textit{structure} of the problem. And now we can say: if we have a \textbf{linear operator} $\mathbf{A}$, then it acts on the vector $\vec{x}$ and the result is the vector $\vec{b}$. \subsection{Reminder: basic operations} Before going further, let us fix explicitly the operations we shall constantly use. \textbf{A vector as a column matrix.} A vector $\vec{x} \in \C^N$ is in fact a matrix of size $N \times 1$. Therefore all the rules of matrix multiplication apply automatically to vectors too. \textbf{Transposition.} The operation $\mathbf{A}^T$ swaps rows and columns: $(\mathbf{A}^T)_{ij} = A_{ji}$. For complex matrices the \textbf{Hermitian conjugate} $\mathbf{A}^H = \overline{\mathbf{A}^T}$ (transposition + complex conjugation of all elements) is more often used. \textbf{Scalar product.} For two vectors $\vec{x}, \vec{y} \in \C^N$ the scalar product is defined as: \[ \langle \vec{x}, \vec{y} \rangle = \vec{x}^H \vec{y} = \sum_{k=1}^N \bar{x}_k y_k \] For real vectors this is simply $\vec{x}^T \vec{y} = \sum x_k y_k$. \textbf{The most important property:} the scalar product of a vector with itself equals the \textbf{square of its Euclidean norm}: \[ \langle \vec{x}, \vec{x} \rangle = \|\vec{x}\|_2^2 = \sum_{k=1}^N |x_k|^2 \] This is the bridge between algebra and geometry: the length of a vector is the root of its ``scalar square''. \textbf{Multiplication of a matrix by a vector.} If $\mathbf{A} \in \C^{M \times N}$ and $\vec{x} \in \C^N$, then the result $\vec{y} = \mathbf{A}\vec{x} \in \C^M$ is computed by the rule: \[ y_i = \sum_{j=1}^N A_{ij} x_j \] That is, the $i$-th component of the result is the scalar product of the $i$-th row of the matrix with the vector $\vec{x}$. \textbf{Multiplication of matrices.} If $\mathbf{A} \in \C^{M \times K}$ and $\mathbf{B} \in \C^{K \times N}$, then the product $\mathbf{C} = \mathbf{A}\mathbf{B} \in \C^{M \times N}$: \[ C_{ij} = \sum_{k=1}^K A_{ik} B_{kj} \] Important: matrix multiplication is, generally speaking, \textbf{not commutative} --- $\mathbf{A}\mathbf{B} \neq \mathbf{B}\mathbf{A}$. This is one of the most important differences from the multiplication of ordinary numbers! All these operations --- scalar product, matrix-vector multiplication, matrix-matrix multiplication --- have their complete analogues in the world of functions and integral operators. Only everywhere where we had a sum $\sum$, an integral $\int$ appears, and transposition is replaced by a more complicated conjugation operator. But that is already a matter of technique. \subsection{The four main problems of linear algebra} In fact, as long as our operators are linear (that is, they can be written in the form (\ref{IntegralOperator}) or $\mathbf{A}\vec{x}$), the whole world of possible problems revolves around a few typical questions. Let us list them. \textbf{Problem 1. Compute a scalar product.} Given two vectors --- find the number $\langle \vec{x}, \vec{y} \rangle$. This is the basic operation from which, like bricks, everything else is built. \textbf{Problem 2. Compute the action of an operator.} Given $\mathbf{A}$ and $\vec{x}$ --- find $\mathbf{A}\vec{x}$. Or given $\mathbf{A}$ and $\mathbf{B}$ --- find $\mathbf{A}\mathbf{B}$. This is a direct computation. \textbf{Problem 3. Solve a linear system of equations.} Given $\mathbf{A}$ and $\vec{b}$ --- find $\vec{x}$ such that $\mathbf{A}\vec{x} = \vec{b}$. This is the \textit{direct} problem: the operator and the result are known, we must find what it acted on. \textbf{Problem 4. Minimise the residual (least squares method).} But what if the system $\mathbf{A}\vec{x} = \vec{b}$ has no exact solution? (For example, there are more equations than unknowns --- an overdetermined system, which is typical for processing experimental data.) Then we look for $\vec{x}$ that \textit{minimises} the residual: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_p \] Almost always $p = 2$ is used as the norm --- this is the famous \textbf{least squares method (LS)}. It has a deep geometric interpretation: we project the vector $\vec{b}$ onto the subspace spanned by the columns of the matrix $\mathbf{A}$. \subsection{And what if the right-hand side is zero?} And now --- a cunning turn. Imagine that $\vec{b} = \vec{0}$. Then the problem $\min_{\vec{x}} \|\mathbf{A}\vec{x}\|_2$ has the trivial solution $\vec{x} = \vec{0}$. But that is uninteresting! Let us impose an additional condition: we look for a \textbf{nonzero} solution, for example with the constraint $\|\vec{x}\|_2 = 1$. That is, we look for the direction in which the matrix $\mathbf{A}$ ``compresses'' space most strongly (or, conversely, stretches it most strongly --- that is already a question of sign). And here \textbf{eigenvalues and eigenvectors} enter the stage. We look for such special vectors $\vec{v}$ and numbers $\lambda$ for which the action of the matrix reduces simply to stretching: \[ \mathbf{A}\vec{v} = \lambda \vec{v} \] Geometrically: an eigenvector is a direction which the matrix \textit{does not rotate}, but only stretches or compresses by a factor of $\lambda$. The problem of minimising $\|\mathbf{A}\vec{x}\|_2$ subject to $\|\vec{x}\|_2 = 1$ is solved precisely via the eigenvalues of $\mathbf{A}^H \mathbf{A}$ or the singular values of $\mathbf{A}$: the minimum is attained on the eigenvector of the matrix $\mathbf{A}^H\mathbf{A}$ corresponding to the \textit{smallest in modulus} eigenvalue. And the maximum --- on the vector corresponding to the largest. \subsection{And if the operator itself is unknown?} You may ask: ``And what if the operator is unknown to us, and we know only the input functions and the results?'' Here the answer is somewhat ambiguous. Think: in matrix form we have only $N$ input parameters (the vector $\vec{x}$), and we want to find $N^2$ parameters of the matrix $\mathbf{A}$. Nature is rarely arranged such that a small number of inputs determines a larger number of parameters well. This is a classical \textbf{underdetermined} problem. But here too there is a solution! If we say that the same matrix $\mathbf{A}$ acts on several different vectors $\vec{x}_1, \dots, \vec{x}_K$ with known results $\vec{b}_1, \dots, \vec{b}_K$, then we can formulate the problem as: \begin{equation} \forall k = 1, \dots, K: \quad \mathbf{A}\vec{x}_k = \vec{b}_k, \quad \text{where } \mathbf{A} \text{ is unknown.} \label{MatApprox} \end{equation} If we collect the vectors into matrices $\mathbf{X} = (\vec{x}_1, \dots, \vec{x}_K)$ and $\mathbf{B} = (\vec{b}_1, \dots, \vec{b}_K)$, the problem is rewritten as: \[ \mathbf{A}\mathbf{X} \approx \mathbf{B} \] And this, attention, is \textbf{the same minimisation problem}, only now we minimise not over $\vec{x}$ but over $\mathbf{A}$: \[ \min_{\mathbf{A}} \|\mathbf{A}\mathbf{X} - \mathbf{B}\|_F \] (Here $\|\cdot\|_F$ is the Frobenius norm for matrices, which we discussed in the previous chapter.) Let us transpose for convenience: $\mathbf{X}^T \mathbf{A}^T \approx \mathbf{B}^T$. And we have again obtained a problem of the form ``minimise $\|\mathbf{M}\vec{z} - \vec{c}\|_2$'', only for each column of $\mathbf{A}^T$ separately. This is classical \textbf{linear regression} --- the basis of chemometrics, QSAR, instrument calibration and much else. \subsection{Singular value decomposition: a bridge into big worlds} And here one of the most beautiful constructions in all of linear algebra enters the stage --- the \textbf{singular value decomposition (SVD)}. Any matrix $\mathbf{A} \in \C^{M \times N}$ can be represented as: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] where $\mathbf{U} \in \C^{M \times M}$ and $\mathbf{V} \in \C^{N \times N}$ are unitary matrices (their columns are orthonormal vectors), and $\mathbf{\Sigma}$ is a diagonal matrix with non-negative numbers $\sigma_1 \ge \sigma_2 \ge \dots \ge 0$ on the diagonal. These numbers are called \textbf{singular values}. The geometric meaning of SVD is astonishing: any matrix is simply a rotation ($\mathbf{V}^H$), then a stretching along the axes ($\mathbf{\Sigma}$), then another rotation ($\mathbf{U}$). That is all! Matrices can do nothing more. SVD is a universal tool. Through it one solves: \begin{itemize} \item least squares problems, \item finding the pseudoinverse matrix, \item data compression and principal component extraction (PCA --- Principal Component Analysis, the basis of chemometrics!), \item regularisation of ill-conditioned problems. \end{itemize} There is an important point here: unitary matrices have one very important property, namely that always $\mathbf{V}^H \mathbf{V} = \mathbf{I}$, where $\mathbf{I}$ is the identity matrix, that is, one with ones on the main diagonal and everything else filled with zeros. \subsection{A bridge into functional analysis: Fredholm theory} And now --- the most interesting thing. Do you remember the integral operator (\ref{MatApprox})? Well, there is an analogue of SVD for it too! It is called the \textbf{singular value decomposition of a compact operator} or, in a more classical formulation, \textbf{Fredholm theory}. The essence is as follows: the kernel of the integral operator $K(x,y)$ can be expanded in a series in ``singular functions'' $u_k(x)$ and $v_k(y)$: \[ K(x,y) = \sum_{k=1}^\infty \sigma_k \, u_k(x) \overline{v_k(y)} \] where $\sigma_k$ are the singular values (decreasing to zero), and $u_k, v_k$ are orthonormal systems of functions. This is exactly the analogue of SVD for matrices, only in an infinite-dimensional space! And all the ideas we understood on matrices carry over here almost word for word. But let us not overload our minds with Fredholm and functional analysis right now. The good news is that \textbf{most of what we need in chemistry and signal processing can be done at the matrix level}. And where it no longer works (for example, in quantum mechanics or in the rigorous theory of integral equations), we simply remember that ``somewhere out there is Fredholm'', and if necessary ask AI for details --- it will gladly tell us about Fredholm theory, about Hilbert--Schmidt kernels and about the spectra of compact operators. The main thing is to understand the \textit{structure}. And the structure is the same everywhere: operator, space, action, expansion in ``basis directions''. %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Conditioning: when a matrix ``wants'' to deceive us} \subsection{When a solution exists, but seems not to} So, let us have a square matrix $\mathbf{A}$ and a linear system $\mathbf{A}\vec{x} = \vec{b}$. What can we say about it? We remember from the school course that sometimes a system of equations is such that, when we start substituting, either there are no solutions, or we end up with two or more identical equations --- and then there are infinitely many solutions. In chemistry, and indeed in any industrial problem, the initial data for matrix equations often come from measurements with instruments. And here we may be unlucky in three ways at once: \begin{enumerate} \item Our system will be degenerate (no unique solution). \item Our system will have many solutions. \item \textbf{The most insidious case:} solving in machine arithmetic, we find ourselves \textit{very close} to the bad case, but do not notice it. And we get an answer that looks reasonable, but in fact is complete nonsense. \end{enumerate} When can this happen? Let us sort it out with a vivid example. \subsection{An example from chromatography: the trap of the large and the small} Suppose we have two chromatograms in which two substances were recorded against the background of a solvent. The solvent peak came out very broad and ``crept'' onto the peaks of the substances sought, while the responses of the substances sought turned out to be very weak. In fact, in the peak of the substances sought we have: \begin{itemize} \item First chromatogram: $(a + b_1)$, where $a$ is a huge response from the solvent, and $b_1$ is a very small number from the substance sought. \item Second chromatogram: $(a + b_2)$, but since the temperature became a little higher, this measurement is slightly inaccurate, say, it increased by some value $\varepsilon$ \textit{relative to the solvent peak}, that is $(1+\varepsilon)(a+b_2)$. \end{itemize} If we want to subtract the second from the first to compute $b_1 - b_2$, then instead we get: \[ (a+b_1) - (1+\varepsilon)(a+b_2) = b_1 - b_2 - \varepsilon(a+b_2) \approx b_1 - b_2 - \varepsilon a \] That is, if $\varepsilon a$ is comparable in order of magnitude to $b_1 - b_2$, we can get not just a large error, but also the \textbf{opposite sign}! This is the essence of poor conditioning: when one number contains simultaneously a very large and a very small component, and we try to extract the small one by subtraction. \begin{warningbox}[Rule for subtracting close numbers] Never subtract two numbers of similar magnitude if you need relative accuracy. Loss of significant digits is not a bug of the computer, it is a fundamental property of arithmetic. \end{warningbox} \subsection{How SVD reveals the mechanism of the catastrophe} The same thing happens when solving systems of equations. And we have a nice way to \textit{estimate in advance} what to do. Suppose we have a system $\mathbf{A}\vec{x} = \vec{b}$, where we seek $\vec{x}$ and everything else is given. Recall the singular value decomposition $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$. Then the solution can be rewritten in several stages: \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*} Here we use the property that solving a system with unitary matrices $\mathbf{U}, \mathbf{V}$ is very simple --- it suffices to multiply by $\mathbf{U}^H$ or $\mathbf{V}^H$, since $\mathbf{U}^H \mathbf{U} = \mathbf{I}$. And solving a system with a diagonal matrix is altogether trivial. Let us draw how this looks for a $3 \times 3$ matrix: \[ \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} \] Do you see? Each equation is simply a division of one component by the corresponding singular value. And here is the most important thing. Imagine that we have computed the singular values of the original matrix and noticed that the largest singular value $\sigma_1$ is very, very much larger than the smallest $\sigma_N$. Then at the step $\vec{t}_2 = \mathbf{\Sigma}^{-1} \vec{t}_1$ the following happens: \begin{itemize} \item In the vector $\vec{t}_1$ the measurement error is distributed more or less uniformly over all components. \item After multiplication by $\mathbf{\Sigma}^{-1}$ the component corresponding to the \textit{small} $\sigma_N$ grows by a factor of $\sigma_1/\sigma_N$! \item And if $\sigma_N$ is altogether zero? Then we divide by zero --- and the solution does not exist. \end{itemize} It is precisely because of this that people decided to denote the ratio $\sigma_1 / \sigma_N$ as the \textbf{condition number} of the matrix: \[ \cond(\mathbf{A}) = \frac{\sigma_{\max}}{\sigma_{\min}} \] In fact, this is a characteristic of \textit{how much we can spoil the solution with such a matrix}. If $\cond(\mathbf{A}) \approx 10^k$, then when solving the system we lose approximately $k$ significant digits of accuracy. For double precision (16 significant digits) this means that at $\cond(\mathbf{A}) \approx 10^{16}$ the answer will consist entirely of noise. \begin{tipbox}[Important remark] The condition number is a characteristic of the \textit{matrix}, not of the algorithm. No algorithm, however ingenious, can solve an ill-conditioned system more accurately than $\cond(\mathbf{A})$ permits. This is a fundamental limitation of the problem, not of our computational methods. At the same time, conditioning \textit{does not affect} the computation of eigenvectors of symmetric/Hermitian matrices or the finding of singular vectors --- these problems are, as a rule, much more stable. \end{tipbox} \subsection{Where does all this occur in chemistry and beyond?} We have had a lot of dry theory, and it would be good to recall where all this is applied in chemistry. \textbf{Quantum chemistry: the Schrödinger equation.} Yes, the Schrödinger wave equation and all its approximations --- the Hartree--Fock equations, Density Functional Theory (DFT) --- all these are eigenvalue problems: finding vectors $\vec{x}$ and values $\lambda$ satisfying the equation $\mathbf{A}\vec{x} = \lambda \vec{x}$. The matrices here --- Fock matrices or Hamiltonian matrices --- often have dimensions of thousands and tens of thousands, and their conditioning directly determines whether we can obtain a physically meaningful answer at all. \textbf{Mass spectrometry and spectroscopy.} Finding a recorded spectrum from a database of responses in a mass spectrometer --- these are algorithms based on solving linear systems. If you have a mixture of substances and want to understand what it consists of, you solve the system $\mathbf{A}\vec{x} = \vec{b}$, where $\mathbf{A}$ is the matrix of reference spectra, $\vec{b}$ is the measured spectrum of the mixture, and $\vec{x}$ are the concentrations of the components. And if the spectra of the components are similar --- the matrix is ill-conditioned, and the concentrations will ``jump''. \textbf{Cryo-electron microscopy.} Reconstruction of the 3D structure of proteins from thousands of 2D projections is a gigantic sparse system of equations. Conditioning there is a key factor determining the resolution of the resulting structure. \textbf{NMR spectroscopy.} Deconvolution (inverse convolution) is a classical ill-conditioned problem. That is precisely why NMR uses Tikhonov regularisation and the Fourier transform --- they turn an ill-conditioned problem into a diagonal one. \textbf{Calibration of analytical instruments.} PLS regression (Partial Least Squares) is, in essence, SVD with regularisation. The conditioning of the matrix of spectra directly affects the accuracy of concentration prediction. \textbf{Internet search (PageRank).} Even the very first Google search was arranged so that all previously found documents were classified into a huge matrix, and its singular value decomposition was computed (or, equivalently, the principal eigenvector of the transition matrix was found), which helped to find the desired document instantly by a matching combination of words. In modern AI, SVD is used very, very often --- from compressing neural network weights to dimensionality reduction in recommender systems. \textbf{Signal and image processing.} Deconvolution of blurred images in microscopy, noise suppression in ECG/EEG, compression of audio (MP3) and video --- all these are problems in which conditioning plays a key role. Once upon a time these three problems --- solving linear systems, finding eigenvalues and SVD --- were a stumbling block in solving large industrial problems. There were even firms that developed \textit{only} such solvers --- and nothing else. True, that was back in the 1990s. \subsection{How much does it cost? Computational complexity} Roughly speaking, if we have a square $N \times N$ matrix and we either solve a linear system with it, or find its eigenvalues or singular values and vectors, then we must spend approximately $\mathcal{O}(N^3)$ arithmetic operations, if we do not specially use the structure or properties of this matrix. This is the ``price'' of the full problem. But if the matrix has special properties, the price can be greatly reduced --- sometimes to $\mathcal{O}(N)$ or $\mathcal{O}(N \log N)$. And that is discussed below. % === LINK BETWEEN SVD AND EIGENVALUES === \subsection{The connection between singular value decomposition and eigenvalues} And now --- one of the most beautiful facts in all of linear algebra. Suppose we have the singular value decomposition: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] Let us see what happens if we multiply $\mathbf{A}$ by its Hermitian conjugate $\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*} Here we used that $\mathbf{V}^H \mathbf{V} = \mathbf{I}$ (unitarity of $\mathbf{V}$), and $\mathbf{\Sigma}^T = \mathbf{\Sigma}$ (diagonal matrix). Similarly: \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{What does this mean?} \begin{itemize} \item The singular vectors of the matrix $\mathbf{A}$ are \textit{exactly} the eigenvectors of the matrices $\mathbf{A}\mathbf{A}^H$ and $\mathbf{A}^H\mathbf{A}$. \item The singular values $\sigma_k$ of the matrix $\mathbf{A}$ are the square roots of the eigenvalues $\lambda_k$ of the matrices $\mathbf{A}\mathbf{A}^H$ or $\mathbf{A}^H\mathbf{A}$: \[ \sigma_k = \sqrt{\lambda_k} \] \end{itemize} This is a deep connection between two seemingly different problems: singular value decomposition and eigenvalue finding. \begin{warningbox}[Caution: the square of the condition number!] But there is one insidious detail. The condition number of the matrices $\mathbf{A}\mathbf{A}^H$ and $\mathbf{A}^H\mathbf{A}$ is the \textbf{square} of the condition number of the original matrix: \[ \cond(\mathbf{A}\mathbf{A}^H) = \cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2 \] If $\cond(\mathbf{A}) = 10^8$, then $\cond(\mathbf{A}^H\mathbf{A}) = 10^{16}$ --- and we will lose all 16 significant digits of accuracy! Therefore, passing from the singular value decomposition problem to the eigenvalue problem for $\mathbf{A}^H\mathbf{A}$ must be done \textbf{with extreme caution}. In modern libraries (LAPACK) SVD is computed directly, without explicitly forming $\mathbf{A}^H\mathbf{A}$ --- precisely in order to avoid this catastrophe. \end{warningbox} %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Classification of matrices: who is who} Now that it is clear that practically everything is based on these ``simple'' mathematical problems, let us systematise a little our knowledge of such solutions. After all, we always have a matrix, and much depends on the properties of this matrix. \subsection{By shape} \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Type} & \textbf{Size} & \textbf{Description} \\ \midrule Square & $N \times N$ & Equal number of rows and columns. Only for such matrices are the determinant, eigenvalues and inverse matrix defined. \\ ``Tall'' & $M \times N, \; M > N$ & More equations than unknowns. Usually has no exact solution --- we solve via least squares. \\ ``Wide'' & $M \times N, \; M < N$ & More unknowns than equations. Infinitely many solutions --- we seek the one of minimal norm. \\ \bottomrule \end{tabularx} \subsection{By symmetry and structure of elements} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Type} & \textbf{Definition and properties} & \textbf{Complexity} \\ \midrule Symmetric real & $\mathbf{A} = \mathbf{A}^T$, elements $\in \R$. Eigenvalues \textbf{real}, eigenvectors --- orthogonal. & $\mathcal{O}(N^3/3)$ \\ Hermitian complex & $\mathbf{A} = \mathbf{A}^H$, elements $\in \C$. Eigenvalues \textbf{real}, eigenvectors --- orthonormal. & $\mathcal{O}(4N^3/3)$ \\ Skew-symmetric & $\mathbf{A} = -\mathbf{A}^T$. Eigenvalues --- purely imaginary or zero. & $\mathcal{O}(N^3/3)$ \\ Orthogonal / Unitary & $\mathbf{A}^T\mathbf{A} = \mathbf{I}$ / $\mathbf{A}^H\mathbf{A} = \mathbf{I}$. Preserves lengths and angles. Eigenvalues equal 1 in modulus. & $\mathcal{O}(N^2)$ for multiplication \\ Identity & $\mathbf{A} = \mathbf{I}$. Diagonal of ones, the rest --- zeros. & $\mathcal{O}(N)$ \\ Diagonal & $A_{ij} = 0$ for $i \neq j$. Stored as a vector of $N$ numbers. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} \subsection{By positive definiteness} This is a subclass of symmetric/Hermitian matrices, and it is critically important for optimisation and statistics. \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Type} & \textbf{Definition and properties} \\ \midrule Positive definite ($\mathbf{A} \succ 0$) & $\forall \vec{x} \neq 0: \; \vec{x}^H \mathbf{A} \vec{x} > 0$. All eigenvalues strictly positive. Conditioning --- the ratio of the largest and smallest eigenvalues. Solved by the Cholesky method in $\mathcal{O}(N^3/3)$. \\ Positive semi-definite ($\mathbf{A} \succeq 0$) & $\forall \vec{x}: \; \vec{x}^H \mathbf{A} \vec{x} \ge 0$. All eigenvalues $\ge 0$. May be degenerate. \\ \bottomrule \end{tabularx} \subsection{By the structure of the arrangement of nonzero elements} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Type} & \textbf{Definition and properties} & \textbf{Complexity} \\ \midrule Banded & Nonzero elements only near the main diagonal: $A_{ij} = 0$ for $|i-j| > k$. Stored as $N \times (2k+1)$. & $\mathcal{O}(N k^2)$ \\ Toeplitz & $A_{ij}$ depends only on $i-j$. Each diagonal is a constant. Arises in problems with stationary processes. & $\mathcal{O}(N^2)$, via FFT --- $\mathcal{O}(N \log N)$ \\ Circulant & Toeplitz + periodicity: each row is a cyclic shift of the previous one. Diagonalised by the Fourier transform. & $\mathcal{O}(N \log N)$ via FFT \\ Block & Consists of blocks, each of which is a matrix. Allows recursive algorithms. & Depends on the block structure \\ Sparse & Most elements are zeros. Stored in special formats (CSR, CSC). The basis of large computations. & $\mathcal{O}(\text{nnz})$, where nnz is the number of nonzeros \\ \bottomrule \end{tabularx} \subsection{Degeneracy} A \textbf{degenerate (singular) matrix} is a matrix whose determinant is zero, or, equivalently, which has at least one singular value equal to zero. For such a matrix: \begin{itemize} \item there is no inverse matrix, \item the system $\mathbf{A}\vec{x} = \vec{b}$ either has no solutions or has infinitely many, \item $\rank(\mathbf{A}) < N$. \end{itemize} In practice matrices are almost never \textit{exactly} degenerate --- but they can be \textit{almost} degenerate, that is, with a very small $\sigma_{\min}$. And this is a much more insidious case, because formally a solution exists, but it is entirely determined by noise in the data. % === SUPPLEMENT TO THE CLASSIFICATION OF MATRICES === \subsection{Additional important types of matrices} \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Type} & \textbf{Definition and properties} & \textbf{Complexity} \\ \midrule Upper triangular & $A_{ij} = 0$ for $i > j$. All nonzero elements above or on the main diagonal. Solving the system --- back substitution. & $\mathcal{O}(N^2)$ \\ Lower triangular & $A_{ij} = 0$ for $i < j$. All nonzero elements below or on the main diagonal. Solving the system --- forward substitution. & $\mathcal{O}(N^2)$ \\ Permutation matrix & In each row and each column exactly one unit, the rest --- zeros. Multiplication by such a matrix is a permutation of rows/columns. $\mathbf{P}^T = \mathbf{P}^{-1}$. & $\mathcal{O}(N)$ \\ \bottomrule \end{tabularx} Permutation matrices are not just an abstraction. They are critically important for the numerical stability of algorithms (for example, in LU decomposition with pivoting), and we shall meet them again. \subsection{Summary table: what to solve and how} \begin{longtable}{l c c c} \toprule \textbf{Problem} & \textbf{General case} & \textbf{Symmetric} & \textbf{Sparse} \\ \midrule Linear system $\mathbf{A}\vec{x} = \vec{b}$ & LU, $\mathcal{O}(N^3)$ & Cholesky, $\mathcal{O}(N^3/3)$ & Iterative, $\mathcal{O}(\text{nnz} \cdot k)$ \\ Least squares $\min \|\mathbf{A}\vec{x} - \vec{b}\|_2$ & QR, $\mathcal{O}(MN^2)$ & --- & LSQR, iterative \\ Eigenvalues & $\mathcal{O}(N^3)$ & $\mathcal{O}(N^3/3)$, all real & Lanczos, $\mathcal{O}(\text{nnz} \cdot k)$ \\ SVD & $\mathcal{O}(MN^2)$ & Via eigenvalues of $\mathbf{A}^T\mathbf{A}$ & Truncated SVD, iterative \\ \bottomrule \end{longtable} Here $k$ is the number of iterations, which depends on the desired accuracy and the conditioning. \begin{tipbox}[Main conclusion] Before solving a problem, \textit{look at the matrix}. Its properties --- symmetry, sparsity, banded structure --- can reduce the complexity from cubic to linear. And its conditioning will tell you whether the answer obtained is worth trusting at all. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{How linear algebra problems are solved: from theory to practice} \subsection{There is no single solution for all cases} So how are these linear algebra problems solved? Surely there is one, or perhaps a couple of the simplest and most reliable solutions? Unfortunately, even now there is no single solution for all cases of life, and none is in sight. The author of this script once participated in writing more than a hundred different such algorithms. Yes, there are very many such algorithms, and one must at least briefly, even without ever implementing it, understand when and what can be used. \subsection{Classification of solution methods} Solution methods can be classified roughly as follows: \subsubsection{Direct decomposition methods} \textbf{LU decomposition:} $\mathbf{A} = \mathbf{L}\mathbf{U}$, where $\mathbf{L}$ is lower triangular and $\mathbf{U}$ is upper triangular. With such a factorisation, solving the system $\mathbf{A}\vec{x} = \vec{b}$ reduces to two simple substitutions: \begin{enumerate} \item Solve $\mathbf{L}\vec{y} = \vec{b}$ (forward substitution, $\mathcal{O}(N^2)$) \item Solve $\mathbf{U}\vec{x} = \vec{y}$ (back substitution, $\mathcal{O}(N^2)$) \end{enumerate} But in general $\cond(\mathbf{L}) \cdot \cond(\mathbf{U}) \ge \cond(\mathbf{A})$, and if no permutations are made, the product of condition numbers may become substantially larger than $\cond(\mathbf{A})$. That is, we can spoil our original problem! Therefore in practice one almost always uses \textbf{LUP decomposition}: $\mathbf{P}\mathbf{A} = \mathbf{L}\mathbf{U}$, where $\mathbf{P}$ is a permutation matrix. The permutations are chosen so that the diagonal of $\mathbf{U}$ contains the largest possible elements (pivoting). This guarantees numerical stability. \textbf{For special matrices:} \begin{itemize} \item If $\mathbf{A}$ is banded, then $\mathbf{L}$ and $\mathbf{U}$ do not go beyond the band. Complexity $\mathcal{O}(N k^2)$, where $k$ is the band width. \item If $\mathbf{A}$ is sparse, then $\mathbf{L}$ and $\mathbf{U}$ may contain substantially more nonzero elements (fill-in), and this can be a big problem. \item For a symmetric matrix --- $\mathbf{A} = \mathbf{L}\mathbf{D}\mathbf{L}^H$ or $\mathbf{P}\mathbf{A}\mathbf{P}^H = \mathbf{L}\mathbf{D}\mathbf{L}^H$. \item For a positive definite one --- $\mathbf{A} = \mathbf{L}\mathbf{L}^H$, the so-called \textbf{Cholesky method}. But the problem of growth of the condition number exists even for positive definite matrices. \end{itemize} \textbf{QR decomposition:} $\mathbf{A} = \mathbf{Q}\mathbf{R}$, where $\mathbf{Q}$ is unitary ($\mathbf{Q}^H\mathbf{Q} = \mathbf{I}$) and $\mathbf{R}$ is upper triangular. This method \textbf{guarantees no increase in the condition number} of the solution, since unitary transformations preserve the lengths of vectors. It is good for dense matrices without structure. But if the original matrix is sparse, then $\mathbf{Q}$ will almost always be dense (only for very special cases and methods), which strongly limits the applicability of this method to solving huge problems, when the matrix itself fits in memory but its dense version $N^2$ does not. \textbf{Schur decomposition:} $\mathbf{A} = \mathbf{Q}\mathbf{G}\mathbf{Q}^H$, where $\mathbf{G}$ is upper triangular (for real matrices --- upper quasi-triangular with $2 \times 2$ blocks on the diagonal). It is used for finding eigenvalues and eigenvectors, which are found by solving the eigenvalue problem for the already triangular matrix $\mathbf{G}$ (and for a triangular matrix the eigenvalues are simply the diagonal elements!). \subsubsection{Iterative methods} These are methods for solving a linear system or an eigenvalue problem when we can efficiently multiply by the matrix because of its structure, and we build the solution iteratively. There are many methods here, but it is important to note: \textbf{the convergence of iterative methods depends on the condition number}. The larger it is, the more slowly they converge. (There are exceptions: if the matrix has only a few very large and very small singular values, then convergence can be substantially faster.) All these methods can also be divided into subclasses: \begin{enumerate} \item \textbf{Special methods for positive definite matrices} with small additional memory costs and a reasonable number of iterations. The best known is the \textbf{conjugate gradient method (CG)}. \item \textbf{Reduction of a nonsymmetric problem to a symmetric one:} instead of solving $\mathbf{A}\vec{x} = \vec{b}$, solve $\mathbf{A}^H\mathbf{A}\vec{x} = \mathbf{A}^H\vec{b}$. If the condition number $\cond(\mathbf{A}^H\mathbf{A})$ is not so enormous compared to $\cond(\mathbf{A})$, then this is reasonable. But remember: $\cond(\mathbf{A}^H\mathbf{A}) = \cond(\mathbf{A})^2$, so this is a double-edged sword. \item \textbf{Preconditioners:} special matrices $\mathbf{P} \approx \mathbf{A}^{-1}$ such that instead of solving $\mathbf{A}\vec{x} = \vec{b}$ we solve $\mathbf{P}\mathbf{A}\vec{x} = \mathbf{P}\vec{b}$, assuming that $\cond(\mathbf{P}\mathbf{A}) \ll \cond(\mathbf{A})$, which greatly accelerates the solution by such iterative methods. By the way, often for large sparse matrices preconditioners are constructed as an \textbf{incomplete LU decomposition}: that is, for the original matrix one constructs $\mathbf{A} \approx \mathbf{L}\mathbf{U}$, but tries to ``throw out'' new nonzero elements when constructing $\mathbf{L}\mathbf{U}$. In that case, applying this matrix as something similar to the original is still possible, so by solving with it we precondition the original problem and improve convergence. \end{enumerate} \subsection{Two worlds: dense and sparse problems} In fact, most solution algorithms can be divided into two classes: \begin{enumerate} \item When we are prepared to spend $\mathcal{O}(N^3)$ arithmetic operations and we have $\mathcal{O}(N^2)$ memory. \item When we need to squeeze, and the matrix is so huge and has some special structure that we can multiply by it, but taking it in full form is very difficult. \end{enumerate} In the first case --- these are \textbf{decomposition methods} (LU, QR, Cholesky). In the second --- these are \textbf{iterative methods}. Moreover, because of the specific structure of modern processors and memory, iterative methods have somewhat worse performance than full decomposition methods. Therefore applying iterative methods to very small matrices is practically never justified. \subsection{Examples from quantum chemistry: DFT} \textbf{Example 1: A small system.} A small system of a few electrons by the DFT method, and the size of the Hamiltonian is about $1000 \times 1000$. We fit in memory (that is only $\sim 8$ MB for double precision). We can find all eigenvalues and the set of eigenvectors we are interested in, and for this the obvious choice is a direct method (for example, Schur decomposition or the QR algorithm). \textbf{Example 2: A large chemical structure.} A rather large chemical structure by the DFT method, in which there are about a thousand electron pairs in the outer orbitals. For this we have constructed the Hamiltonian of the Schrödinger equation, and we must find one eigenvector for each electron pair. And the Hamiltonian is given such that we have taken about a hundred basis functions for each orbital, that is, the dimension of this Hamiltonian is about $100\,000 \times 100\,000$. This seems not so much, but the matrix of such a Hamiltonian alone already occupies about \textbf{80 GB} in RAM, and the solution itself will take almost a gigabyte more. Here an iterative solution is much more obvious (for example, the Lanczos method or Davidson), which finds only the few eigenvectors needed, without working with the full matrix. \subsection{Do not reinvent the wheel: libraries} Nowadays in many programming languages there are well and for years polished libraries for solving such systems of equations, finding singular values and vectors, and eigenvectors and eigenvalues. It is enough to look in the documentation for \textbf{NumPy} for Python, and in \textbf{LAPACK/BLAS} for C/C++ --- and everything will immediately become clear. Moreover, the code is very well optimised for modern computers and rather complicated. For example, the modern version of SVD in LAPACK contains about \textbf{half a million lines of code}, and repeating it or even doing better is really an almost impossible task. \begin{tipbox}[Practical advice] Never write your own solvers for linear algebra, unless it is a study task. Use: \begin{itemize} \item \textbf{Python:} NumPy (\texttt{numpy.linalg}), SciPy (\texttt{scipy.linalg}) \item \textbf{C/C++:} LAPACK, BLAS, Eigen \item \textbf{Fortran:} LAPACK (this is its native element) \item \textbf{MATLAB:} built-in functions (\texttt{eig}, \texttt{svd}, \texttt{lu}, \texttt{qr}) \end{itemize} These libraries have gone through decades of optimisation and testing. They know about the processor cache, about SIMD instructions, about multithreading --- all that you cannot repeat by hand. \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Nonlinear operators: when the world ceases to be straight} \subsection{What is a nonlinear operator?} Until now we have spoken of linear operators --- those that can be written as $\mathbf{A}\vec{x}$ or $\int K(x,y)g(y)\diff y$. But the real world is nonlinear. A \textbf{nonlinear operator} is a mapping $\mathcal{F}$ that acts on a function or vector $\vec{x}$ but which \textit{does not} possess the properties of additivity and homogeneity: \[ \mathcal{F}(\alpha \vec{x} + \beta \vec{y}) \neq \alpha \mathcal{F}(\vec{x}) + \beta \mathcal{F}(\vec{y}) \] Roughly speaking, let us now define that we have some function $f(\vec{x})$ that depends on one or several parameters. And we want either: \begin{itemize} \item To solve the equation $f(\vec{x}) = 0$ (a system of nonlinear equations), \item To find the minimum $\min_{\vec{x}} f(\vec{x})$ (an optimisation problem). \end{itemize} \subsection{An example from chromatography: the plate model} A classical example is the system of equations of balance and dissociation constants on each chromatographic plate. We want to model a chromatograph, or rather its column. Separation occurs in it, and the mobile phase, moving along the column, at each step in fact creates such an equilibrium. You may say --- oh, there are only some 4--5 equations here, and just as many unknowns! Yes, I agree, not many. But: \begin{enumerate} \item There are many plates (say a thousand) --- that is one. \item We split the time of passage of the substance along the column into many time steps --- that is two. \end{enumerate} That is, at each time step we have a thousand times the same system of equations with different parameters of the input concentrations. And there will be substantially more such time steps than the number of plates, since the substance is after all retained by the column. That is, we must solve something like \textbf{ten million times} a system of equations of only 5 equations with 5 unknowns. But if we fail to solve it even once --- we cannot obtain any further solution. Unfortunately, such a system is not linear. The balance equations are linear, but the equations relating the dissociation constants contain products and ratios of the unknowns with respect to each other. And, surprisingly, such a system can sometimes indeed be solved very poorly. \subsection{Classification of nonlinear problems} Before speaking of methods, let us classify the problems themselves: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Property} & \textbf{Description} \\ \midrule \textbf{Smoothness:} & \\ \quad Smooth & The function has continuous derivatives (at least the first, and preferably the second). Example: $f(x) = x^2 + \sin x$. \\ \quad Almost smooth & The function is smooth almost everywhere, but there are break points or discontinuities of derivatives. Example: $f(x) = |x|$. \\ \quad Discontinuous & The function has discontinuities or is very noisy. Example: experimental data with artefacts. \\ \midrule \textbf{Number of minima:} & \\ \quad Unimodal & There is only one global minimum (or maximum). Example: convex functions. \\ \quad Multimodal & There are several local minima, and the task is to find the global one. Example: potentials of complex molecules. \\ \bottomrule \end{tabularx} \subsection{Methods for solving nonlinear problems} Now let us systematise the solution methods. Each method requires something, depends on something, and each has its strengths and weaknesses. \subsubsection{Monte Carlo methods} \textbf{Idea:} Random search. We generate random points in the parameter space, evaluate the function at them, and choose the best. \textbf{Requirements:} None. It can work even with discontinuous functions. \textbf{Pros:} \begin{itemize} \item Does not get stuck in local minima (with a sufficient number of iterations), \item Simple to implement, \item Parallelises perfectly. \end{itemize} \textbf{Cons:} \begin{itemize} \item Very slow convergence: the error decreases as $\mathcal{O}(1/\sqrt{N})$, where $N$ is the number of iterations, \item Does not use information about the structure of the function. \end{itemize} \textbf{When to use:} When the function is very noisy, discontinuous, or when one simply needs to find ``some'' solution to initialise other methods. \subsubsection{Gradient methods (Gradient Descent)} \textbf{Idea:} Move in the direction opposite to the gradient: \[ \vec{x}_{k+1} = \vec{x}_k - \alpha_k \nabla f(\vec{x}_k) \] where $\alpha_k$ is the learning rate. \textbf{Requirements:} The function must be differentiable (smooth). \textbf{Pros:} \begin{itemize} \item Simple to implement, \item Guaranteed convergence to a local minimum (with a correct choice of $\alpha$), \item Works well in high dimensions. \end{itemize} \textbf{Cons:} \begin{itemize} \item Gets stuck in local minima, \item Can converge very slowly in ``ravines'' (when the eigenvalues of the Hessian differ greatly), \item Requires choosing the step $\alpha$ (too large --- diverges, too small --- slow). \end{itemize} \subsubsection{Conjugate direction method (Conjugate Gradient for optimisation)} \textbf{Idea:} Construct search directions so that they are ``conjugate'' with respect to the Hessian of the function. This allows finding the minimum of a quadratic function in $N$ steps (where $N$ is the dimension). \textbf{Important remark:} The conjugate gradient method for solving linear systems iteratively is \textit{only named} the same as the conjugate gradient method for nonlinear minimisation. In fact, it is the same algorithm, simply applied to different problems: \begin{itemize} \item For linear systems $\mathbf{A}\vec{x} = \vec{b}$ --- this is the minimisation of the quadratic function $f(\vec{x}) = \frac{1}{2}\vec{x}^T\mathbf{A}\vec{x} - \vec{b}^T\vec{x}$, \item For nonlinear optimisation --- this is a generalisation of the same idea to arbitrary functions. \end{itemize} \textbf{Requirements:} A smooth function, preferably with a continuous gradient. \textbf{Pros:} \begin{itemize} \item Converges faster than ordinary gradient descent, \item Does not require storing the Hessian (unlike Newton's method), \item Memory --- $\mathcal{O}(N)$. \end{itemize} \textbf{Cons:} \begin{itemize} \item Gets stuck in local minima, \item For non-quadratic functions a ``reset'' of directions is required. \end{itemize} \subsubsection{Newton's methods} \textbf{Idea:} Use not only the gradient but also the second derivatives (the Hessian). We expand the function in a Taylor series up to second order: \[ 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} \] and find the minimum of this quadratic approximation. We obtain the iteration: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{H}^{-1}(\vec{x}_k) \nabla f(\vec{x}_k) \] where $\mathbf{H}$ is the matrix of second derivatives (the Hessian). \textbf{Requirements:} The function must be twice differentiable. The Hessian must be computable. \textbf{Pros:} \begin{itemize} \item \textbf{Quadratic convergence:} if we start close to the solution, the number of correct digits doubles at each iteration, \item Not sensitive to ``ravines'' (unlike gradient descent). \end{itemize} \textbf{Cons:} \begin{itemize} \item Requires computing the Hessian --- $\mathcal{O}(N^2)$ elements, \item Requires solving a linear system with the Hessian --- $\mathcal{O}(N^3)$ operations, \item May diverge if we start far from the solution, \item The Hessian may be degenerate or not positive definite. \end{itemize} \subsubsection{The Newton--Raphson method and block Newton for quantum mechanics} The \textbf{Newton--Raphson method} is a variant of Newton's method for solving systems of nonlinear equations $\vec{F}(\vec{x}) = \vec{0}$: \[ \vec{x}_{k+1} = \vec{x}_k - \mathbf{J}^{-1}(\vec{x}_k) \vec{F}(\vec{x}_k) \] where $\mathbf{J}$ is the Jacobian matrix (the matrix of first derivatives $\partial F_i / \partial x_j$). \textbf{Block Newton for quantum mechanics:} In quantum chemistry problems (for example, in the Hartree--Fock method or DFT) a system of equations often arises that can be split into blocks: \begin{itemize} \item Equations for the orbitals (expansion coefficients), \item Equations for the density, \item Equations for the energy. \end{itemize} The block Newton method takes this structure into account: the Hessian is represented in block form, and each block is processed separately. This allows: \begin{itemize} \item Efficient use of the structure of the problem, \item Parallelisation of computations, \item Applying different methods to different blocks (for example, Newton for one block, gradient descent for another). \end{itemize} \subsubsection{BFGS and L-BFGS methods} \textbf{BFGS (Broyden--Fletcher--Goldfarb--Shanno)} is a quasi-Newton method. The idea: not to compute the Hessian explicitly, but to approximate it at each iteration using information about the change of the gradient. \textbf{Requirements:} A smooth function with a continuous gradient. \textbf{Pros:} \begin{itemize} \item Convergence almost like Newton's, but without computing the Hessian, \item Guaranteed positive definiteness of the Hessian approximation, \item Works well in practice --- one of the most popular methods. \end{itemize} \textbf{Cons:} \begin{itemize} \item Requires storing an $N \times N$ matrix (the Hessian approximation) --- $\mathcal{O}(N^2)$ memory. \end{itemize} \textbf{L-BFGS (Limited-memory BFGS)} is a modification of BFGS for large problems. The idea: not to store the full matrix, but to store only the last $m$ pairs of vectors $(\vec{s}_k, \vec{y}_k)$, where: \[ \vec{s}_k = \vec{x}_{k+1} - \vec{x}_k, \quad \vec{y}_k = \nabla f_{k+1} - \nabla f_k \] \textbf{Pros:} \begin{itemize} \item Memory --- $\mathcal{O}(mN)$, where $m \sim 5\text{--}20$ (does not depend on $N^2$!), \item Ideal for large problems ($N > 1000$), \item Very popular in machine learning and quantum chemistry. \end{itemize} \textbf{Cons:} \begin{itemize} \item Converges slightly more slowly than full BFGS, \item Requires tuning the parameter $m$. \end{itemize} \subsubsection{Simplex methods (Nelder--Mead)} \textbf{Idea:} We construct a simplex (in $N$-dimensional space these are $N+1$ points), evaluate the function at the vertices, and iteratively reflect, contract or expand the simplex, moving towards the minimum. \textbf{Requirements:} None! The function can be nonsmooth, noisy, even discontinuous. \textbf{Pros:} \begin{itemize} \item Does not require gradients, \item Simple to implement, \item Works well for small dimensions ($N < 10$). \end{itemize} \textbf{Cons:} \begin{itemize} \item Very slow for large $N$, \item May get stuck, \item No theoretical guarantees of convergence. \end{itemize} \textbf{When to use:} When the function is very noisy or when $N$ is small and one is too lazy to derive gradients. \subsubsection{Simulated Annealing} \textbf{Idea:} Inspired by the physical process of annealing of metals. We start with a high ``temperature'' $T$, which allows the algorithm to ``jump over'' local minima. We gradually lower $T$, and the algorithm ``freezes'' in the global minimum. At each step: \begin{enumerate} \item We propose a random change $\vec{x} \to \vec{x}'$, \item We compute $\Delta f = f(\vec{x}') - f(\vec{x})$, \item If $\Delta f < 0$ --- we accept the change, \item If $\Delta f > 0$ --- we accept with probability $P = \exp(-\Delta f / T)$. \end{enumerate} \textbf{Requirements:} None. \textbf{Pros:} \begin{itemize} \item Can find the global minimum (with a correct ``annealing schedule''), \item Does not get stuck in local minima in the early stages, \item Simple to implement. \end{itemize} \textbf{Cons:} \begin{itemize} \item Very slow, \item Requires tuning the schedule $T(t)$, \item No guarantees of convergence in reasonable time. \end{itemize} \textbf{When to use:} When the problem is multimodal (many local minima) and the global one must be found. \subsection{Computing gradients} All gradient methods require computing $\nabla f(\vec{x})$. How is this done? \subsubsection{Finite differences} \textbf{Idea:} We approximate the derivative: \[ \frac{\partial f}{\partial x_j} \approx \frac{f(\vec{x} + h \vec{e}_j) - f(\vec{x})}{h} \] where $\vec{e}_j$ is the $j$-th basis vector. \textbf{Problems:} \begin{enumerate} \item \textbf{Expensive:} To compute the gradient in $N$-dimensional space we need $N+1$ evaluations of the function (or $2N$ for the central difference). \item \textbf{Unstable:} If $h$ is too large --- a large approximation error. If $h$ is too small --- loss of accuracy due to subtracting close numbers. The optimal $h \sim \sqrt{\epsilon_{\text{mach}}}$, where $\epsilon_{\text{mach}}$ is the machine precision. \item \textbf{Noise:} If the function is computed with noise (for example, experimental data), then finite differences amplify the noise. \end{enumerate} \subsubsection{Automatic differentiation (AD): analytic computation of the gradient} And now --- magic! It turns out that one can compute the gradient of a function \textit{analytically}, using automatic differentiation, and do it in only \textbf{a few times} the cost of computing the function itself! \textbf{Idea:} Any function $f(\vec{x})$ is a composition of elementary operations ($+$, $-$, $\times$, $\div$, $\sin$, $\cos$, $\exp$, $\log$, etc.). We can: \begin{enumerate} \item Represent the computation of the function as a \textbf{computational graph}, \item Apply the \textbf{chain rule} in reverse order (reverse mode). \end{enumerate} \textbf{Result:} The gradient is computed at $\mathcal{O}(1) \times$ the cost of the function (in practice --- 3--5 times more expensive, but not $N$ times!). \textbf{How it works in a nutshell:} \begin{enumerate} \item Forward pass: we compute $f(\vec{x})$, saving all intermediate results. \item Reverse pass: we go through the graph in the reverse direction, applying the chain rule: \[ \frac{\partial f}{\partial x_j} = \sum_{\text{paths}} \frac{\partial f}{\partial u_1} \frac{\partial u_1}{\partial u_2} \cdots \frac{\partial u_k}{\partial x_j} \] \end{enumerate} \textbf{Where it is used:} \begin{itemize} \item \textbf{PyTorch, TensorFlow:} The basis of all modern neural networks is precisely reverse-mode automatic differentiation. \item \textbf{Quantum chemistry:} Programs like PySCF, Psi4 use AD to compute energy gradients with respect to nuclear coordinates. \item \textbf{Optimisation:} All modern libraries (SciPy, JAX) use AD. \end{itemize} \begin{successbox}[Main conclusion] Never compute gradients by finite differences if automatic differentiation can be used! It is faster, more accurate and more stable. \end{successbox} \subsection{A live example: force fields and conformers} Another live example from chemistry is molecular modelling. For modelling small molecules one often uses a representation in which the total energy is written as a sum of simple interactions: \begin{enumerate} \item \textbf{Van der Waals interactions} between non-bonded atoms (Lennard-Jones potential): \[ 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{Harmonic bonds} between bonded atoms: \[ E_{\text{bond}} = \sum_{\text{bonds}} \frac{1}{2} k_b (r - r_0)^2 \] \item \textbf{Angle interactions} between three consecutive atoms: \[ E_{\text{angle}} = \sum_{\text{angles}} \frac{1}{2} k_\theta (\theta - \theta_0)^2 \] \item \textbf{Torsional (dihedral) interactions} between four consecutive atoms: \[ E_{\text{torsion}} = \sum_{\text{dihedrals}} \frac{V_n}{2} [1 + \cos(n\phi - \gamma)] \] \end{enumerate} There will not be many of them, but about the square of the number of atoms (for van der Waals). Each such function is some simple formula: sometimes with sines, sometimes with exponentials, sometimes with powers. And in total a molecule has $3N$ degrees of freedom, since each atom has 3 spatial coordinates. \subsubsection{The conformer problem} It would seem the task is clear: find the minimum of this function $E(\vec{x})$, where $\vec{x} \in \R^{3N}$ are the coordinates of all atoms. But here the most interesting begins. The energy function of a molecule is a \textbf{multimodal landscape} with a huge number of local minima. Each local minimum corresponds to its own \textbf{conformation} (conformer) of the molecule --- its own way of twisting the molecule in space. For example, for a protein of 100 amino acids the number of possible conformers can reach $10^{100}$ --- this is the famous \textbf{Levinthal paradox}. And the global minimum (the native structure) is only one of $10^{100}$ possibilities. \subsubsection{Why simple minimisation does not work} If we take a random initial conformation and run ordinary gradient descent or BFGS, we will very quickly (within a few hundred iterations) converge to the \textit{nearest} local minimum. But this will almost certainly \textit{not} be the global minimum --- that is, not the conformation the molecule adopts in reality. \subsubsection{Solution: a combination of methods} Therefore in practice a combination of methods is used: \begin{enumerate} \item \textbf{Simulated Annealing:} We start with a high ``temperature'', which allows the algorithm to jump over energy barriers and explore different regions of conformational space. We gradually lower the temperature, ``freezing'' the system in a low-energy region. \item \textbf{Local minimisation (BFGS):} After simulated annealing has found a ``good'' region, we run a fast local method (BFGS or the conjugate gradient method) for precise convergence to a local minimum. \item \textbf{Repetition:} We run this combination many times with different initial conditions and choose the conformer with the lowest energy. \end{enumerate} \begin{tipbox}[Popular programs] This approach is implemented in popular molecular dynamics programs: \begin{itemize} \item \textbf{AMBER, GROMACS, NAMD:} Use force fields (AMBER, CHARMM, OPLS) and methods like simulated annealing + local minimisation. \item \textbf{Rosetta:} For protein structure prediction it uses a fragment approach + Monte Carlo + local minimisation. \end{itemize} \end{tipbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{The fast Fourier transform: when $O(N^2)$ turns into $O(N \log N)$} \subsection{The pattern search problem} Sometimes we need to find a certain fragment in a huge sequence. If it is a very short fragment, then a simple enumeration of the whole sequence and comparison with this short fragment is a rather good strategy. This is often done when searching for short fragments in protein structures. But sometimes the dimensions of what must be compared are almost equal. For example, we have some periodic crystal structure, and we know that there is some deviation in it, and we even understand that the deviation is only partially similar to what we want to compare with. Formally: the initial data are $v_i$, $i = 1, \dots, N$, and what must be compared is $w_j$, $j = 1, \dots, M$, where $M < N$. We can explicitly compute in a loop: \[ \forall k = 0, \dots, N-M: \quad s_k = \sum_{j=1}^M w_j v_{j+k} \] and find the maximum over $s$. But the computational complexity of such a search will be quadratic in $N$: if $M \simeq N/2$, then we must perform $N^2/2$ arithmetic operations. For $N = 10^6$ this is $5 \times 10^{11}$ operations --- too many! \subsection{Matrix formulation} Here people came up with a way to find such $s_k$ much faster. Let us represent the computation of these $s$ in matrix form. Suppose we have the matrix: \[ \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} \] of size $(N-M+1) \times N$. If we multiply it by our vector $\vec{v}$, we obtain precisely our coveted $s_k$: \[ \vec{s} = \mathbf{W} \vec{v} \] But how do we do this quickly? \subsection{Circulant matrix} If we supplement this matrix from below with another $M-1$ rows of the form: \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*} then after multiplication we obtain the same vector $\vec{s}$, but at the end $M-1$ more numbers will be appended to it, which we can discard. But this extended matrix is a \textbf{circulant matrix} $\mathbf{C}$! It has a remarkable property: each successive row is obtained by a cyclic shift of the previous one. \subsection{Diagonalisation of a circulant matrix} A circulant matrix has an interesting representation via the Fourier matrix: \[ \mathbf{C} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \] where: \begin{itemize} \item $\vec{c}$ is the first column of the original matrix $\mathbf{C}$, \item $\mathbf{F}$ is the discrete Fourier transform (DFT) matrix, \item $\mathbf{F}^H$ is the Hermitian conjugate (complex conjugate and transposed) matrix. \end{itemize} The Fourier matrix is defined as: \[ \mathbf{F} = \{f_{jk}\}_{j,k=0}^{N-1}, \quad f_{jk} = e^{-2\pi i j k / N} \] \subsection{Fast multiplication via FFT} Now the most interesting thing. If we can quickly multiply the Fourier matrix by a vector, then we can compute this $\vec{s}$ faster: \[ \vec{s} = \mathbf{C} \vec{v} = \frac{1}{N} \mathbf{F}^H \diag(\mathbf{F} \vec{c}) \mathbf{F} \vec{v} \] This computation consists of three steps: \begin{enumerate} \item $\vec{a} = \mathbf{F} \vec{v}$ --- forward DFT of the vector $\vec{v}$, \item $\vec{b} = \diag(\mathbf{F} \vec{c}) \vec{a}$ --- element-wise multiplication (this is $\mathcal{O}(N)$), \item $\vec{s} = \frac{1}{N} \mathbf{F}^H \vec{b}$ --- inverse DFT. \end{enumerate} \subsection{Decomposition of the Fourier matrix} How do we quickly multiply $\mathbf{F}$ by a vector? It turns out that $\mathbf{F}$ can be represented as a product of special matrices: \begin{enumerate} \item One permutation matrix (bit-reversal permutation), \item Several pairs of matrices: \begin{itemize} \item Block matrices of the form $\begin{pmatrix} \mathbf{I} & \mathbf{I} \\ \mathbf{I} & -\mathbf{I} \end{pmatrix}$, \item Diagonal matrices with elements $e^{-2\pi i k / N}$ (the so-called ``twiddle factors''). \end{itemize} \end{enumerate} Let us write out these matrices explicitly for $N = 8$: \textbf{Permutation matrix (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} \] (rows permuted according to bit reversal: $0, 4, 2, 6, 1, 5, 3, 7$) \textbf{Block matrix (first level):} \[ \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{Diagonal matrix (twiddle factors, first level):} \[ \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) \] And so on for $\log_2 N$ levels. Then each such multiplication will require only $N$ and $2N$ arithmetic operations, and the total number of such steps will be $\log_2 N$. \textbf{Total:} multiplying the Fourier matrix by a vector can be done in $\mathcal{O}(N \log_2 N)$ arithmetic operations! And hence all $s_k$ can be computed in the same $\mathcal{O}(N \log_2 N)$, and not in $\mathcal{O}(N^2)$. For $N = 10^6$: \begin{itemize} \item Naive algorithm: $10^{12}$ operations, \item FFT: $2 \times 10^7$ operations (50\,000 times faster!). \end{itemize} \subsection{The fast Fourier transform (FFT)} This is the famous \textbf{fast Fourier transform} (FFT), discovered by Cooley and Tukey in 1965 (although Gauss knew it already in 1805!). It is one of the most in-demand algorithms in modern chemistry. \subsubsection{Application in NMR} With its help, for example, the initial FID (Free Induction Decay) from NMR is transformed into the spectral domain. If one looks at each row of the Fourier matrix, one can see an oscillating function in it: \begin{itemize} \item The closer the row is to the centre of the matrix, the greater the oscillation, \item The first row has none at all --- it is simply a constant, \item The very bottom (or second) row has an oscillation period equal to the length of the row. \end{itemize} That is precisely why, having multiplied the Fourier matrix by the FID, we obtain maxima where the resonant frequencies are: it is simply that such a row coincided in oscillation frequency with the resonant frequency. \begin{tipbox}[Why FFT is so important for NMR] Modern FIDs are very long --- they are often digitised into hundreds of thousands or even millions of numbers. If the fast Fourier transform did not exist, multiplication by such a matrix would last on modern computers even longer than the NMR spectra acquisition itself --- that is, really minutes and tens of minutes even on modern workstations. Thanks to FFT the transformation takes \textbf{fractions of a second}. \end{tipbox} \subsubsection{Other applications of FFT in chemistry} \begin{itemize} \item \textbf{Cryo-electron microscopy:} Reconstruction of 3D structures from 2D projections. \item \textbf{X-ray structure analysis:} Transformation of the diffraction pattern into electron density. \item \textbf{Molecular dynamics:} Computation of electrostatic interactions via the PPPM method (Particle-Particle Particle-Mesh). \item \textbf{Signal processing:} Noise filtering, data compression. \item \textbf{Chemometrics:} Fast convolution and correlation of spectra. \end{itemize} \begin{successbox}[Summary] FFT is one of those algorithms that changed the world. Without it there would be no modern spectroscopy, signal processing, audio and video compression, and much else. It is a must-know for any chemist working with data. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{When FFT does not see the obvious: spectral leakage and the Prony method} \subsection{The paradox: the eye sees a sinusoid, but FFT does not} Although FFT is a must-have algorithm almost everywhere in experimental chemistry, it too has its own features that one must definitely pay attention to. Let us consider an example. We have recorded a spectrum, the signal oscillates, and this is directly visible to the eye in the graph. But we managed to record only slightly more than one period. When we applied FFT, we obtained something very strange: the peak seems to be there, but also seems not to be --- it turned out somehow very smeared. And if we had a lot of other noise, then we obtained in another part of the spectrum something very similar, and yes, accidentally the noise turned out ``brighter'' than the signal, and our program produced a completely wrong result. But we see this sinusoid with our eyes! Why does FFT not see it? \subsection{Spectral leakage} FFT is good then, and \textit{only then}, when we have applied it to a signal in which \textbf{an exact integer number of periods} of this signal fits. Why? Because FFT implicitly assumes that our finite signal \textit{is periodically continued} beyond the measurement window. If exactly $k$ full periods fit in the window, then the periodic continuation will be smooth, and FFT will give one clear peak. But if the signal contains, say, $1.4$ periods, then in the periodic continuation a \textbf{discontinuity} arises at the boundary of the window. FFT is forced to approximate this discontinuity by a multitude of harmonics --- and energy ``leaks'' from the main peak into all other frequencies. This phenomenon is called \textbf{spectral leakage}. Mathematically: if the true frequency of the signal falls \textit{between} two neighbouring frequency bins of the FFT, then instead of one delta function we obtain a smeared ``hump'' (a function of the form $\sin(x)/x$, the so-called \textit{sinc}). \begin{warningbox}[Rule of FFT] FFT sees only those frequencies that fit into the measurement window an integer number of times. Everything else is leakage. \end{warningbox} \subsection{Zero-padding: a simple but not ideal solution} What should we do then? The simplest option is to add many zeros at the end of the signal (zero-padding) and hope that on the whole we will be substantially luckier. After all, if, for example, we add zeros so that the original data are stretched 5 times, then we will definitely get 7 periods (instead of 1.4). And if we increase, say, 7 times, then there will be 9.8 periods --- also very close to an integer. People often add a different number of zeros, and then try to combine the resulting spectra, guessing which peaks turned out more accurate in the corresponding extended spectra. This really helps! \begin{tipbox}[Important nuance of zero-padding] Zero-padding \textit{does not increase the real frequency resolution} --- it only interpolates the spectrum, making it smoother. The resolution is determined by the \textit{length of the original signal}, and not by the length of the zero-padded one. But in practice zero-padding helps to ``hit'' an integer number of periods and reduce leakage. \end{tipbox} \subsection{Frequency drift: when the signal ``drifts away''} But there is another problem that zero-padding cannot cope with. Sometimes we measure some spectrum for a long and tedious time, and it ``slightly'' took and drifted. That is, at the beginning we had frequency $\omega$, and towards the end $\omega + \varepsilon$, and this small correction is enough for us to completely miss with the FFT and obtain instead of a clear single peak something very smeared. But after all, one can see with the eye --- there is the sinusoid, both here and there, both at the beginning and at the end! Only FFT again does not see it. The reason is that FFT assumes \textbf{stationarity} of the signal --- that is, that the frequencies do not change with time. If the frequency drifts, then no harmonic of the FFT can approximate the whole signal well. What should we do? \subsection{The Prony method: the idea of shifts} Several centuries ago Gaspard de Prony (1795) answered this question for us. More precisely, he invented how to solve his own problem (expansion of a signal into a sum of exponentials), but it precisely solves our problem too --- when the spectrum drifts slightly in time. \textbf{Idea:} Take the signal and store it in a vector $\vec{s} = (s_0, s_1, \dots, s_{N-1})^T$. Then shift this vector down by one sample, then again, and put it side by side again. Do this $L$ times. What do we have here? If we had a periodic signal with some unknown sinusoid $\sin(\omega t)$, then after a shift by $p$ samples we obtain $\sin(\omega t + \omega p)$. Recalling the addition formulas from Chapter 1: \[ \sin(\omega t + \omega p) = \sin(\omega t)\cos(\omega p) + \cos(\omega t)\sin(\omega p) \] Since $\cos(\omega p)$ and $\sin(\omega p)$ are \textit{constants} for the whole vector (they do not depend on $t$), the set of such shifted vectors will have only \textbf{as many nonzero singular values as twice the number of distinct oscillating harmonics}. Moreover, if some frequency has ``drifted'' a little with time, this \textit{does not affect at all} such an approximation, but completely destroys the FFT result, as we noted above. \subsection{Prony method, variant 1: linear prediction} Let our signal be modelled as a sum of $K$ complex exponentials: \[ s_n = \sum_{m=1}^{K} c_m z_m^n, \quad n = 0, 1, \dots, N-1 \] where $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ are complex ``frequencies'' (containing both the damping $\alpha_m$ and the oscillation $\omega_m$), and $c_m$ are complex amplitudes. \textbf{Key fact:} If the signal consists of $K$ exponentials and there is no noise, then it satisfies a \textbf{linear recurrence relation} of order $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 \] This means that the $(n+1)$-th sample can be \textit{exactly predicted} from $K$ previous ones! The coefficients $a_1, \dots, a_K$ are the coefficients of the characteristic polynomial whose roots are the $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{How to find the coefficients?} Let us write the recurrence relation for all available samples in matrix form: \[ \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}}} \] This is an overdetermined system (if $N > 2K$), and we solve it by the least squares method: \[ \vec{a} = -(\mathbf{S}^H \mathbf{S})^{-1} \mathbf{S}^H \vec{s}_{\text{future}} \] After finding $\vec{a}$ we look for the roots of the polynomial $P(z)$ --- these are our $z_m$, from which the frequencies $\omega_m$ and dampings $\alpha_m$ are extracted. \subsection{Prony method, variant 2: SVD of the shift matrix} The second variant is more robust to noise and more elegant. \textbf{Step 1: Construct the Hankel matrix from the signal shifts.} \[ \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} \] of size $(N-L+1) \times L$, where $L > K$ (we choose $L$ definitely larger than the expected number of harmonics). \textbf{Step 2: Do SVD.} \[ \mathbf{H} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \] If the signal consists of $K$ exponentials and there is no noise, then $\rank(\mathbf{H}) = K$, and the last $L - K$ singular values are zero: \[ \sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_K > 0, \quad \sigma_{K+1} = \cdots = \sigma_L = 0 \] In the presence of noise the small singular values are not zero, but they are \textit{substantially smaller} than the first $K$. We discard them (this is called \textbf{truncated SVD} or \textbf{regularisation}). \textbf{Step 3: Extract the polynomial from the null space.} The singular vector $\vec{v}_{\min}$ corresponding to the \textit{smallest} singular value (the last column of $\mathbf{V}$) lies in the (approximate) null space of the matrix $\mathbf{H}$. This means: \[ \mathbf{H} \vec{v}_{\min} \approx \vec{0} \] If we write the components $\vec{v}_{\min} = (v_0, v_1, \dots, v_{L-1})^T$, this relation is equivalent to: \[ \sum_{j=0}^{L-1} v_j s_{n+j} \approx 0 \quad \forall n \] That is, $\vec{v}_{\min}$ contains the coefficients of the polynomial: \[ Q(z) = v_0 + v_1 z + v_2 z^2 + \cdots + v_{L-1} z^{L-1} \] \textbf{Step 4: Find the roots of the polynomial.} The roots of $Q(z)$ contain $K$ ``signal'' roots $z_m = e^{(\alpha_m + \i \omega_m)\Delta t}$ (lying inside or on the unit circle) and $L - K$ ``noise'' roots (scattered randomly). From the signal roots we extract: \[ \omega_m = \frac{\arg(z_m)}{\Delta t}, \quad \alpha_m = \frac{\ln|z_m|}{\Delta t} \] \begin{tipbox}[Why is the SVD variant better?] The SVD variant of the Prony method is more robust to noise, because: \begin{enumerate} \item Truncation of small singular values automatically filters the noise. \item We do not solve the overdetermined system directly (which may be ill-conditioned), but use an orthogonal decomposition. \item The number of harmonics $K$ is determined automatically from the ``jump'' in the spectrum of singular values. \end{enumerate} \end{tipbox} \subsection{Limitations of the Prony method} But Prony is not a panacea. As soon as a component that is well approximated both by the function itself and by its shifts (for example, white noise or a very broadband signal) appears in the original signal, we immediately obtain a ``root'' for it, and then it turns out to be \textbf{spurious}. Other limitations: \begin{itemize} \item The method assumes that the signal is a sum of \textit{exponentials} (including sinusoids as a special case). If the signal has a different structure (for example, rectangular pulses), the method will work poorly. \item The number of harmonics $K$ must be known or guessed in advance (although SVD helps with this). \item Finding the roots of a high-degree polynomial ($L > 50$) may itself be numerically unstable. \item The method is sensitive to strong noise: at a low signal-to-noise ratio spurious roots can ``masquerade'' as genuine ones. \end{itemize} At the same time this method is very actively and quite long used in NMR spectroscopy and complements FFT well. In practice a \textbf{hybrid approach} is often used: \begin{enumerate} \item FFT for a quick overview of the spectrum and determination of the number of peaks, \item The Prony method (or its modern variants: MUSIC, ESPRIT, Matrix Pencil) for precise determination of frequencies, dampings and amplitudes of individual peaks. \end{enumerate} \begin{successbox}[Summary] FFT is a powerful but ``blind'' tool: it sees only what fits into its frequency grid. The Prony method ``sees'' frequencies between bins and is robust to drift, but requires more computations and caution. In real NMR spectroscopy they work in tandem. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Least squares, residuals and Tikhonov regularisation} \subsection{Statement of the problem} So, we often solve problems that we call least squares problems. Let us consider some of them to understand more precisely how this can be solved. \subsection{An example from medicine: computed tomography (CT)} Let us take an example from an adjacent field --- from medicine. Suppose we have an X-ray CT scanner. We have placed the patient (or whatever we are investigating) immobile on some platform, and around him on opposite sides we place a transmitter and receivers of X-rays. We perform measurements of how strongly the radiation has been absorbed by the body under investigation, and we place these receivers and transmitters at various angles. Usually the patient lies on a couch, and the receiver and transmitters are located on the surface of a ring. This ring is positioned so that its axis is approximately along the couch, and the ring itself both rotates about its axis and moves along the couch. What is happening here? Soft tissues hardly absorb X-rays, while bones absorb quite strongly. We can mentally divide the whole space where the patient is located into small (usually equal) cubes --- \textbf{voxels} (volume elements). If we know exactly the position of the transmitter and receiver, we can draw a line (the X-ray beam) that intersects our cubes. Let each such cube have one unknown --- the intensity of absorption of X-rays. We number all these cubes and write the intensity of their absorption as an unknown vector $\vec{x} = (x_1, \dots, x_N)^T$. Then one mutual arrangement of transmitter and receiver gives us a rather sparse set of coefficients --- the lengths of passage of the beam through the corresponding cube. Let all this be written in a row of the matrix $\mathbf{A}$, and the number of columns in it equals $N$ (the number of voxels), and the number of rows equals the number of measurements performed (let it be $M$). But the measured values themselves (we may have, for example, one X-ray source and many receivers; then each mutual arrangement of receiver and source is one row) will be written in the vector $\vec{b}$, also of length $M$. Then we can formulate the least squares problem as: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2 \] that is, we want each row of the matrix $\mathbf{A}$, multiplied by $\vec{x}$ (the total absorption intensity), to be equal to what we measure. \subsection{Normal equations} We remember that such a problem can be solved, but let us look at it carefully. The residual function: \[ 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 \] If we find the derivative of $E$ with respect to all $x_j$ and equate it to zero (at the minimum), we obtain the system of equations: \[ 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} \] These are the so-called \textbf{normal equations}. Here everything seems fine: we obtain a positive semi-definite matrix $\mathbf{A}^H \mathbf{A}$ and some right-hand side, with which we can solve this problem. \subsection{The degeneracy problem} But what happens if we have constructed such cubes, but we had no experiment in which the X-ray beam passed through, for example, the $i$-th cube? Obviously, then in the original matrix $\mathbf{A}$ the corresponding $i$-th column will contain only zeros, and the matrix $\mathbf{A}^H \mathbf{A}$ will have zeros both in the $i$-th column and in the $i$-th row, that is, it will guaranteed have one zero singular value. Suppose this cube was somewhere in space, and there is no part of the body under investigation in it, that is --- fortunately for us --- we do not need this cube. But such a formulation \textit{will spoil} the solution: we simply cannot solve this problem, since the matrix $\mathbf{A}^H \mathbf{A}$ will be degenerate. Moreover, we cannot guarantee that even if $\mathbf{A}^H \mathbf{A}$ has no such zero columns and rows, its conditioning will be good. That is, instead of a solution we may obtain completely wrong numbers. \subsection{Solution via SVD and the pseudoinverse matrix} How should we act in this case? Let us return to the singular value decomposition of the matrix $\mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H$. We can notice that: \begin{itemize} \item If $\mathbf{\Sigma}$ has zeros on the diagonal, they relate to those cubes through which the X-ray beam did not pass. \item If there is something very small, it relates to those experiments where the X-ray beam only a few times and, possibly, only touching, hit some cube or set of cubes. \end{itemize} That is, we should refuse to take such data into account, but they are very difficult to filter out already at the stage of forming the matrix $\mathbf{A}$. But since these data are uninformative, let us simply ``throw them out''! That is, let us perform this SVD and zero out the small diagonal values in $\mathbf{\Sigma}$, and reconstruct $\mathbf{A}$. This, of course, is good, but we still do not know how to solve such a problem, since $\mathbf{A}^H \mathbf{A}$ will also contain zero singular values. But if for a square non-degenerate system $\mathbf{A}\vec{x} = \vec{b}$ one can represent the solution as: \[ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^H \quad \Rightarrow \quad \vec{x} = \mathbf{V} \mathbf{\Sigma}^{-1} \mathbf{U}^H \vec{b} \] then our least squares problem too we could probably solve by analogy. True, where $\mathbf{\Sigma}$ has zeros, in $\mathbf{\Sigma}^{-1}$ we will not invert them, but write a zero there too. Then it will be incorrect to call it inverse, and we shall call it \textbf{pseudoinverse}, denoting it as $\mathbf{\Sigma}^+$ (or $\mathbf{A}^+$ for the whole matrix). \subsection{The size problem: when SVD does not fit in memory} All would be fine, except that in real problems the matrix $\mathbf{A}$ can be huge. Often its elements are computed analytically, and there is no need to store them. And the dense matrix $\mathbf{U}$, whose dimension equals that of $\mathbf{A}$, usually does not fit in memory. For example, for an ordinary CT scan cubes of size about 1 mm are used. In a second about 100--500 measurements occur on 100--500 receivers, and the whole measurement lasts about a minute. That is, $N$ and $M$ can easily be of the order of 10--100 million, and the full singular matrix will require 100 petabytes of memory, which is unreasonably much. \subsection{Tikhonov regularisation} One can notice that if to the degenerate matrix $\mathbf{A}^H \mathbf{A}$ we add the identity $\lambda \mathbf{I}$, such that $\lambda$ is in the range of those very small singular values that we zeroed out, then such an addition can be written as: \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*} And we see that we add this lambda to each singular value, ``turning'' a degenerate or ill-conditioned matrix into a well-conditioned one. Moreover, for small $\lambda$ (of the order of the small singular values) we almost do not change the large singular values, but the small ones we simply replace by $\lambda$. Then, solving the system with such a matrix, we simply do not ``spoil'' the large singular values and substantially reduce the response of the small singular values. So, the conditioning of the matrix has become substantially smaller, and we do not ``spoil'' the solution with the matrix. Knowing that the matrix is sparse, we can apply some iterative methods (even the conjugate gradient method) and converge very quickly. In fact, such a $\lambda$ is the so-called \textbf{Tikhonov regulariser}, since instead of the original problem we can solve the following: \[ \min_{\vec{x}} \|\mathbf{A}\vec{x} - \vec{b}\|_2^2 + \lambda \|\vec{x}\|_2^2 \] in fact requiring that both the residual and the norm of the solution be minimised simultaneously, and their mutual proportion is regulated precisely by the value of this $\lambda$. There is much beautiful theory here, once derived in the 1970s by Tikhonov and his followers, but the essence remains simple: such an additional ``regulariser'' helps to solve \textbf{ill-posed} problems. \subsection{Iterative reduction of $\lambda$} Moreover, often at first $\lambda$ is set rather large, then the conjugate gradient converges very quickly (the conditioning of the matrix will be very small), and then $\lambda$ is reduced, each time using the previous solution as the initial value. This allows: \begin{itemize} \item Quickly obtaining a ``rough'' solution with a large $\lambda$, \item Gradually refining it by reducing $\lambda$, \item Avoiding problems with conditioning in the early stages. \end{itemize} \begin{successbox}[Summary] Tikhonov regularisation is a powerful tool for solving ill-conditioned least squares problems. It adds a penalty for the norm of the solution, which stabilises the solution and allows using iterative methods even for degenerate systems. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Basis functions: from finite differences to splines} \subsection{The idea: approximate the unknown through the known} Until now we have worked with discrete vectors and matrices. But real problems --- differential equations, integrals, wave functions --- live in continuous space. How do we reduce an infinite-dimensional problem to a finite-dimensional one? Answer: \textbf{basis functions}. We represent the unknown function $f(x)$ as a linear combination of known basis functions $\phi_k(x)$: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] and look for the coefficients $c_k$. This turns the problem of finding a function into the problem of finding a vector $\vec{c} \in \R^N$. \subsection{The Ritz method (variational method)} One of the most general approaches is the \textbf{Ritz method}. Suppose we have a functional $E[f]$ that we want to minimise (for example, energy in quantum mechanics). We substitute the expansion: \[ f(x) \approx \sum_{k=1}^N c_k \phi_k(x) \] and obtain a function of the coefficients: \[ E(\vec{c}) = E\left[\sum_{k=1}^N c_k \phi_k\right] \] Now we minimise $E(\vec{c})$ over $\vec{c}$ --- this is already an ordinary optimisation problem in $\R^N$. \subsection{Finite differences (Finite Difference Method, FDM)} The simplest choice of basis functions is \textbf{local constants} or \textbf{linear functions} on a uniform grid. \subsubsection{An example from chemistry: diffusion} Suppose we have the diffusion equation for the concentration of a substance $C(x,t)$: \[ \frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} \] We discretise space with step $h$: $x_j = jh$, and time with step $\Delta t$: $t_n = n\Delta t$. We approximate the second derivative by the central difference: \[ \frac{\partial^2 C}{\partial x^2}\bigg|_{x_j} \approx \frac{C_{j+1} - 2C_j + C_{j-1}}{h^2} \] where $C_j^n \approx C(x_j, t_n)$. Then the Euler scheme in time: \[ \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} \] This gives an explicit recurrence formula: \[ 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}[Stability] The explicit scheme is stable only for $\frac{D\Delta t}{h^2} \le \frac{1}{2}$. If the time step is too large --- the solution ``falls apart''. \end{warningbox} \subsection{Finite elements (Finite Element Method, FEM)} A more complicated but more flexible approach is \textbf{finite elements}. We divide the domain into small elements (triangles, tetrahedra) and on each element approximate the function by \textbf{local polynomials} (usually linear or quadratic). The basis functions are \textbf{piecewise linear ``hats''} (hat functions), each of which equals 1 at one node and 0 at all others. Advantages of FEM: \begin{itemize} \item Can work with complex geometries, \item Can adapt the grid (refine where the solution changes rapidly), \item Well theoretically grounded. \end{itemize} \subsection{Basis functions in quantum chemistry: Gaussian orbitals} And now --- the most interesting for chemists. In quantum chemistry (the Hartree--Fock method, DFT) we solve the Schrödinger equation: \[ \hat{H} \Psi = E \Psi \] where $\hat{H}$ is the Hamilton operator, $\Psi$ is the wave function. We represent the molecular orbitals $\psi_i(\vec{r})$ as a linear combination of \textbf{atomic orbitals} (LCAO --- Linear Combination of Atomic Orbitals): \[ \psi_i(\vec{r}) = \sum_{\mu=1}^K c_{\mu i} \phi_\mu(\vec{r}) \] where $\phi_\mu(\vec{r})$ are basis functions centred on the atoms. \subsubsection{Gaussian functions (Gaussian Type Orbitals, GTO)} In most modern programs (Gaussian, ORCA, Q-Chem) \textbf{Gaussian orbitals} are used as basis functions: \[ \phi_\mu(\vec{r}) = (x - X_A)^l (y - Y_A)^m (z - Z_A)^n \exp(-\alpha |\vec{r} - \vec{R}_A|^2) \] where: \begin{itemize} \item $\vec{R}_A = (X_A, Y_A, Z_A)$ are the coordinates of atom $A$, \item $l, m, n$ are the angular quantum numbers ($s, p, d, f$ orbitals), \item $\alpha$ is the exponent (controls the ``size'' of the orbital). \end{itemize} Why Gaussian functions? \begin{itemize} \item \textbf{The product of Gaussians is again a Gaussian:} This is critically important for computing multi-centre integrals (electron-electron repulsion), \item \textbf{Analytic integrals:} All integrals are computed analytically, \item \textbf{Contractions:} Real atomic orbitals are approximated by a sum of several Gaussians (contractions), which improves accuracy. \end{itemize} \subsubsection{The Ritz method in quantum chemistry} We substitute the LCAO expansion into the Hartree--Fock or DFT energy functional and minimise over the coefficients $c_{\mu i}$. This leads to the eigenvalue problem: \[ \mathbf{F} \vec{c}_i = \varepsilon_i \mathbf{S} \vec{c}_i \] where: \begin{itemize} \item $\mathbf{F}$ is the Fock matrix (or the Kohn--Sham matrix in DFT), \item $\mathbf{S}$ is the overlap matrix ($S_{\mu\nu} = \langle \phi_\mu | \phi_\nu \rangle$), \item $\varepsilon_i$ are the orbital energies, \item $\vec{c}_i$ are the expansion coefficients of the $i$-th molecular orbital. \end{itemize} This is a generalised eigenvalue problem, and it is solved iteratively (by the self-consistent field method, SCF). \begin{tipbox}[Basis sets] Popular basis sets in quantum chemistry: \begin{itemize} \item \textbf{STO-3G:} Minimal basis (3 Gaussians per atomic orbital), \item \textbf{6-31G*:} Double-zeta quality with polarisation functions, \item \textbf{cc-pVTZ:} Correlation-consistent triple zeta (high accuracy). \end{itemize} The larger the basis --- the more accurate the result, but the more expensive the computations ($\mathcal{O}(N^4)$ for Hartree--Fock, where $N$ is the number of basis functions). \end{tipbox} \subsection{Splines: smooth piecewise polynomial functions} And now --- about splines. These are basis functions widely used in signal processing, computer graphics and numerical methods. \subsubsection{What is a spline?} A \textbf{spline} is a piecewise polynomial function that is \textbf{smooth} (has continuous derivatives up to some order) at the stitching nodes. The most popular is the \textbf{cubic spline} (degree 3, smoothness $C^2$ --- the function, the first and the second derivatives are continuous). \subsubsection{B-splines (Basis Splines)} \textbf{B-splines} are a special set of basis splines with \textbf{compact support}. Each B-spline is nonzero only on a small interval. Advantages: \begin{itemize} \item \textbf{Locality:} Changing one coefficient affects only a small region, \item \textbf{Numerical stability:} The matrices are well conditioned, \item \textbf{Recursive definition:} B-splines are constructed recursively (de Boor's formula). \end{itemize} \subsubsection{Two-dimensional splines} B-splines are easily generalised to the two-dimensional case via the \textbf{tensor product}: \[ B_{ij}(x,y) = B_i(x) \cdot B_j(y) \] This is used in: \begin{itemize} \item Image processing (smoothing, interpolation), \item Finite elements (higher orders), \item Computer graphics (surfaces). \end{itemize} \subsubsection{Bézier splines} \textbf{Bézier splines} are parametric curves defined by \textbf{control points}. The curve does not pass through the control points, but is ``attracted'' to them. The formula of a Bézier curve of degree $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] \] where $\vec{P}_i$ are the control points, and $\binom{n}{i} (1-t)^{n-i} t^i$ are the \textbf{Bernstein polynomials}. Applications: \begin{itemize} \item Computer graphics (vector graphics, TrueType fonts), \item CAD systems (AutoCAD, SolidWorks), \item Animation and trajectories. \end{itemize} \subsubsection{Connection with finite elements} Interestingly, B-splines with compact support \textbf{smoothly flow} into basis finite elements. If we take B-splines of degree 0 --- these are piecewise constant functions (as in the simplest finite elements). Degree 1 --- piecewise linear ``hats''. Degree 2 and higher --- smoother bases. This allows constructing \textbf{isoparametric finite elements} of higher orders, which give a more accurate approximation with fewer nodes. \begin{successbox}[Summary] Basis functions are the bridge between the continuous world of differential equations and the discrete world of computers. Finite differences are the simplest, finite elements are the most flexible, Gaussian orbitals are specialised for quantum chemistry, and splines are for smooth interpolation and graphics. Understanding their properties is the key to choosing the right method for your problem. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Integration: from analytics to Monte Carlo} \subsection{Why do we need integration?} The integral is one of the most fundamental operations in mathematics and physics. We constantly compute: \begin{itemize} \item Areas and volumes, \item Mathematical expectations and variances, \item Normalisation constants in quantum mechanics, \item Multidimensional integrals in statistical physics, \item Convolutions in signal processing. \end{itemize} And if at school we learned to take integrals analytically, then in real life --- especially in chemistry and physics --- analytics often ends very quickly. \subsection{Analytic integration: when it works} Let us recall the basics. If we have a function $f(x)$ given analytically, and we know its antiderivative $F(x)$, then: \[ \int_a^b f(x)\,dx = F(b) - F(a) \] For example: \[ \int_0^1 x^2\,dx = \frac{x^3}{3}\bigg|_0^1 = \frac{1}{3} \] Beautiful, exact, fast. But: \begin{itemize} \item Not all functions have an elementary antiderivative (for example, $e^{-x^2}$ --- the error integral), \item Not all functions are given analytically --- often we have only a set of points (experimental data), \item In multidimensional cases analytics is almost always powerless. \end{itemize} \subsection{Numerical integration in 1D: Simpson's method} What should we do if the integral cannot be taken analytically? Answer: approximate the function by something simple and integrate that simple thing. In the previous chapter we learned to construct cubic splines. And here --- we also construct polynomials, and integrate with such pieces. \subsubsection{Rectangle method} The simplest approach: divide the interval $[a,b]$ into $N$ equal parts of width $h = (b-a)/N$, and on each interval approximate the function by a constant (the value at the midpoint): \[ \int_a^b f(x)\,dx \approx h \sum_{k=0}^{N-1} f\left(a + \left(k + \frac{1}{2}\right)h\right) \] The error is $\mathcal{O}(h^2)$. \subsubsection{Trapezoidal method} A little better: we approximate the function by a linear polynomial on each interval: \[ \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] \] The error is $\mathcal{O}(h^2)$. \subsubsection{Simpson's method} Even better: we approximate the function by a \textbf{quadratic polynomial} on each pair of intervals: \[ \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] \] where $N$ is an even number. The error is $\mathcal{O}(h^4)$! This is already very good for smooth functions. \begin{tipbox}[Why is Simpson so good?] Simpson's method is, in essence, the integration of piecewise quadratic splines. If the function is smooth (has continuous derivatives up to 4th order), then the error decreases as $h^4$, which is much faster than for trapezoids. For very smooth functions there are even more accurate methods --- Gaussian quadratures, which use optimally chosen nodes and weights. \end{tipbox} \subsection{The curse of dimensionality} But what happens if we integrate a 2-, 3-, 4-, 10-dimensional function? Then we need $N^{10}$ integration points, where $N$ is at least 100... That is somehow a lot! Let us compute. For a 10-dimensional integral with $N = 100$ points along each axis: \[ 100^{10} = 10^{20} \text{ points} \] Even if each point is computed in 1 nanosecond, this will take: \[ 10^{20} \text{ ns} = 10^{11} \text{ s} \approx 3000 \text{ years} \] This is the \textbf{curse of dimensionality}. \subsubsection{Where does this occur in chemistry?} \begin{itemize} \item \textbf{Quantum chemistry:} Computation of many-electron integrals (electron-electron repulsion) --- these are 6-dimensional integrals (3 coordinates per electron). \item \textbf{Statistical mechanics:} The partition function is an integral over the phase space of all particles. For $N$ particles --- this is a $6N$-dimensional integral (3 coordinates + 3 momenta per particle). \item \textbf{Molecular dynamics:} Averaging over configurations --- multidimensional integration. \item \textbf{Machine learning:} Bayesian inference --- integration over the space of model parameters. \end{itemize} \subsection{The Monte Carlo method for integration} And here the Monte Carlo method enters the stage --- one of the most elegant and surprising algorithms in computational mathematics. \subsubsection{The idea in a nutshell} Imagine that we want to compute the integral: \[ I = \int_{[a,b]^d} f(\vec{x})\,d\vec{x} \] where $d$ is the dimension (can be very large). The Monte Carlo method says: ``Let us simply throw $N$ random points $\vec{x}_1, \vec{x}_2, \dots, \vec{x}_N$ uniformly into the integration domain, evaluate the function at them and average''. The formula: \[ I \approx \frac{V}{N} \sum_{k=1}^N f(\vec{x}_k) \] where $V = (b-a)^d$ is the volume of the integration domain, and $\vec{x}_k$ are random points uniformly distributed in $[a,b]^d$. \subsubsection{Why does this work?} This is simply the law of large numbers. The mathematical expectation of the random variable $f(\vec{X})$, where $\vec{X}$ is uniformly distributed in $[a,b]^d$, equals: \[ \E[f(\vec{X})] = \frac{1}{V} \int_{[a,b]^d} f(\vec{x})\,d\vec{x} = \frac{I}{V} \] By the law of large numbers, the mean value $\frac{1}{N}\sum f(\vec{x}_k)$ converges to $\E[f(\vec{X})]$ as $N \to \infty$. \subsubsection{Rate of convergence} And here --- magic! The error of the Monte Carlo method decreases as: \[ \text{Error} \sim \frac{\sigma}{\sqrt{N}} \] where $\sigma$ is the standard deviation of the function $f(\vec{x})$ in the integration domain. \textbf{Important:} this rate of convergence \textit{does not depend on the dimension} $d$! For comparison: \begin{itemize} \item Simpson's method in 1D: error $\sim h^4 \sim N^{-4}$, \item Simpson's method in the $d$-dimensional case: error $\sim N^{-4/d}$ (catastrophically slow for large $d$), \item Monte Carlo method in any dimension: error $\sim N^{-1/2}$. \end{itemize} Yes, Monte Carlo converges more slowly than Simpson in 1D. But in the 10-dimensional case: \begin{itemize} \item Simpson: $N^{-4/10} = N^{-0.4}$, \item Monte Carlo: $N^{-0.5}$. \end{itemize} Monte Carlo is \textit{faster}! \begin{warningbox}[The price of Monte Carlo] The Monte Carlo method converges as $N^{-1/2}$ --- this means that to increase accuracy 10 times, 100 times more points are needed. This is slow by absolute standards, but it is the \textit{only} method that works in high dimensions. \end{warningbox} \subsection{An example from quantum chemistry: computation of orbitals} Let us see how the Monte Carlo method is applied in real quantum chemistry. \subsubsection{Problem: normalisation of the wave function} In quantum mechanics the wave function $\Psi(\vec{r}_1, \vec{r}_2, \dots, \vec{r}_N)$ describes the state of $N$ electrons. It must be normalised: \[ \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 \] This is a $3N$-dimensional integral! For a water molecule ($N = 10$ electrons) --- this is a 30-dimensional integral. \subsubsection{Application of Monte Carlo} We generate $M$ random configurations of electrons: \[ \{\vec{r}_1^{(k)}, \vec{r}_2^{(k)}, \dots, \vec{r}_{10}^{(k)}\}, \quad k = 1, \dots, M \] where each coordinate $\vec{r}_i^{(k)} = (x_i^{(k)}, y_i^{(k)}, z_i^{(k)})$ is chosen from some distribution (for example, Gaussian). We compute $|\Psi|^2$ in each configuration and average: \[ \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)} But that is not all! In quantum chemistry there is the \textbf{Variational Monte Carlo (VMC)} method, which uses Monte Carlo for the \textit{minimisation} of energy. We choose a trial wave function $\Psi_T(\vec{r}_1, \dots, \vec{r}_N; \vec{\alpha})$ with parameters $\vec{\alpha}$ and compute the energy: \[ 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} \] Both integrals are $3N$-dimensional, and we compute them by Monte Carlo. Then we minimise $E(\vec{\alpha})$ over the parameters $\vec{\alpha}$ --- and obtain an approximate wave function. \subsubsection{Quantum Monte Carlo (QMC)} There are even more advanced methods --- \textbf{Diffusion Monte Carlo (DMC)}, which solve the Schrödinger equation directly, modelling the ``diffusion'' of electrons in imaginary time. These methods give \textit{almost exact} solutions for small molecules and are used as a benchmark for testing other methods. \begin{tipbox}[Where is Monte Carlo used in chemistry?] \begin{itemize} \item \textbf{Quantum Monte Carlo (QMC):} Exact solution of the Schrödinger equation for small systems, \item \textbf{Monte Carlo molecular dynamics:} Sampling of configurations in statistical mechanics, \item \textbf{Integration in DFT:} Numerical integration of the exchange-correlation functional, \item \textbf{Bayesian optimisation:} Search for the global minimum of the energy of a molecule. \end{itemize} \end{tipbox} \subsection{Improvements of the Monte Carlo method} The basic Monte Carlo method is simply a uniform random walk. But there are many improvements: \subsubsection{Importance sampling} Instead of a uniform distribution we use a distribution proportional to $|f(\vec{x})|$. Then the points will more often fall into regions where the function is large, and the error will decrease. \subsubsection{The Metropolis method (Metropolis-Hastings)} For complicated multidimensional integrals we use Markov chains: each next point depends on the previous one. This allows efficient exploration of the space, even if it is very large. \subsubsection{Quasi-Monte Carlo} Instead of random points we use \textbf{deterministic} low-discrepancy sequences (Sobol, Halton sequences). They cover the space ``more uniformly'' than random points and converge faster: error $\sim N^{-1} (\log N)^d$. \begin{successbox}[Summary] The Monte Carlo method is the only practical way to compute multidimensional integrals. Its rate of convergence does not depend on the dimension, which makes it indispensable in quantum chemistry, statistical physics and machine learning. Although it converges more slowly than deterministic methods in low dimensions, in high dimensions it simply has no competition. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Statistics: how to separate information from noise} Until now we have mostly considered mathematical problems in which the numbers are assumed to be known. But experimental chemistry is arranged differently. If we measure the concentration of a substance, the position of a spectral line or the intensity of a signal, the instrument never reports the ``true'' value with infinite precision. Every measurement differs slightly from the others. Therefore, instead of a single number, we have to deal with a \textbf{random variable}. Suppose we measure the same quantity several times. We obtain: $(x_1,x_2,\ldots,x_N)=\vec{x}$, where $\vec{x}=\mu+\vec{\varepsilon},$ and \begin{itemize} \item $\mu$ is the unknown true or mean value; \item $\vec{\varepsilon}$ is the random measurement error. \end{itemize} This simple equation is one of the fundamental ideas of statistics. Statistics begins where we understand: \textit{we observe data, but we are interested in hidden quantities.} \subsection{The mean} The simplest estimate of $\mu$ is the arithmetic mean: % $$ \mu = \frac{1}{N} \sum_{i=1}^N x_i. $$ % Why exactly this? If the errors have a mean value near zero, then the arithmetic mean minimises the sum of squared deviations: $\min_{\mu} \sum_{i=1}^N |x_i - \mu|^2$, and the solution is given precisely by the formula above. By the way, this already explains why repeating an experiment usually increases the accuracy of the result. \subsection{Variance and standard deviation} The mean tells us where the centre of the data is, but says nothing about how strongly the results are scattered. For this one introduces the variance --- a measure of how strongly our measurements are ``scattered'' around the true value (to which the mean tends): $$ \sigma^2 = \frac{1}{N}\sum_{i=1}^N |x_i - \mu|^2 $$ Yes, we have just minimised it, have we not? Correct, but the value of the minimum will be equal to zero if and only if all the values $x_i$ were the same; otherwise --- this value is precisely the variance. And in order to measure it in the same units, we simply recall that the second norm contains a root, and write $$ \sigma = \sqrt{\frac{1}{N}\sum_{i=1}^N |x_i - \mu|^2} $$ That is, the variance is, in essence, the \textbf{square of the length of the vector of deviations from the mean}. Statistics again turns into linear algebra. \subsection{The normal distribution} In many physical and chemical measurements random errors are approximately described by the normal distribution: $$ p(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x-\mu)^2}{2\sigma^2} \right). $$ It is not necessary to memorise this formula. It is far more important to understand the meaning of the two parameters: $\mu$ and $\sigma$. The first determines the position of the centre of the distribution, the second --- its width, that is, the larger $\sigma$, the stronger the scatter of the measurements. \begin{tipbox}[Intuition] \begin{itemize} \item The mean answers the question: \textit{``Where is the result?''} \item The standard deviation answers the question: \textit{``How strongly are the results scattered?''} \end{itemize} \end{tipbox} \subsection{Why does the mean become more accurate when the experiment is repeated?} Suppose we have made $N$ independent measurements. For independent errors the variance of the mean is $\displaystyle \frac{\sigma^2}{N}$, consequently, \textbf{the error of the mean decreases as} $\frac{\sigma}{\sqrt N}$. Therefore increasing the number of measurements increases accuracy, but not linearly (as in Monte Carlo!). \begin{tipbox} To reduce the random error by a factor of $10$, in the idealised case approximately $100$ times more independent measurements are needed. \end{tipbox} \subsection{Random error and systematic error} Here a very important distinction must be made. If the instrument gives $x=\mu+\varepsilon$, where $\varepsilon$ fluctuates randomly around zero, repeated measurements can reduce the influence of this error. But let us imagine that the instrument systematically overestimates the result: $x=\mu+b+\varepsilon$, where $b$ is a constant offset. Then averaging gives $\bar{x}\approx\mu+b$. And no matter how many times we repeat the experiment, the offset $b$ will not disappear. \begin{warningbox}[Important] Repeating measurements reduces random error, but does not eliminate systematic error. Statistics cannot correct a wrongly calibrated instrument. \end{warningbox} This is one of the reasons why calibration, control samples and independent measurement methods are so important in chemistry. \subsection{Two quantities can be related} Suppose two quantities are measured simultaneously: $x_i$ and $y_i$. For example: % \begin{itemize} \item the concentration of a substance and the intensity of a signal; \item temperature and reaction rate; \item pressure and volume; \item two spectral characteristics. \end{itemize} We may be interested in the question: \textit{do $x$ and $y$ change in a coordinated way?} For this one uses the covariance: $\displaystyle \operatorname{Cov}(x,y) = \frac{1}{N}(x-\mu_x)^T (y-\mu_y)$ If large values of $x$ usually correspond to large values of $y$, the covariance is positive. If large values of $x$ correspond to small $y$, it is negative. If there is no linear relation, the covariance may be close to zero. \subsection{Correlation} The covariance depends on the units of measurement. Therefore one often uses the normalised quantity: $\displaystyle \rho_{xy} = \frac{\operatorname{Cov}(x,y)} {\sigma_x\sigma_y}.$ For it $\displaystyle -1\leq\rho_{xy}\leq1$. A value near $1$ means a strong positive linear relation, a value near $-1$ --- a strong negative one, and a value near $0$ means the absence of a pronounced linear correlation. But here one very important thing must be remembered: \begin{warningbox}[Correlation does not mean causation] If two quantities are correlated, this does not yet mean that one causes the other. Correlation speaks of a statistical relation between the data. To establish a causal relation, additional physical, chemical or experimental arguments are needed. \end{warningbox} \subsection{The covariance matrix} If we have not two but $p$ measured quantities, $x_1,x_2,\ldots,x_p$, then the covariances of all pairs can be collected into a single matrix: $$ \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}. $$ This matrix is symmetric: $\Sigma=\Sigma^T$, and it does not merely store statistical information. It is a mathematical object of linear algebra. That is precisely why eigenvalues and eigenvectors, which we have already encountered before, again arise in statistics. \subsection{PCA: statistics meets linear algebra} Let us imagine that we have $N$ chemical samples, for each of which $p$ characteristics have been measured. We obtain a data matrix: $X\in\mathbb{R}^{N\times p}$. Some characteristics may be strongly related to each other. For example, if two measured quantities almost always change together, the information about them is partially redundant. It would be desirable to find new coordinates in which: \begin{itemize} \item the first coordinate contains the maximum possible variation of the data; \item the second --- the maximum possible remaining variation; \item the third --- the next, and so on. \end{itemize} This is the main idea of \textbf{Principal Component Analysis} (PCA). Mathematically PCA is closely related to the eigenvectors of the covariance matrix: $\displaystyle \Sigma v_i=\lambda_i v_i$. The eigenvectors $v_i$ define new directions, and the eigenvalues $\lambda_i$ show how much of the variation of the data falls on the corresponding direction. If the first few eigenvalues are significantly larger than the rest, $$ \lambda_1,\lambda_2,\ldots,\lambda_k \gg \lambda_{k+1},\ldots,\lambda_p, $$ then most of the information can be represented by just a few coordinates. For example, $$ 5000\text{ measured parameters} \quad\longrightarrow\quad 20\text{ principal components}. $$ This does not mean that we have ``thrown away'' chemical information. We have found a more compact mathematical representation of the data. And here again the SVD appears, which we have already encountered: $$ X=U\Sigma V^T. $$ PCA and SVD turn out to be two sides of the same linear-algebraic idea. \subsection{Regression: when we want to build a model} Suppose we have measured the concentration $c_i$ and the corresponding signal $y_i$. We assume a linear model: $$ y_i=ac_i+b+\varepsilon_i. $$ Since the measurements contain errors, the points usually do not lie exactly on a single straight line. Therefore we look for $a$ and $b$ that minimise the sum of squared errors: $$ \min_{a,b} \sum_{i=1}^N \left(y_i-ac_i-b\right)^2. $$ This is the \textbf{least squares method}. Thus regression is not something completely new. We already know this mathematical construction: $$ \boxed{ \text{data} \rightarrow \text{model} \rightarrow \text{residual} \rightarrow \text{minimisation} } $$ That is precisely why linear algebra and optimisation are so important for statistics. \subsection{Overfitting} Now a more complicated problem arises. If there is little data, we can choose a very complicated function that will pass practically through every experimental point. For example, instead of a straight line one can use a polynomial of high degree: $$ y=a_0+a_1x+a_2x^2+\cdots+a_kx^k. $$ For a sufficiently large $k$ one can describe the available measurements almost perfectly. But this does not necessarily mean that we have found the correct physical law. We may simply have fitted the noise. This is called \textbf{overfitting}. That is precisely why in statistics and machine learning we are interested not only in the question: \textit{``How well does the model describe the known data?''} but also: \textit{``How well does it work on new data?''} \subsection{Training, validation and testing} If we build a model from the available data, it is useful to split the data into parts: $$ \boxed{ \text{training} \quad+\quad \text{validation} \quad+\quad \text{test} } $$ On the training set the model is trained. The validation set is used to choose the parameters of the model. The test set must remain independent and is used for the final evaluation of the model's ability to work on new data. This idea will be especially important when we reach machine learning. \subsection{What statistics is really trying to do} One can say that statistics deals not so much with ``counting averages'' as with a more general task: \textbf{to extract reliable information from imperfect data.} We observe: $$ \text{data} = \text{structure} + \text{noise}. $$ Our task is to reconstruct from the data the structure that interests us and at the same time to estimate how much we can trust it. In the simplest case this looks like this: $$ x_1,x_2,\ldots,x_N \longrightarrow \bar{x}\pm\text{uncertainty}. $$ In a more complicated case: $$ X \longrightarrow \text{PCA} \longrightarrow \text{low-dimensional representation}. $$ And even more complicated: $$ X \longrightarrow \text{model} \longrightarrow \text{prediction} \longrightarrow \text{validation on new data}. $$ It is precisely from here that modern machine learning naturally grows. \begin{successbox}[What is important to remember] \begin{enumerate} \item A real measurement contains random and, possibly, systematic error. \item The mean describes the central value of the data. \item The standard deviation describes their scatter. \item The random error of the mean decreases as $1/\sqrt N$. \item Systematic error is not eliminated by averaging. \item Covariance describes the joint variation of quantities. \item PCA looks for the most informative directions in multidimensional data. \item Regression builds a mathematical model from measurements. \item A good model must work not only on known data but also on new data. \end{enumerate} \end{successbox} \subsection{And again linear algebra} And, perhaps, the most pleasant observation for us is that statistics has turned out to be not such a foreign subject after all. We started with measurements: $\displaystyle x_1,x_2,\ldots,x_N,$ moved to vectors: $\displaystyle \vec x,$ then to norms: $\displaystyle \|\vec x\|_2,$ to covariance matrices: $\displaystyle \Sigma,$ to eigenvalues: $\displaystyle \Sigma v=\lambda v,$ to SVD: $\displaystyle X=U\Sigma V^T,$ and finally to optimisation: $\displaystyle \min \|Ax-b\|_2^2. $ That is, statistics does not destroy our mathematical picture. It shows how the very same mathematics begins to work with \textbf{real, imperfect and noisy data}. And this is precisely the situation that a modern chemist encounters almost every day. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{From SVD to compressed sensing: when there is less data than needed} \subsection{A practical problem: separation of spectra} Let us consider a practical problem. A mass of colourful organic matter was brought to our laboratory, and the only thing at hand was a VIS/IR spectrometer for fluorescence spectra. We were told that the samples contain several reaction products with presumably differing spectra, and we were asked to determine both the concentrations and the spectra of the pure substances. We chose a not very strong UV source so that the fluorescence was very bright and gave different spectra for different substances. This seems to be a problem for SVD! After all, if we put each such spectrum in its own column of the matrix $\mathbf{A}$, assume that we hypothetically have a matrix of pure spectra of these substances $\mathbf{P}$, and for each sample there are coefficients of concentration of the substances themselves $\mathbf{Q}$, then: \[ \mathbf{A} = \mathbf{P}\mathbf{Q} \] and we can obtain $\mathbf{P}$ as the left singular vectors of the matrix $\mathbf{A}$, and $\mathbf{Q}$ as the right singular vectors multiplied on the left by the diagonal. We can always also write the singular value decomposition as: \[ \forall i, j: \quad a_{ij} = \sum_{r=1}^R d_r u_{ir} v_{rj} \] Moreover, if we had a single substance ($R=1$), this is exactly how it would be. \subsection{But something went wrong...} Instead of beautiful pure spectra with positive values, we obtained some strange graphs in $\mathbf{P}$ --- with negative values, oscillations, ``ghost'' peaks. In fact --- this is exactly how it should have been, the spectra \textit{cannot} be obtained this way. The spectra themselves are positively defined graphs, so the scalar product of any pair, although it can be very close to zero, is not necessarily so. In the singular value decomposition we require that all columns of $\mathbf{U}$ be \textbf{orthogonal} to each other. And real spectra of substances are not orthogonal! They are simply ``similar'' or ``not similar''. We need something else to pull these spectra apart so that our mathematics sees them separately. \subsection{The magic of the third dimension: Kruskal's theorem} If we set up the experiment so that we obtain 3D data instead of 2D as in the example above, then it can be proved mathematically that a \textbf{non-orthogonal} decomposition also exists. For example, if we can record fluorescence spectra at several different excitation wavelengths. Then we will already have three-dimensional data: \[ \forall i, j, k: \quad a_{ijk} = \sum_{r=1}^R \alpha_r \, b_{ir} \, c_{jr} \, d_{kr} \] and, by \textbf{Kruskal's theorem} (Kruskal, 1977), we have no requirements on the orthogonality of each of these matrices! \subsubsection{What Kruskal's theorem says} Kruskal's theorem is one of the cornerstones of multidimensional data analysis. It gives conditions for the \textbf{uniqueness} of a tensor decomposition (the so-called CANDECOMP/PARAFAC decomposition). For a third-order tensor $\mathcal{A} = \sum_{r=1}^R \vec{a}_r \circ \vec{b}_r \circ \vec{c}_r$ (where $\circ$ is the outer product of vectors), the decomposition is \textbf{unique} up to permutation and scaling of components if: \[ k_A + k_B + k_C \geq 2R + 2 \] where $k_A$ is the \textbf{k-rank} (Kruskal rank) of the matrix $\mathbf{A}$, that is, the maximum number $k$ such that any subset of $k$ columns of the matrix $\mathbf{A}$ is linearly independent. \begin{tipbox}[Why is this so important?] In ordinary SVD we have \textit{infinitely many} decompositions $\mathbf{A} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^H$ --- any orthogonal transformation inside gives a new valid decomposition. Therefore SVD cannot isolate the ``physical'' components --- only mathematically convenient (orthogonal) ones. Kruskal's theorem says: for tensors of the 3rd and higher order, when the condition on the k-ranks is satisfied, the decomposition is \textit{unique}! This means that we can recover precisely those physical components that were in the data --- without requiring orthogonality. \end{tipbox} \subsubsection{Solution methods: PARAFAC and ALS} In practice the tensor decomposition is sought by iterative methods. The most popular is \textbf{ALS} (Alternating Least Squares): \begin{enumerate} \item We fix $\mathbf{B}$ and $\mathbf{C}$ and solve the least squares problem for $\mathbf{A}$, \item We fix $\mathbf{A}$ and $\mathbf{C}$ and solve for $\mathbf{B}$, \item We fix $\mathbf{A}$ and $\mathbf{B}$ and solve for $\mathbf{C}$, \item We repeat until convergence. \end{enumerate} This is a classical nonlinear minimisation problem (we discussed it in the chapter on nonlinear operators), and it may get stuck in local minima. Therefore in practice ALS is run many times with different initial approximations. \begin{orangebox}[Applicability in chemistry] PARAFAC/CANDECOMP is widely used in: \begin{itemize} \item Fluorescence spectroscopy (EEM --- Excitation-Emission Matrix), \item Chromatography-mass spectrometry, \item NMR spectroscopy, \item Analysis of metabolic data (metabolomics). \end{itemize} \end{orangebox} \subsection{Multidimensional NMR spectra} And now --- the most interesting for structural chemists. In NMR spectroscopy of proteins one can construct very multidimensional spectra: 2D, 3D, 4D and even 5D. Each dimension is its own frequency (the chemical shift of a certain nucleus: $^1$H, $^{13}$C, $^{15}$N). Cross-peaks in such spectra show which nuclei are close to each other in space, which allows reconstructing the 3D structure of the protein. But here is the problem: if we want to record a 4D spectrum with resolution $1024 \times 256 \times 256 \times 64$ points, then the total number of points is $\sim 4 \times 10^9$. And each measurement takes time (to wait for relaxation), and in total this is \textbf{weeks} of spectrometer operation! There was a time when, in order to determine all cross-peaks of ubiquitin (a peptide of 76 amino acids), one had to record a spectrum for about \textbf{three weeks}. The protein could degrade during that time! \subsection{Sparse sampling: less means more} But one can record only randomly sparse one-dimensional spectra --- and even not 10\%, and not 1\%, but really very, very few, and apply the same decomposition, only for sparse data, and obtain equally reliable results! Applying such methods 20 years ago, the author with his colleagues accelerated the recording of such a multidimensional spectrum without loss of information from three weeks to \textbf{15 minutes}. The idea is simple: if we know that the spectrum is \textit{sparse} in some basis (that is, consists of a small number of peaks), then we do not need to measure all points. It is enough to measure a random subset, and then \textit{recover} the missing points through optimisation. \subsection{Compressed Sensing: a revolution in measurements} This brings us to one of the most beautiful ideas in modern mathematics and signal processing --- \textbf{Compressed Sensing} (compressive sampling). \subsubsection{Classical theory: Nyquist--Shannon} The traditional theory of discretisation (Nyquist--Shannon) says: to reconstruct a signal with maximum frequency $f_{\max}$, one must measure it at a frequency of at least $2f_{\max}$. For a 4D NMR spectrum with resolution $1024 \times 256 \times 256 \times 64$ this means that we are \textit{obliged} to measure all $4 \times 10^9$ points. Fewer --- impossible, otherwise there will be aliasing (frequency folding). \subsubsection{The revolution: compressed sensing} But in 2004--2006 Donoho, Candès, Tao and other mathematicians showed: if a signal is \textbf{sparse} in some basis, then it can be reconstructed from a \textbf{substantially smaller} number of measurements! Formally: let the signal $\vec{x} \in \R^N$ have a sparse representation $\vec{x} = \mathbf{\Psi}\vec{s}$, where $\vec{s}$ has only $K \ll N$ nonzero elements. Then we can measure $\vec{y} = \mathbf{\Phi}\vec{x}$, where $\mathbf{\Phi} \in \R^{M \times N}$ is the measurement matrix with $M \sim K \log(N/K) \ll N$, and reconstruct $\vec{x}$ by solving the problem: \[ \min_{\vec{s}} \|\vec{s}\|_1 \quad \text{subject to} \quad \vec{y} = \mathbf{\Phi}\mathbf{\Psi}\vec{s} \] \begin{tipbox} \textbf{Why the L1 norm and not L0?} The idea is that we want to find the \textit{sparsest} solution (minimise $\|\vec{s}\|_0$ --- the number of nonzero elements). But minimising the L0 norm is an NP-hard problem (enumeration of all combinations). The magic of compressed sensing is that under certain conditions on the matrix $\mathbf{\Phi}$ (the so-called Restricted Isometry Property, RIP), minimising the L1 norm gives \textit{exactly the same result} as minimising L0! And L1 minimisation is a convex problem, which is solved efficiently (it is linear programming). \end{tipbox} \subsubsection{Conditions of applicability} For successful reconstruction two conditions are needed: \begin{enumerate} \item \textbf{Sparsity:} The signal must be sparse in some basis (for example, an NMR spectrum is a set of delta functions, that is, very sparse in the frequency domain). \item \textbf{Incoherence:} The measurement matrix $\mathbf{\Phi}$ must be ``incoherent'' with the basis $\mathbf{\Psi}$. In practice this means that the measurement points must be chosen \textbf{randomly} (or pseudorandomly). \end{enumerate} \subsection{Examples of Compressed Sensing in chemistry and beyond} \subsubsection{Fast MRI (Magnetic Resonance Imaging)} One of the best-known applications is the acceleration of MRI in medicine. Classical MRI requires lengthy scanning (k-space is filled line by line). With compressed sensing one can measure only \textbf{random} lines of k-space and reconstruct the full image via L1 minimisation. This reduces the scanning time by 5--10 times --- critical for patients who cannot lie still for long. \subsubsection{Sparse NMR spectroscopy} In multidimensional NMR of proteins compressed sensing allows: \begin{itemize} \item Measuring only a random subset of points in the multidimensional k-space, \item Reconstructing the full spectrum via L1 minimisation, \item Reducing the experiment time from weeks to hours or minutes. \end{itemize} This is precisely what the author of the script did 20 years ago, and what has now become the standard in modern NMR spectrometers (sparse sampling, Poisson-gap sampling, etc.). \subsubsection{Single-Pixel Camera (Rice camera)} An amazing example: a camera that has \textbf{no pixel matrix}, but only a single detector. The image is reconstructed via compressed measurements using a DMD (Digital Micromirror Device). This works because real images are sparse in the wavelet basis. \subsubsection{Mass spectrometry} In mass spectrometry compressed sensing is used to accelerate FT-ICR (Fourier Transform Ion Cyclotron Resonance) and Orbitrap --- one can measure the FID for a shorter time and still obtain high resolution. \subsubsection{Data compression} JPEG uses the discrete cosine transform (a close relative of Fourier), and images are sparse in this basis --- that is why JPEG compresses so well. Compressed sensing is the next step: we do not just compress already measured data, but immediately measure less. \subsection{Connection with previous chapters} Let us see how compressed sensing connects everything we have covered: \begin{itemize} \item \textbf{Linear algebra:} We solve an overdetermined system $\mathbf{\Phi}\mathbf{\Psi}\vec{s} = \vec{y}$ via L1 norm minimisation. \item \textbf{Conditioning:} The matrix $\mathbf{\Phi}$ must be well conditioned and satisfy the Restricted Isometry Property (RIP). \item \textbf{SVD and tensors:} For multidimensional data we use tensor decompositions (PARAFAC), which give uniqueness without orthogonality. \item \textbf{Nonlinear optimisation:} L1 minimisation is a convex problem, but not smooth (there is no derivative at zero). We use methods like ISTA (Iterative Shrinkage-Thresholding Algorithm) or ADMM. \item \textbf{FFT:} The operator $\mathbf{\Psi}$ is often the Fourier transform, and we use FFT for fast multiplication. \end{itemize} \begin{successbox}[Main conclusion] Compressed sensing is not just another algorithm. It is a \textit{philosophy of measurement}: if we know that the signal is simple (sparse), then we can measure it substantially less than classical theory requires. It has revolutionised MRI, NMR spectroscopy, mass spectrometry and many other fields. And all this is based on beautiful mathematics: convex optimisation, probability theory (random matrices) and linear algebra. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Information compression: from archivers to JPEG} \subsection{Two worlds of compression} Information compression is one of the most practical areas of applied mathematics. We constantly compress and decompress data: archives, pictures, video, music. But not all compression is the same. There are two fundamentally different approaches: \begin{itemize} \item \textbf{Lossless compression:} We can reconstruct the original data \textit{bit for bit}. Examples: ZIP, PNG, FLAC. \item \textbf{Lossy compression:} We sacrifice part of the information for the sake of greater compression. Examples: JPEG, MP3, MP4. \end{itemize} \subsection{Shannon's information entropy} Before speaking of methods, let us understand \textit{how much} data can be compressed at all. Claude Shannon in 1948 introduced the concept of \textbf{information entropy}. If we have an alphabet of $n$ symbols, and symbol $i$ occurs with probability $p_i$, then the entropy (the average amount of information per symbol) is: \[ H = -\sum_{i=1}^n p_i \log_2 p_i \quad \text{[bits]} \] This is the \textbf{theoretical limit} of lossless compression. We cannot compress data more strongly than $H$ bits per symbol. \begin{tipbox}[Example] In English text the letter ``e'' occurs most often ($p \approx 0.13$), and ``z'' --- very rarely ($p \approx 0.001$). If we encoded all letters equally (8 bits per symbol), we would spend too many bits on rare letters. The entropy of English text is about 4.2 bits per symbol, so theoretically we can compress it about 2 times. \end{tipbox} \subsection{Lossless compression: Huffman coding} The idea is simple: we notice that some symbols and sequences of 2--3 symbols occur more often, we set up a table of such symbols, and the more often a symbol or sequence occurs, the fewer bits we use to encode it. \subsubsection{The Huffman algorithm} \begin{enumerate} \item We count the frequencies of all symbols in the text. \item We build a binary tree: we merge the two rarest symbols into a node with the total frequency, and repeat until we obtain one tree. \item The code of a symbol is the path from the root to the leaf (0 --- left, 1 --- right). \end{enumerate} Frequent symbols end up closer to the root --- their codes are shorter. Rare ones --- farther, codes longer. \begin{tipbox}[Example] Let us have the symbols A (frequency 0.5), B (0.25), C (0.125), D (0.125). The Huffman tree will give the codes: \begin{itemize} \item A: 0 (1 bit) \item B: 10 (2 bits) \item C: 110 (3 bits) \item D: 111 (3 bits) \end{itemize} Average length: $0.5 \times 1 + 0.25 \times 2 + 0.125 \times 3 + 0.125 \times 3 = 1.75$ bits per symbol --- close to the entropy! \end{tipbox} \subsubsection{Dictionary methods: LZ77, LZ78, LZW} Algorithms of the LZ family (Lempel--Ziv) go further: they look not only for frequent symbols but also for \textbf{repeating sequences}. Idea: if we have already seen the string ``abracadabra'', then when it appears again we do not write it anew, but refer to the previous occurrence: ``(back 11, length 11)''. This is the basis of the ZIP, GZIP, PNG formats. \subsection{A chemical example: searching protein sequences} And now --- a live example from biology. We have a database of protein sequences (millions of amino acids), and we need to find whether it contains homologues (similar proteins) to our new protein. Direct comparison ``our protein vs every protein in the database'' is too slow. We need a smart algorithm for compression and search. \subsubsection{BLAST (Basic Local Alignment Search Tool)} BLAST is one of the most cited algorithms in biology. The idea: \begin{enumerate} \item We split our protein into short ``words'' (k-mers, usually $k = 3$ for proteins). \item For each word we look for similar words in the database (with small mutations). \item We extend the matches in both directions until the similarity falls below a threshold. \end{enumerate} This is, in essence, a dictionary compression method + fast search by a hash table. BLAST allows finding homologues in seconds, although full alignment would take hours. \begin{orangebox}[Why is this important?] BLAST and its variants (PSI-BLAST, BLASTP, BLASTN) are the basis of modern bioinformatics. Without them there would be no genome decoding, no drug development, no evolutionary studies. \end{orangebox} \subsection{Lossy compression: JPEG and low-rank approximation} Now --- lossy compression. The idea: we sacrifice part of the information that the human eye (or ear) will not notice anyway. \subsubsection{JPEG via the discrete cosine transform (DCT)} JPEG works like this: \begin{enumerate} \item We split the picture into $8 \times 8$ pixel blocks. \item For each block we apply the \textbf{discrete cosine transform} (DCT) --- a close relative of Fourier. \item DCT decomposes the block into 64 frequency components (from low frequencies --- the general background, to high frequencies --- fine details). \item We quantise the coefficients: divide by a quantisation matrix and round. High frequencies (fine details) are quantised more coarsely --- we ``lose'' them. \item The remaining nonzero coefficients are compressed losslessly (Huffman). \end{enumerate} \subsubsection{Low-rank approximation via SVD} An alternative view of image compression is via SVD. Let us have the image matrix $\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 \] We can approximate $\mathbf{A}$ by the first $r$ singular values: \[ \mathbf{A}_r = \sum_{k=1}^r \sigma_k \vec{u}_k \vec{v}_k^H \] This is the \textbf{low-rank approximation} of rank $r$. By the Eckart--Young theorem, this is the \textit{best} approximation of rank $r$ in the sense of the L2 norm. \begin{tipbox}[Compression ratio] The original matrix: $M \times N$ numbers. Approximation of rank $r$: $r(M + N + 1)$ numbers. If $r \ll \min(M,N)$, the compression is enormous! For example, for a $1000 \times 1000$ picture and $r = 50$: compression by $1000 \times 1000 / (50 \times 2001) \approx 10$ times. \end{tipbox} \subsection{Wavelets: local frequencies} Both DCT and SVD are global transformations: they decompose the whole picture into frequencies. But pictures have \textit{local} features: edges of objects, textures, points. \textbf{Wavelets} are basis functions that are localised both in space and in frequency. They allow analysing a signal at different scales. \subsubsection{The wavelet transform} Instead of sinusoids (as in Fourier) we use ``little waves'' (wavelets) --- functions that decay quickly. The most popular is the \textbf{Haar wavelet}: \[ \psi(t) = \begin{cases} 1, & 0 \leq t < 0.5 \\ -1, & 0.5 \leq t < 1 \\ 0, & \text{otherwise} \end{cases} \] The wavelet transform decomposes a signal by scales (frequencies) and positions. This allows: \begin{itemize} \item Compressing pictures better than JPEG (the JPEG2000 format uses wavelets), \item Detecting edges (edges of objects), \item Removing noise while preserving important details. \end{itemize} \begin{successbox}[Summary] Information compression is a balance between mathematical theory (entropy, SVD, wavelets) and psychophysics (what the eye sees, what the ear hears). From archivers to JPEG --- everywhere one idea: find structure in the data and use it for a compact representation. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Image comparison: from convolution to neural networks} \subsection{The problem: find an object in a picture} How do we compare two pictures and understand that they contain the same object? We have already compared two sequences by constructing a circulant matrix (the chapter on FFT). Probably pictures can be compared in the same way? Let them be 2D. But not everything here is as good as it seems. If one picture is simply \textbf{shifted} relative to the other along one or both axes --- yes, everything will work. We can use 2D convolution (via FFT) and find the maximum --- this will be the position of the object. But if there is a \textbf{rotation}? Or a \textbf{stretch}? Or a \textbf{change of illumination}? Simple convolution will no longer help. \subsection{Edge detection} The first step to understanding the content of a picture is the detection of \textbf{edges}. Edges are places where the pixel intensity changes sharply. \subsubsection{The Sobel operator} The simplest way is to compute the intensity gradient. The Sobel operator uses two $3 \times 3$ masks to compute derivatives with respect to $x$ and $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 \] where $*$ is convolution, $I$ is the image. Gradient modulus: $G = \sqrt{G_x^2 + G_y^2}$. Where $G$ is large --- there is an edge. \subsubsection{Canny edge detector} A more advanced algorithm: \begin{enumerate} \item We smooth the picture with a Gaussian filter (remove noise). \item We compute the gradient (as with Sobel). \item We suppress non-maxima (keep only the strongest edges). \item We apply hysteresis: if an edge is above the upper threshold --- keep, below the lower --- remove, between --- keep if connected to a strong edge. \end{enumerate} The result --- thin, connected lines of edges. \subsection{Keypoints: key points of an image} Edges are good, but we need something more ``point-like'' to compare pictures. For this one looks for \textbf{keypoints} (interest points) --- places that are easy to find and that are stable under transformations. \subsubsection{Harris corner detector} Idea: we look for corners --- places where an edge changes direction. For each pixel we compute the matrix of second moments of the gradient: \[ \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} \] where $w(x,y)$ is a window (for example, Gaussian). If both eigenvalues of $\mathbf{M}$ are large --- this is a corner. If one is large, the other small --- this is an edge. If both are small --- this is a flat region. \subsubsection{SIFT (Scale-Invariant Feature Transform)} SIFT is one of the most popular algorithms (Lowe, 1999). It finds keypoints that are stable under: \begin{itemize} \item Scale (the object may be closer or farther), \item Rotation, \item Change of illumination, \item Small changes of viewpoint. \end{itemize} Algorithm: \begin{enumerate} \item We build a \textbf{scale pyramid}: we blur the picture with a Gaussian filter with different $\sigma$ and reduce the size. \item We look for \textbf{extrema} in the difference of Gaussians (DoG --- Difference of Gaussians) --- this is an approximation of the Laplacian. \item For each keypoint we compute a \textbf{descriptor}: a histogram of gradients in a $16 \times 16$ pixel neighbourhood. This is a vector of 128 numbers. \item We compare the descriptors between pictures (Euclidean distance). \end{enumerate} \subsubsection{SURF, ORB, AKAZE} There are many improvements of SIFT: \begin{itemize} \item \textbf{SURF} (Speeded-Up Robust Features) --- faster, uses integral images. \item \textbf{ORB} (Oriented FAST and Rotated BRIEF) --- very fast, patent-free. \item \textbf{AKAZE} --- uses nonlinear diffusion scales. \end{itemize} \subsection{Finding the transformation: RANSAC} We have found keypoints in two pictures and matched them (feature matching). But not all matches are correct --- there are outliers. How do we find the transformation (rotation, scale, shift) that maps one picture into the other, if we have outliers? \subsubsection{RANSAC (Random Sample Consensus)} RANSAC is an elegant algorithm for robust parameter estimation: \begin{enumerate} \item We randomly choose the minimum number of points (for an affine transformation --- 3 pairs of points). \item From these points we compute the parameters of the transformation. \item We count how many \textit{of all} points agree with this transformation (inliers) --- distance less than a threshold. \item We repeat $N$ times, choose the transformation with the maximum number of inliers. \item Finally we refine the parameters over all inliers (least squares method). \end{enumerate} \begin{tipbox}[Why does RANSAC work?] If the fraction of outliers is not too large (say, $< 50\%$), then the probability of randomly choosing only inliers is nonzero. Repeating enough times, we will almost certainly find the ``right'' transformation. \end{tipbox} \subsection{Particle swarm methods and optimisation} Sometimes we have no explicit keypoints, but we know that one picture is a transformed version of the other. How do we find the parameters of the transformation? This is an optimisation problem: we maximise a \textbf{similarity metric} (for example, mutual information) over the parameters of the transformation. \subsubsection{Particle Swarm Optimization (PSO)} The particle swarm method is a heuristic optimisation method inspired by the behaviour of a flock of birds or a school of fish: \begin{enumerate} \item We initialise a ``swarm'' of $N$ particles --- each particle is a set of transformation parameters (for example, rotation angle, scale, shifts along $x$ and $y$). \item Each particle has a ``velocity'' and a ``position''. \item At each step the particle ``remembers'' its best position (personal best) and knows the best position of the whole swarm (global best). \item The velocity of the particle is updated as a weighted sum: inertia (continue moving in the same direction), cognitive component (attraction to personal best) and social component (attraction to global best). \item We repeat until convergence. \end{enumerate} PSO is especially useful when: \begin{itemize} \item The similarity function is nonsmooth or has many local maxima, \item The parameter space is large (for example, a 3D transformation with 12 parameters), \item We cannot compute the gradient (for example, mutual information is not differentiable). \end{itemize} \subsection{Neural networks: from hand-crafted features to learned} All the methods we discussed above (Sobel, Harris, SIFT) are \textbf{hand-crafted features}. We ourselves invent what an ``edge'', a ``corner'', an ``interesting point'' is, and how to build a descriptor from them. But what if we let the computer \textit{itself} learn to find features? \subsubsection{Convolutional neural networks (CNN)} \textbf{Convolutional Neural Networks} are a special type of neural networks created for working with images. The idea: \begin{enumerate} \item \textbf{Convolutional layers:} We apply a set of trainable filters (like Sobel masks, but the parameters are learned from data). Each filter extracts its own feature --- from simple edges in the first layers to complex textures and parts of objects in the deep layers. \item \textbf{Pooling layers:} We reduce the size of the feature map (max pooling --- take the maximum in a $2 \times 2$ window). This gives invariance to small shifts. \item \textbf{Fully connected layers:} At the end --- an ordinary neural network that classifies or regresses. \end{enumerate} \subsubsection{Hierarchy of features} What is surprising is that CNNs themselves learn the same hierarchy that we constructed by hand: \begin{itemize} \item \textbf{First layers:} Simple edges, gradients (like Sobel), \item \textbf{Middle layers:} Textures, corners, simple shapes (like SIFT descriptors), \item \textbf{Deep layers:} Parts of objects (eyes, wheels, leaves), \item \textbf{Last layers:} Whole objects (face, car, tree). \end{itemize} This is \textbf{feature learning} instead of hand-crafting. \subsection{Pretrained models and transfer learning} Training a CNN from scratch requires millions of labelled images and days of computation on a GPU. But there is a trick --- \textbf{transfer learning}. Idea: we take a network pretrained on a huge database (for example, ImageNet --- 1.4 million images, 1000 classes) and use it as a \textbf{feature extractor}. \subsubsection{Popular architectures} \begin{itemize} \item \textbf{VGG (2014):} Simple, deep (16--19 layers), good features. \item \textbf{ResNet (2015):} Very deep (up to 152 layers) with ``skip connections'' --- solves the problem of vanishing gradients. \item \textbf{EfficientNet (2019):} Optimised in all parameters (accuracy, speed, size). \item \textbf{Vision Transformer (ViT, 2020):} Uses the transformer architecture (as in NLP) for images. \end{itemize} \subsubsection{How to use a pretrained model} \begin{enumerate} \item \textbf{Feature extraction:} We take a pretrained network, cut off the last classification layer, and use the penultimate layer as a feature extractor. We obtain a vector of, say, 2048 numbers for each image. We compare the vectors (cosine distance). \item \textbf{Fine-tuning:} We take a pretrained network and fine-tune it on our data (for example, on microscopic images of cells). The first layers are ``frozen'' (they already work well), the last ones are trained. \item \textbf{Siamese networks:} Two identical networks with shared weights, trained so that similar images have close feature vectors, and different ones --- distant. \end{enumerate} \subsection{Applications in chemistry and biology} \subsubsection{Cryo-electron microscopy (cryo-EM)} In cryo-EM we obtain thousands of 2D projections of protein molecules in random orientations. The task: \begin{enumerate} \item Find all molecules on the micrographs (particle picking) --- this is an object detection task, \item Determine the orientation of each projection --- this is an image comparison task, \item Reconstruct the 3D structure --- this is an inverse tomography problem. \end{enumerate} Modern programs (RELION, cryoSPARC) use CNNs for particle picking and classification. This has allowed achieving a ``resolution revolution'' --- determining protein structures with atomic resolution. \subsubsection{Microscopy and cell analysis} In biological research one needs to: \begin{itemize} \item Count cells in microscopic images, \item Classify cell types, \item Track cell movement over time, \item Find anomalies (cancer cells). \end{itemize} CNN (especially U-Net for segmentation) is the standard in the field. \subsubsection{Chemometrics and drug discovery} In drug discovery one needs to compare molecular structures. But molecules are not pictures! However, they can be represented as 2D images (molecular graphs drawn on a grid) and CNN can be used to predict activity. \subsection{Connection with previous chapters} Let us see how image comparison connects everything we have covered: \begin{itemize} \item \textbf{Linear algebra:} Convolution is multiplication of a matrix by a vector (or tensor). CNN is a sequence of linear operations with nonlinearities. \item \textbf{SVD and low-rank approximation:} CNN can be compressed via SVD (low-rank factorisation of convolutional kernels). \item \textbf{FFT:} Convolution via FFT is $\mathcal{O}(N \log N)$ instead of $\mathcal{O}(N^2)$. In CNN this is critical for speed. \item \textbf{Nonlinear optimisation:} Training a CNN is the minimisation of a loss function (cross-entropy, MSE) via backpropagation + SGD/Adam (variants of gradient descent). \item \textbf{Compressed sensing:} In MRI and cryo-EM we use compressed measurements + CNN for reconstruction. \end{itemize} \begin{successbox}[Main conclusion] Image comparison is from simple to complex: \begin{enumerate} \item Convolution (for shifts) --- via FFT, \item Keypoints (SIFT) + RANSAC (for rotations and scale), \item Optimisation (PSO) for complex transformations, \item Neural networks (CNN) --- for semantic similarity. \end{enumerate} Each method is a compromise between universality, speed and accuracy. And all of them are based on the mathematics we have covered: linear algebra, optimisation, Fourier analysis, probability theory. \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Machine learning: from the perceptron to AlphaFold} \subsection{What is machine learning?} Machine Learning (ML) is a branch of artificial intelligence where the computer \textit{itself learns} to solve problems based on data, rather than by rigidly prescribed rules. Instead of writing a program ``if temperature > 100, then water boils'', we give the computer thousands of examples ``temperature --- state of water'' and ask it to find the pattern. \subsubsection{Three types of learning} \begin{itemize} \item \textbf{Supervised learning:} There are labelled data (input $\to$ output). The task is to learn to predict the output for new inputs. Examples: image classification, prediction of molecular properties. \item \textbf{Unsupervised learning:} Data without labels. The task is to find structure, clusters, hidden patterns. Examples: clustering of spectra, PCA. \item \textbf{Reinforcement learning:} An agent interacts with the environment and receives rewards/penalties. The task is to learn a strategy that maximises the reward. Examples: playing chess, controlling a robot. \end{itemize} \subsection{From the perceptron to neural networks} \subsubsection{The perceptron (1958, Rosenblatt)} The simplest ``neural network'' is the perceptron. It computes: \[ y = \sigma\left(\sum_{i=1}^n w_i x_i + b\right) \] where $x_i$ are the inputs, $w_i$ are the weights, $b$ is the bias, $\sigma$ is the activation function (for example, a step or a sigmoid). The perceptron can solve only \textbf{linearly separable} problems (for example, logical OR). For XOR (exclusive OR) one perceptron is not enough --- a multilayer network is needed. \subsubsection{Multilayer perceptron (MLP)} A network of several layers: an input layer, one or more hidden layers, an output layer. Each neuron is connected to all neurons of the next layer. Formula for the $l$-th layer: \[ \vec{h}^{(l)} = \sigma\left(\mathbf{W}^{(l)} \vec{h}^{(l-1)} + \vec{b}^{(l)}\right) \] where $\mathbf{W}^{(l)}$ is the weight matrix, $\vec{b}^{(l)}$ is the bias vector, $\sigma$ is a nonlinear activation function (ReLU, tanh, sigmoid). \subsubsection{Training: backpropagation} How to train a neural network? Minimise the loss function $L$ (for example, MSE for regression or cross-entropy for classification) via gradient descent. \textbf{Backpropagation} is an efficient algorithm for computing gradients via the chain rule. We go from the output to the input, computing $\partial L / \partial w_{ij}$ for each weight, and update the weights: \[ w_{ij} \leftarrow w_{ij} - \eta \frac{\partial L}{\partial w_{ij}} \] where $\eta$ is the learning rate. \subsection{Convolutional neural networks (CNN)} For images MLP is inefficient: too many parameters (each pixel is connected to all neurons). CNN use \textbf{convolutional layers}: a small filter (for example, $3 \times 3$) slides over the image and computes a convolution. This gives: \begin{itemize} \item \textbf{Locality of connections:} Each neuron looks only at a small region, \item \textbf{Weight sharing:} The same filter is applied to the whole image, \item \textbf{Shift invariance:} The object will be found anywhere. \end{itemize} Popular architectures: LeNet (1998), AlexNet (2012), VGG (2014), ResNet (2015), EfficientNet (2019). \subsection{Recurrent neural networks (RNN)} For sequences (text, time series) we need networks with \textbf{memory}. RNN pass the hidden state from one step to another: \[ \vec{h}_t = \sigma\left(\mathbf{W}_h \vec{h}_{t-1} + \mathbf{W}_x \vec{x}_t + \vec{b}\right) \] The problem of ordinary RNN is \textbf{vanishing gradients}: the network learns poorly on long sequences. Solutions: \begin{itemize} \item \textbf{LSTM} (Long Short-Term Memory, 1997): Adds ``gates'' that control the flow of information. \item \textbf{GRU} (Gated Recurrent Unit, 2014): A simplified version of LSTM. \end{itemize} \subsection{The revolution: transformers and attention} In 2017 the paper ``Attention Is All You Need'' (Vaswani et al.) was published, which turned NLP and much else upside down. \subsubsection{The attention mechanism} Idea: instead of compressing the whole sequence into one vector (as in RNN), we allow each element to ``look'' at all other elements and weight them by importance. For an input sequence $\vec{x}_1, \dots, \vec{x}_n$ we compute: \[ \text{Attention}(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^T}{\sqrt{d_k}}\right)\mathbf{V} \] where $\mathbf{Q}$ (query), $\mathbf{K}$ (key), $\mathbf{V}$ (value) are linear projections of the inputs. \subsubsection{Transformer} The Transformer architecture completely abandons recurrence and convolutions. Instead, only the attention mechanism (self-attention) and fully connected layers are used. Key components: \begin{itemize} \item \textbf{Multi-Head Attention:} Several parallel attention mechanisms, each of which learns to focus on different aspects of the data. \item \textbf{Positional Encoding:} Since the Transformer has no recurrence, it does not know the order of elements. We add positional encodings (sines/cosines or trainable vectors). \item \textbf{Layer Normalization:} Normalisation of activations for training stability. \item \textbf{Residual Connections:} ``Skip connections'' to combat vanishing gradients (as in ResNet). \end{itemize} Transformer became the basis for: \begin{itemize} \item \textbf{BERT} (2018): Bidirectional Encoder Representations from Transformers --- pretraining on masked words. \item \textbf{GPT} (2018--2023): Generative Pre-trained Transformer --- autoregressive text generation. GPT-3 (175 billion parameters), GPT-4 (multimodal). \item \textbf{T5, BART}: Encoder-decoder architectures for translation, summarisation. \end{itemize} \subsection{Machine learning in chemistry: evolution} \subsubsection{The QSAR era (1990s -- 2010s)} QSAR (Quantitative Structure-Activity Relationship) is a classical approach: we compute \textbf{molecular descriptors} (numerical characteristics of a molecule) and build a regression/classification. Descriptors: \begin{itemize} \item Physicochemical: molar mass, logP (lipophilicity), polar surface area, \item Topological: connectivity indices, shapes of the molecular graph, \item Electronic: atomic charges, orbital energies, \item 3D descriptors: moments of inertia, radii of gyration. \end{itemize} Methods: PLS (Partial Least Squares), Random Forest, SVM (Support Vector Machines). \begin{tipbox}[Example] Predicting the toxicity of a molecule: we compute 200 descriptors, collect a database of 10\,000 molecules with known toxicity, train a Random Forest. Accuracy --- 80--85\%. \end{tipbox} \subsubsection{Graph neural networks (2015 -- present)} A molecule is a \textbf{graph}: atoms are nodes, bonds are edges. Graph Neural Networks (GNN) work directly with graphs. \textbf{Message Passing Neural Networks (MPNN):} \begin{enumerate} \item Each atom has an initial representation (a feature vector: atom type, charge, hybridisation). \item At each step atoms ``exchange messages'' with neighbours through bonds. \item After $K$ steps each atom ``knows'' about its neighbourhood of radius $K$ bonds. \item The final representation of the molecule is an aggregation of all atomic representations. \end{enumerate} Popular architectures: \begin{itemize} \item \textbf{GCN} (Graph Convolutional Network): An analogue of convolution for graphs. \item \textbf{GAT} (Graph Attention Network): Attention between atoms. \item \textbf{SchNet, DimeNet, SphereNet}: Take into account 3D coordinates of atoms and angles. \end{itemize} \begin{successbox}[Why are GNN so good for chemistry?] \begin{itemize} \item \textbf{Permutation invariance:} The order of atoms does not matter, \item \textbf{Structure awareness:} GNN see bonds and geometry, \item \textbf{Interpretability:} One can look at which atoms/bonds are important for the prediction. \end{itemize} \end{successbox} \subsubsection{Prediction of molecular properties} Modern must-have models for a chemist: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Task} & \textbf{Model} & \textbf{Accuracy} \\ \midrule Molecular energy & SchNet, DimeNet++ & MAE $\sim$1 kcal/mol \\ Solvation & 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{The AlphaFold revolution: protein structure prediction} \subsubsection{The problem} Predicting the 3D structure of a protein from its amino acid sequence is one of the grand challenges of biology. Experimental methods (X-ray, cryo-EM, NMR) are expensive and slow. \subsubsection{AlphaFold (2020, DeepMind)} AlphaFold 2 is a breakthrough that solved the problem of protein structure prediction with atomic accuracy. Architecture: \begin{enumerate} \item \textbf{Evoformer:} A Transformer that processes: \begin{itemize} \item Multiple sequence alignment (MSA) --- evolutionary information, \item Pair representation --- pairwise distances between residues. \end{itemize} \item \textbf{Structure Module:} Iteratively builds 3D coordinates of atoms, minimising the loss function. \item \textbf{Confidence Estimation:} Predicts pLDDT (per-residue confidence) --- how confident the model is about each residue. \end{enumerate} Results on CASP14 (Critical Assessment of Structure Prediction): \begin{itemize} \item Average GDT\_TS (Global Distance Test) --- 92.4 (experimental accuracy --- $\sim$90), \item For 2/3 of proteins --- accuracy $< 1$ Å (atomic accuracy). \end{itemize} \subsubsection{AlphaFold 3 (2024)} Extension to: \begin{itemize} \item Protein--ligand complexes, \item DNA/RNA structures, \item Post-translational modifications, \item Ionic interactions. \end{itemize} \subsection{A foundation model for conformers} \subsubsection{The conformational space problem} A molecule is not a static structure. It constantly ``breathes'', rotates around bonds, changes conformation. For drug design it is critically important to know not one structure, but an \textbf{ensemble} of low-energy conformers. Classical methods: \begin{itemize} \item \textbf{Molecular Dynamics (MD):} We model the motion of atoms in time, but it is slow (nanoseconds --- microseconds). \item \textbf{Monte Carlo:} Random walk over conformational space, but inefficient for large molecules. \item \textbf{Systematic/Rotamer search:} Enumeration of all combinations of torsion angles, but exponential growth. \end{itemize} \subsubsection{ML approach: conformer prediction} Modern models learn to predict the distribution of conformers directly from data. \textbf{GeoLDM (Geometric Latent Diffusion Models, 2023):} \begin{enumerate} \item \textbf{Encoder:} Encodes a 3D conformation into a latent space. \item \textbf{Diffusion Process:} Adds noise to the latent vectors (as in DALL-E, Stable Diffusion). \item \textbf{Decoder:} Generates new conformations from noise through the reverse diffusion process. \item \textbf{Energy Model:} Filters the generated conformations by energy. \end{enumerate} \textbf{ConfGF (Conformation Graph Factor, 2022):} \begin{itemize} \item Uses GNN to predict 3D coordinates, \item Generates an ensemble of conformers via sampling, \item Takes into account the Boltzmann distribution (low-energy conformers are more probable). \end{itemize} \subsubsection{Foundation model: Uni-Mol (2023)} Uni-Mol is a ``foundation model'' for molecules, pretrained on millions of conformations. Architecture: \begin{enumerate} \item \textbf{3D Transformer:} Works with atomic coordinates and features. \item \textbf{Pre-training tasks:} \begin{itemize} \item Prediction of distances between atoms, \item Prediction of angles and dihedral angles, \item Masked atom prediction (as in BERT), \item Contrastive learning (similar conformers --- close). \end{itemize} \item \textbf{Fine-tuning:} Fine-tuning on specific tasks (energy, properties, conformers). \end{enumerate} Results: \begin{itemize} \item Energy prediction: MAE $\sim$0.5 kcal/mol (better than DFT for some classes), \item Conformer generation: RMSD $< 0.5$ Å for 90\% of molecules, \item Property prediction: state-of-the-art on most benchmarks. \end{itemize} \subsection{Must-have tools for the modern chemist} \subsubsection{Software} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Tool} & \textbf{Purpose} \\ \midrule RDKit & Chemoinformatics: descriptors, fingerprints, similarity \\ DeepChem & ML for drug discovery, GNN \\ PyTorch Geometric & Graph neural networks \\ OpenMM & Molecular dynamics on GPU \\ ASE (Atomic Simulation Environment) & Quantum chemistry + ML \\ SchNetPack & SchNet and other GNN for chemistry \\ AlphaFold (ColabFold) & Protein structure prediction \\ Uni-Mol & Foundation model for molecules \\ \bottomrule \end{tabularx} \subsubsection{Typical workflow} \begin{enumerate} \item \textbf{Data collection:} PubChem, ChEMBL, PDB (Protein Data Bank). \item \textbf{Preprocessing:} RDKit for descriptors, generation of conformers. \item \textbf{Modelling:} GNN for property prediction, AlphaFold for structure. \item \textbf{Validation:} Cross-validation, external test set, experimental verification. \item \textbf{Interpretation:} Attention weights, SHAP values for understanding predictions. \end{enumerate} \subsection{Connection with previous chapters} Let us see how ML connects everything we have covered: \begin{itemize} \item \textbf{Linear algebra:} Neural networks are a sequence of matrix multiplications $\mathbf{W}\vec{x} + \vec{b}$. Backpropagation is the chain rule for matrices. \item \textbf{SVD and tensors:} Model compression via low-rank factorisation. Tensor decompositions for efficient computations. \item \textbf{FFT:} Convolutional layers use FFT for acceleration. Attention via FFT (Linear Attention). \item \textbf{Nonlinear optimisation:} Training neural networks is the minimisation of a loss via SGD, Adam (adaptive methods), L-BFGS (for fine-tuning). \item \textbf{Compressed sensing:} Sparse training, pruning of neural networks. \item \textbf{Graphs and tensors:} GNN work with molecular graphs. 3D models use tensors of coordinates. \item \textbf{Monte Carlo:} Diffusion models are stochastic processes. Variational inference --- a Bayesian approach. \end{itemize} \begin{successbox}[Main conclusion] Machine learning is not a replacement for classical chemistry, but a \textbf{powerful tool} that: \begin{itemize} \item Accelerates computations thousands of times (property prediction vs DFT), \item Opens new possibilities (protein structure prediction, molecule generation), \item Requires understanding of mathematics (linear algebra, optimisation, probability theory). \end{itemize} The modern chemist must know: \begin{itemize} \item The basics of ML (neural networks, GNN, transformers), \item The tools (RDKit, PyTorch, DeepChem), \item When ML works and when it does not (physical limitations, interpretability). \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Modern computing hardware: from the transistor to the supercomputer} \subsection{Numbers worth knowing} A modern computer --- a workstation --- is: \begin{itemize} \item 128--1024 GB of RAM, \item $\sim$500 Gflop/s of peak performance (billions of floating-point operations per second), \item Several processors (from 4 to 128 cores). \end{itemize} But here is a paradox: the speed of memory access is \textbf{hundreds of times} slower than the speed of arithmetic operations. And if access is random (not sequential read-write), then another \textbf{hundred times} slower. Why is that? Let us sort it out. \subsection{How a processor is arranged} Everyone says that a processor has ``many transistors''. But what do they do? In fact, a processor is an entity that accesses \textbf{instructions} and accesses \textbf{data}. Both take up space in memory. \subsubsection{Processor commands} The processor reads commands, which usually allow: \begin{itemize} \item Accessing data in memory and placing them in \textbf{registers} --- a temporary, very, very fast memory, \item Performing an operation: addition, subtraction, multiplication, comparison, \item Making a \textbf{jump}: by default the processor executes command after command, but if it encounters a jump instruction, it can leap forward or backward. \end{itemize} Jumps are: \begin{itemize} \item \textbf{Unconditional:} The jump always happens, \item \textbf{Conditional:} If we have previously compared numbers and have the answer ``yes'' or ``no'', then the jump happens only if the condition is satisfied. \end{itemize} \subsubsection{SIMD: one command --- many data} A modern processor command can work with several numbers at once. For example, one AVX-512 instruction can perform 16 pairs of additions simultaneously! This is called \textbf{SIMD} (Single Instruction, Multiple Data) --- one instruction, many data. If the data are already on the processor itself, this is quite easy to do. \begin{tipbox}[Where did SIMD come from?] People noticed that one often has to perform identical operations: for example, transform pixels identically to draw a picture, process array elements identically. Why execute the same instruction 16 times if one can do it once --- but for 16 numbers at once? \end{tipbox} But to perform such operations with \textit{different} data with each command, we would need a very powerful system for pumping data from memory to the processor and back. \subsubsection{The memory wall problem} Unfortunately, such a system would occupy hundreds and thousands of times more transistors than are needed for an arithmetic block. And people noticed that in programs we often work with some small set of numbers, and then switch to another set, and so on. Making a small block of memory with fast access is not so costly in terms of the number of transistors. The only trouble is that forcing the programmer to drag data back and forth into this block from RAM and back would be difficult. Therefore people automated this --- and the \textbf{processor cache} appeared. It is fast memory, there is not much of it, but the processor itself ``remembers'' what you are using now, and for some time keeps these data in its cache, expecting that later they may be needed again. \subsection{Memory hierarchy} A modern computer is a multilayer memory structure: \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Level} & \textbf{Size} & \textbf{Speed} & \textbf{Description} \\ \midrule Registers & 32--128 $\times$ 64 bits & $\sim$0.3 ns & Cells inside the processor that it works with directly \\ L1 cache & $\sim$32--64 KB & $\sim$1 ns & Separate for instructions and data, own for each core \\ L2 cache & $\sim$256 KB -- 1 MB & $\sim$3--5 ns & Own for each core or shared by a pair of cores \\ L3 cache & $\sim$8--64 MB & $\sim$10--20 ns & Shared by all cores of the processor \\ RAM (DRAM) & 128--1024 GB & $\sim$50--100 ns & Main memory \\ Disk (SSD) & 1--10 TB & $\sim$10--100 $\mu$s & Long-term storage \\ \bottomrule \end{tabularx} Let us compare how much a modern workstation can on average do in one nanosecond: \begin{itemize} \item \textbf{Arithmetic:} All its processors, taking into account long SIMD instructions, can perform about \textbf{1000 floating-point operations}. \item \textbf{L1 cache:} One can read and write 2--3 numbers on each processor, in total $\sim$200 numbers. \item \textbf{L2/L3 cache:} This quantity falls by tens of times. \item \textbf{RAM:} Only a few reads and writes per nanosecond, and only if this happens in large blocks of about 512 bytes approximately sequentially. \item \textbf{Random access:} If each time we address different blocks, the speed falls by another tens of times --- several nanoseconds may be needed per number. \end{itemize} \begin{warningbox}[Memory Wall] Random access to memory is \textbf{tens of thousands of times} slower than the real computational power of a modern processor! This is the main bottleneck of modern computations. That is precisely why: \begin{itemize} \item Matrix algorithms work faster if the data are located ``densely'' in memory, \item Sparse matrices are slow (many random accesses), \item Cache-oriented algorithms are a separate science. \end{itemize} \end{warningbox} \subsection{Graphics processing units (GPU)} But that was not enough for people. When they drew 3D objects of computer graphics on the screen, they noticed that all this can be done almost entirely \textbf{in parallel}. Roughly: if we had 1024 processors, we would split the whole screen into a grid of $32 \times 32$ squares, and each processor would draw in its own square. Thus graphics cards appeared --- they have a host of small specialised processors that work fully in parallel. \subsubsection{From graphics to computations} In the late 1990s such graphics cards were used to accelerate the display of graphics. Around 2005 people began to produce such cards so that small pieces of code could be executed on thousands of such small processors --- thus appeared \textbf{CUDA} (NVIDIA) and \textbf{OpenCL} (an open standard). By 2010 all numerical methods began to migrate to graphics cards. Gradually, to single and double floating-point arithmetic, people added the so-called \textbf{half} (FP16) and \textbf{quarter} (FP8, INT8) floating points --- because with them it was faster to perform computations and they occupied less memory. In modern NVIDIA graphics cards (H100, B200) it is possible to perform about \textbf{a million operations} with such small objects in the same one nanosecond --- and this has allowed using such graphics cards for modern artificial intelligence systems. \subsubsection{GPU architecture} A modern GPU (for example, NVIDIA H100) is: \begin{itemize} \item \textbf{Streaming Multiprocessors (SM):} $\sim$132 multiprocessors, \item \textbf{CUDA cores:} $\sim$16\,896 cores for FP32 arithmetic, \item \textbf{Tensor cores:} $\sim$528 specialised cores for matrix operations (FP16, FP8, INT8), \item \textbf{Memory:} 80 GB HBM3 (High Bandwidth Memory) with a bandwidth of $\sim$3 TB/s. \end{itemize} \subsubsection{Threads and warps} A GPU works with \textbf{threads}. Thousands of threads execute in parallel, but they are organised into groups: \begin{itemize} \item \textbf{Warp:} A group of 32 threads that execute \textbf{synchronously} --- all 32 threads execute the same instruction at the same moment in time. \item \textbf{Block:} A group of warps (up to 1024 threads) that can exchange data through a shared \textbf{shared memory}. \item \textbf{Grid:} All blocks launched on the GPU. \end{itemize} \begin{orangebox}[Rule of GPU efficiency] For a GPU to work efficiently, the programs on each multiprocessor for each thread must \textbf{not diverge}. This means: \begin{itemize} \item All 32 threads in a warp must execute \textit{the same} instruction (without conditional jumps that split the threads --- ``warp divergence''), \item All threads must access \textit{sequential} memory addresses (coalesced memory access), \item There must be enough threads to ``hide'' the memory latencies. \end{itemize} If these conditions are not met --- the GPU works several times more slowly. \end{orangebox} \subsection{Parallel computing: Flynn's taxonomy} In 1966 Michael Flynn proposed a classification of parallel architectures by instruction and data streams: \begin{tabularx}{\textwidth}{l X l} \toprule \textbf{Type} & \textbf{Description} & \textbf{Example} \\ \midrule \textbf{SISD} & Single Instruction, Single Data. One instruction stream, one data stream. Classical sequential processor. & Old CPUs \\ \textbf{SIMD} & Single Instruction, Multiple Data. One instruction applied to many data. & Vector instructions (AVX), GPU \\ \textbf{MISD} & Multiple Instruction, Single Data. Different instructions on the same data. Rarely encountered. & Some specialised architectures \\ \textbf{MIMD} & Multiple Instruction, Multiple Data. Many processors, each executing its own instructions on its own data. & Modern multiprocessors, clusters \\ \bottomrule \end{tabularx} \subsubsection{Shared and distributed memory} In MIMD systems there are two approaches: \begin{itemize} \item \textbf{Shared memory:} All processors have access to a single memory. Example: a multi-core CPU. Easy to program (all see the same data), but difficult to scale (memory access conflicts). \item \textbf{Distributed memory:} Each processor has its own local memory, data exchange --- through a network. Example: clusters, supercomputers. More difficult to program (one must explicitly transfer data --- MPI), but it scales to millions of cores. \end{itemize} \subsection{TOP500: the most powerful supercomputers in the world} The TOP500 list is a ranking of the 500 most powerful supercomputers in the world, updated twice a year (June and November). Performance is measured in the LINPACK benchmark --- solving a large system of linear equations. \subsubsection{Modern leaders (2024--2025)} \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Name} & \textbf{Rank} & \textbf{Performance} & \textbf{Architecture} \\ \midrule Frontier & \#1 & $\sim$1.2 Eflop/s & AMD EPYC + AMD Instinct MI250X \\ Aurora & \#2 & $\sim$1 Eflop/s & Intel Xeon + Intel Max PVC \\ Eagle & \#3 & $\sim$561 Pflop/s & AMD EPYC + NVIDIA H100 \\ Fugaku & \#4 & $\sim$442 Pflop/s & Fujitsu A64FX (ARM) \\ LUMI & \#5 & $\sim$231 Pflop/s & AMD EPYC + AMD MI250X \\ \bottomrule \end{tabularx} \begin{tipbox}[Scale] 1 Eflop/s = $10^{18}$ operations per second. Frontier is $\sim$1.2 million billion operations per second. For comparison: your laptop is $\sim$0.1 Tflop/s = $10^{11}$ operations per second. The difference is 10 million times! \end{tipbox} \subsubsection{How a modern supercomputer is arranged} Typical architecture: \begin{enumerate} \item \textbf{Nodes:} $\sim$10\,000--100\,000 servers, each with 2--4 CPUs + 4--8 GPUs, \item \textbf{Network:} High-speed interconnect (InfiniBand, Slingshot) with a bandwidth of 200--400 Gbit/s per node, \item \textbf{Storage:} A parallel file system (Lustre, GPFS) with a bandwidth of $\sim$1--10 TB/s, \item \textbf{Cooling:} Liquid cooling (water or even two-phase), since the power is 20--50 MW. \end{enumerate} \subsection{Digitisers: a bridge between the analogue and digital worlds} But the processor --- after all, it must get data for computation from somewhere. And here \textbf{analog-to-digital converters} (ADC) enter the stage. They are determined by two parameters: \begin{itemize} \item \textbf{Accuracy (bit depth):} How many bits are used to encode one sample, \item \textbf{Speed (sampling rate):} How many samples per second. \end{itemize} \subsubsection{The trade-off: accuracy vs speed} There is a fundamental trade-off: the higher the accuracy, the lower the speed. \begin{tabularx}{\textwidth}{l l l X} \toprule \textbf{Bit depth} & \textbf{Speed} & \textbf{Minimum bit} & \textbf{Application} \\ \midrule 12 bit & $\sim$10 Gsamp/s & $\sim$500 $\mu$V & Oscilloscopes, fast imaging \\ 16 bit & $\sim$1 Gsamp/s & $\sim$10--50 $\mu$V & Audio, NMR, precise measurements \\ 24 bit & $\sim$4 Msamp/s & $\sim$100--500 nV & High-resolution audio, seismography \\ \bottomrule \end{tabularx} \begin{warningbox}[Why are nanovolts not measured directly?] 24 bits give a range of tens and hundreds of nanovolts. But nobody measures nanovolts directly --- radar interference from everything around is about tens of microvolts! That is, one usually measures \textit{currents} there, not volts --- or uses special shielded chambers. \end{warningbox} \subsection{The physics of switching: from the transistor to megavolts} In a modern processor everything works by turning transistors on and off at about 0.7 V. The speed of such switching is tens of gigahertz. But the currents there are very small (microamperes per transistor). \subsubsection{The trade-off: voltage vs speed} If the voltage is increased, the switching speed begins to fall: \begin{tabularx}{\textwidth}{l l X} \toprule \textbf{Voltage} & \textbf{Switching time} & \textbf{Technology} \\ \midrule 0.7 V & $\sim$0.05 ns (20 GHz) & Modern CMOS transistors \\ 40 V & $\sim$3--5 ns & Power MOSFETs (ordinary ratings) \\ 1000 V & $\sim$3--5 ns & Power MOSFETs (top ratings) \\ 4.5--5 kV & $\sim$20--30 ns & IGBT, high-voltage MOSFETs \\ $>$10 kV & $\sim$100 ns -- 1 $\mu$s & A single semiconductor can no longer cope \\ 1 MV & $\sim$100--500 $\mu$s & Tubes, semiconductor assemblies, multipliers \\ \bottomrule \end{tabularx} \textit{Although modern CMOS transistors switch at frequencies around 20 GHz, to execute one processor cycle one must have a margin of several such switchings, therefore the clock frequency of processors is in the range of about 3-4 GHz.} \subsubsection{The limit of voltage rise} Above $\sim$5 kV a single semiconductor can no longer cope --- one must switch to: \begin{itemize} \item \textbf{Good old tubes} (!) --- vacuum triodes, thyratrons, \item \textbf{Semiconductor assemblies} --- series connection of several transistors, \item \textbf{Voltage multipliers} --- cascade circuits with tubes or spark gaps. \end{itemize} Creating even 1 megavolt is not very difficult (Marx generator, Cockcroft--Walton cascades). Another matter is to turn this megavolt on and off in a nanosecond. That really requires hundreds of microseconds. \begin{successbox}[Fundamental limit] In fact $\sim 10^{12}$ V/s is the modern limit of voltage rise practically over the entire range of possible voltages. This is a physical limitation associated with: \begin{itemize} \item The speed of motion of charge carriers in semiconductors (limit $\sim 10^7$ cm/s), \item Parasitic capacitances and inductances, \item The speed of propagation of electromagnetic waves in a medium. \end{itemize} \end{successbox} \subsection{Connection with previous chapters} Let us see how everything we have covered is connected with the ``hardware'': \begin{itemize} \item \textbf{Linear algebra:} Matrix operations are the basis of all computations. GPUs are optimised precisely for them (Tensor Cores). \item \textbf{FFT:} The fast Fourier transform is $\mathcal{O}(N \log N)$ operations, and it is critical for NMR, signal processing, data compression. GPUs accelerate FFT by tens of times. \item \textbf{Iterative methods:} The conjugate gradient method, GMRES --- these are the basis for solving large sparse systems. They work on supercomputers with distributed memory (MPI + OpenMP + CUDA). \item \textbf{Machine learning:} Training neural networks is billions of matrix multiplications. Without GPUs and Tensor Cores modern models (GPT-4, AlphaFold) would be impossible. \item \textbf{Compressed sensing:} L1 minimisation is an iterative method that requires thousands of iterations. GPUs accelerate it 100--1000 times. \item \textbf{Digitisers:} ADC is the bridge between the analogue world (spectra, signals) and the digital one (processor). The accuracy of the ADC determines what mathematics we can apply. \end{itemize} \begin{successbox}[Main conclusion] Modern computing hardware is: \begin{itemize} \item \textbf{Memory hierarchy:} Registers $\to$ L1 $\to$ L2 $\to$ L3 $\to$ RAM $\to$ disk. Each step is 10--100 times slower, but 10--1000 times larger. \item \textbf{Parallelism:} SIMD (vector instructions), multi-core (CPU), massive parallelism (GPU), clusters (supercomputers). \item \textbf{Specialisation:} Tensor Cores for AI, FPGA for specific tasks, ASIC for mining. \item \textbf{Physical limitations:} Memory wall, power wall, speed of light --- all this limits the growth of performance. \end{itemize} For a chemist this means: \begin{itemize} \item Understanding where the ``bottleneck'' is (memory or computations), \item Being able to write code that efficiently uses the cache and GPU, \item Knowing when a problem is solved on a workstation, and when a supercomputer is needed. \end{itemize} \end{successbox} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Afterword: How to use this knowledge} Previously, after reading a similar book, when you encountered a real problem, there arose a need to implement a solution. For this one had to be quite good at one or several programming languages and understand how to use third-party software systems and packages. In the modern age of artificial intelligence development all this can be done with the support of AI systems --- and it will be significantly more effective. But for this one must understand a little the key applications of programming languages and the so-called \textbf{workflow} --- how program development happens. \subsection{Two worlds of programming languages} Although there are a huge number of programming languages, they can be divided into two fundamentally different classes: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Compiled} & \textbf{Interpreted / bytecode} \\ \midrule C, C++, CUDA, Fortran, Rust & Python, JavaScript, Java, C\#, R \\ \midrule Code is compiled into machine instructions for a specific processor & Code is interpreted on the fly or compiled into bytecode \\ \midrule Maximum performance, direct access to the ``hardware'' & Flexibility, cross-platform, fast development \\ \midrule \textbf{When to use:} & \textbf{When to use:} \\ \quad --- Heavy computations (DFT, MD, ML) & \quad --- Prototyping, data analysis \\ \quad --- GPU acceleration (CUDA) & \quad --- Visualisation, reporting \\ \quad --- Embedded systems & \quad --- Web interfaces, scripts \\ \bottomrule \end{tabularx} \subsubsection{Why interpreted languages are convenient} Interpreted languages allow achieving greater flexibility when moving from one platform to another. For example, we open the same web page on a desktop, on a smartphone and on a tablet --- and see the same content. This content is usually generated using the interpreted language JavaScript, and if we had to send our code for each of the processors every time, the server would have to have all variants of such code in advance. This is sometimes even impossible to foresee --- since even different Android smartphones may have different processor architectures. \subsubsection{Typical chemist's workflow} In order to quickly display some results, it is often easier to use interpreted programming languages (Python --- the absolute standard in science). But when the task is complex and requires many computational resources, it becomes necessary to use compiled programming languages. A typical workflow looks like this: \begin{enumerate} \item \textbf{Prototype in Python:} Quickly check the idea, visualise the data, understand whether the algorithm works. \item \textbf{Optimisation:} If the code works but is slow --- rewrite the ``bottlenecks'' in C++/CUDA or use ready-made libraries (NumPy, SciPy, PyTorch). \item \textbf{Production:} If the task is solved regularly --- package it, add tests, documentation. \end{enumerate} \subsection{Working with AI assistants} Nowadays practically any algorithm can be implemented by writing a correct prompt with the help of chats. But it is worth noting: each chat is almost like a person. And if you give it a task that is not fully thought through, it may also muddle things, or, not understanding, program something not quite right. \subsubsection{Golden rules of working with AI} \begin{enumerate} \item \textbf{Break the task into small blocks.} Do not ask ``write a program for solving the Schrödinger equation''. Ask: ``write a function that computes the Fock matrix for a given basis''. \item \textbf{Ask for tests.} For each block, ask the chat to write test verification programs, run such tests and make sure everything is in order. \item \textbf{Combine gradually.} Only after each block works, combine them to obtain the final result. \item \textbf{Demand explanations.} If the chat has written something you do not understand --- ask it to explain. If the explanation is unclear --- ask for something simpler. This is your script, you must understand every line. \item \textbf{Check on known examples.} If you are solving a system of equations --- check on a system for which you know the answer. If you are doing a Fourier transform --- check on a sinusoid with a known frequency. \end{enumerate} \begin{warningbox}[Typical mistakes] \begin{itemize} \item \textbf{Blind trust:} AI may write code that ``looks right'' but contains subtle errors (for example, incorrect normalisation, wrong array bounds, memory leaks). \item \textbf{Ignoring context:} AI does not know your specific task --- you must explain it in detail, with examples, with constraints. \item \textbf{Absence of tests:} Code without tests is code that will break at the most inopportune moment. \end{itemize} \end{warningbox} \subsection{Must-have tools for a chemist} Here is the minimum set of tools you should know: \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Language} & \textbf{Application} \\ \midrule \textbf{Python} & The absolute standard: data analysis, ML, visualisation, scripts \\ \textbf{R} & Statistical analysis, bioinformatics \\ \textbf{MATLAB} & Engineering computations, signal processing (an alternative to Python) \\ \textbf{C++/CUDA} & High-performance computing, GPU \\ \textbf{Bash/Shell} & Automation, working with clusters \\ \bottomrule \end{tabularx} \begin{tabularx}{\textwidth}{l X} \toprule \textbf{Library} & \textbf{Purpose} \\ \midrule \textbf{NumPy, SciPy} & Linear algebra, optimisation, integration \\ \textbf{Pandas} & Working with tabular data \\ \textbf{Matplotlib, Seaborn} & Visualisation \\ \textbf{PyTorch, TensorFlow} & Machine learning, neural networks \\ \textbf{RDKit} & Chemoinformatics, molecules \\ \textbf{ASE, PySCF} & Quantum chemistry \\ \textbf{OpenMM, GROMACS} & Molecular dynamics \\ \bottomrule \end{tabularx} \subsection{Final advice} Mathematics is not a set of formulas to be memorised. It is a \textbf{way of thinking}. It is the ability to see structure in chaos, to find patterns, to build models, to test hypotheses. When you encounter a real problem --- whether it is protein structure prediction, NMR spectrum analysis, synthesis optimisation or something completely new --- remember this script. Remember that: \begin{itemize} \item Any problem is either a system of equations, or an optimisation problem, or an approximation problem, \item For any problem there are proven mathematical methods, \item Modern tools (Python, GPU, AI) allow solving problems that were impossible 20 years ago, \item The main thing is to understand the essence of the problem, and the rest is technique. \end{itemize} We believe in you. You will succeed. And if someday this script helps you solve a problem that changes the world --- we will be incredibly proud. \begin{successbox}[Last word] Do not be afraid to make mistakes. Do not be afraid to ask. Do not be afraid to experiment. Mathematics is not an exam, it is an adventure. And you are at the beginning of the most interesting path. Good luck! \end{successbox} \vfill \begin{center} \textit{With love, \\ Papa and Mama} \\[0.5cm] \textit{September 2026} \end{center} \end{document}