Expand knowledge base

1 files changed, 158 insertions, 0 deletions

diff --git a/content/know/concept/pulay-mixing/index.pdc b/content/know/concept/pulay-mixing/index.pdc
new file mode 100644
index 0000000..9102c0e
--- /dev/null
+++ b/content/know/concept/pulay-mixing/index.pdc
@@ -0,0 +1,158 @@
+---
+title: "Pulay mixing"
+firstLetter: "P"
+publishDate: 2021-03-02
+categories:
+- Numerical methods
+
+date: 2021-03-02T19:11:51+01:00
+draft: false
+markup: pandoc
+---
+
+# Pulay mixing
+Some numerical problems are most easily solved *iteratively*,
+by generating a series $\rho_1$, $\rho_2$, etc.
+converging towards the desired solution $\rho_*$.
+**Pulay mixing**, also often called
+**direct inversion in the iterative subspace** (DIIS),
+is an effective method to speed up convergence,
+which also helps to avoid periodic divergences.
+
+The key concept it relies on is the **residual vector** $R_n$
+of the $n$th iteration, which in some way measures the error of the current $\rho_n$.
+Its exact definition varies,
+but is generally along the lines of the difference between
+the input of the iteration and the raw resulting output:
+
+$$\begin{aligned}
+	R_n
+	= R[\rho_n]
+	= \rho_n^\mathrm{new}[\rho_n] - \rho_n
+\end{aligned}$$
+
+It is not always clear what to do with $\rho_n^\mathrm{new}$.
+Directly using it as the next input ($\rho_{n+1} = \rho_n^\mathrm{new}$)
+often leads to oscillation,
+and linear mixing ($\rho_{n+1} = (1\!-\!f) \rho_n + f \rho_n^\mathrm{new}$)
+can take a very long time to converge properly.
+Pulay mixing offers an improvement.
+
+The idea is to construct the next iteration's input $\rho_{n+1}$
+as a linear combination of the previous inputs $\rho_1$, $\rho_2$, ..., $\rho_n$,
+such that it is as close as possible to the optimal $\rho_*$:
+
+$$\begin{aligned}
+	\boxed{
+		\rho_{n+1}
+		= \sum_{m = 1}^n \alpha_m \rho_m
+	}
+\end{aligned}$$
+
+To do so, we make two assumptions.
+Firstly, the current $\rho_n$ is already close to $\rho_*$,
+so that such a linear combination makes sense.
+Secondly, the iteration is linear,
+such that the raw output $\rho_{n+1}^\mathrm{new}$
+is also a linear combination with the *same coefficients*:
+
+$$\begin{aligned}
+	\rho_{n+1}^\mathrm{new}
+	= \sum_{m = 1}^n \alpha_m \rho_m^\mathrm{new}
+\end{aligned}$$
+
+We will return to these assumptions later.
+The point is that $R_{n+1}$ is also a linear combination:
+
+$$\begin{aligned}
+	R_{n+1}
+	= \rho_{n+1}^\mathrm{new} - \rho_{n+1}
+	= \sum_{m = 1}^n \alpha_m \rho_m^\mathrm{new} - \sum_{m = 1}^n \alpha_m \rho_m
+	= \sum_{m = 1}^n \alpha_m R_m
+\end{aligned}$$
+
+The goal is to choose the coefficients $\alpha_m$ such that
+the norm of the error $|R_{n+1}| \approx 0$,
+subject to the following constraint to preserve the normalization of $\rho_{n+1}$:
+
+$$\begin{aligned}
+	\sum_{m=1}^n \alpha_m = 1
+\end{aligned}$$
+
+We thus want to minimize the following quantity,
+where $\lambda$ is a [Lagrange multiplier](/know/concept/lagrange-multiplier/):
+
+$$\begin{aligned}
+	\braket{R_{n+1}}{R_{n+1}} + \lambda \sum_{m = 1}^n \alpha_m
+	= \sum_{m=1}^n \alpha_m \Big( \sum_{k=1}^n \alpha_k \braket{R_m}{R_k} + \lambda \Big)
+\end{aligned}$$
+
+By differentiating the right-hand side with respect to $\alpha_m$,
+we get a system of equations that we can write in matrix form,
+which is relatively cheap to solve numerically:
+
+$$\begin{aligned}
+	\begin{bmatrix}
+		\braket{R_1}{R_1} & \cdots & \braket{R_1}{R_n} & 1 \\
+		\vdots & \ddots & \vdots & \vdots \\
+		\braket{R_n}{R_1} & \cdots & \braket{R_n}{R_n} & 1 \\
+		1 & \cdots & 1 & 0
+	\end{bmatrix}
+	\begin{bmatrix}
+		\alpha_1 \\ \vdots \\ \alpha_n \\ \lambda
+	\end{bmatrix}
+	=
+	\begin{bmatrix}
+		0 \\ \vdots \\ 0 \\ 1
+	\end{bmatrix}
+\end{aligned}$$
+
+This method is very effective.
+However, in practice, the earlier inputs $\rho_1$, $\rho_2$, etc.
+are much further from $\rho_*$ than $\rho_n$,
+so usually only the most recent $N$ inputs $\rho_{n - N}$, ..., $\rho_n$ are used:
+
+$$\begin{aligned}
+	\rho_{n+1}
+	= \sum_{m = N}^n \alpha_m \rho_m
+\end{aligned}$$
+
+You might be confused by the absence of all $\rho_m^\mathrm{new}$
+in the creation of $\rho_{n+1}$, as if the iteration's outputs are being ignored.
+This is due to the first assumption,
+which states that $\rho_n^\mathrm{new}$ are $\rho_n$ are already similar,
+such that they are interchangeable.
+
+Speaking of which, about those assumptions:
+while they will clearly become more accurate as $\rho_n$ approaches $\rho_*$,
+they might be very dubious in the beginning.
+A consequence of this is that the early iterations might get "trapped"
+in a suboptimal subspace spanned by $\rho_1$, $\rho_2$, etc.
+To say it another way, we would be varying $n$ coefficients $\alpha_m$
+to try to optimize a $D$-dimensional $\rho_{n+1}$,
+where in general $D \gg n$, at least in the beginning.
+
+There is an easy fix to this problem:
+add a small amount of the raw residual $R_m$
+to "nudge" $\rho_{n+1}$ towards the right subspace,
+where $\beta \in [0,1]$ is a tunable parameter:
+
+$$\begin{aligned}
+	\boxed{
+		\rho_{n+1}
+		= \sum_{m = N}^n \alpha_m (\rho_m + \beta R_m)
+	}
+\end{aligned}$$
+
+In other words, we end up introducing a small amount of the raw outputs $\rho_m^\mathrm{new}$,
+while still giving more weight to iterations with smaller residuals.
+
+Pulay mixing is very effective:
+it can accelerate convergence by up to one order of magnitude!
+
+
+
+## References
+1.  P. Pulay,
+    [Convergence acceleration of iterative sequences. The case of SCF iteration](https://doi.org/10.1016/0009-2614(80)80396-4),
+    1980, Elsevier.