---
title: "Pulay mixing"
sort_title: "Pulay mixing"
date: 2021-03-02
categories:
- Numerical methods
layout: "concept"
---

Some numerical problems are most easily solved *iteratively*,
by generating a series of "solutions" $$f_1$$, $$f_2$$, etc.
converging towards the true $$f_\infty$$.
**Pulay mixing**, also often called
**direct inversion in the iterative subspace** (DIIS),
can speed up the convergence for some types of problems,
and also helps to avoid periodic divergences.

The key concept it relies on is the **residual vector** $$R_n$$
of the $$n$$th iteration, which measures the error of the current $$f_n$$.
Its exact definition can vary,
but it is generally the difference between
the input $$f_n$$ of the $$n$$th iteration
and the raw resulting output $$f_n^\mathrm{new}$$:

$$\begin{aligned}
	R_n
	\equiv R[f_n]
	\equiv f_n^\mathrm{new}[f_n] - f_n
\end{aligned}$$

It is not always clear what to do with $$f_n^\mathrm{new}$$.
Directly using it as the next input ($$f_{n+1} = f_n^\mathrm{new}$$)
often leads to oscillation,
and linear mixing ($$f_{n+1} = (1\!-\!c) f_n + c f_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 $$f_{n+1}$$
as a linear combination of the previous inputs $$f_1$$, $$f_2$$, ..., $$f_n$$,
such that it is as close as possible to the optimal $$f_\infty$$:

$$\begin{aligned}
	\boxed{
		f_{n+1}
		= \sum_{m = 1}^n \alpha_m f_m
	}
\end{aligned}$$

To do so, we make two assumptions.
First, that the current $$f_n$$ is already close to $$f_\infty$$,
so such a linear combination makes sense.
Second, that the iteration is linear,
so the raw output $$f_{n+1}^\mathrm{new}$$
is also a linear combination *with the same coefficients*:

$$\begin{aligned}
	f_{n+1}^\mathrm{new}
	= \sum_{m = 1}^n \alpha_m f_m^\mathrm{new}
\end{aligned}$$

We will revisit these assumptions later.
The point is that $$R_{n+1}$$ can now also be written
as a linear combination of old residuals $$R_m$$:

$$\begin{aligned}
	R_{n+1}
	= f_{n+1}^\mathrm{new} - f_{n+1}
	= \sum_{m = 1}^n \alpha_m f_m^\mathrm{new} - \sum_{m = 1}^n \alpha_m f_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 $$f_{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}
	\inprod{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 \inprod{R_m}{R_k} + \lambda \Big)
\end{aligned}$$

By differentiating the right-hand side with respect to $$\alpha_m^*$$
and demanding that the result is zero,
we get a cheap-to-solve system of equations, in matrix form:

$$\begin{aligned}
	\begin{bmatrix}
		\inprod{R_1}{R_1} & \cdots & \inprod{R_1}{R_n} & 1 \\
		\vdots & \ddots & \vdots & \vdots \\
		\inprod{R_n}{R_1} & \cdots & \inprod{R_n}{R_n} & 1 \\
		1 & \cdots & 1 & 0
	\end{bmatrix}
	\cdot
	\begin{bmatrix}
		\alpha_1 \\ \vdots \\ \alpha_n \\ \lambda
	\end{bmatrix}
	=
	\begin{bmatrix}
		0 \\ \vdots \\ 0 \\ 1
	\end{bmatrix}
\end{aligned}$$

From this, we can also see that the Lagrange multiplier
$$\lambda = - \inprod{R_{n+1}}{R_{n+1}}$$,
where $$R_{n+1}$$ is the *predicted* residual of the next iteration,
subject to the two assumptions.
This fact makes $$\lambda$$ a useful measure of convergence.

In practice, the earlier inputs $$f_1$$, $$f_2$$, etc.
are much further from $$f_\infty$$ than $$f_n$$,
so usually only the most recent $$N\!+\!1$$ inputs $$f_{n - N}, ..., f_n$$ are used.
This also keeps the matrix small:

$$\begin{aligned}
	f_{n+1}
	= \sum_{m = n-N}^n \alpha_m f_m
\end{aligned}$$

You might be confused by the absence of any $$f_m^\mathrm{new}$$
in the creation of $$f_{n+1}$$, as if the iteration's outputs are being ignored.
This is due to the first assumption,
which states that $$f_n^\mathrm{new}$$ and $$f_n$$ are already similar,
such that they are basically interchangeable.

Although those assumptions will clearly become more realistic as $$f_n \to f_\infty$$,
they might be very dubious at first.
Consequently, the early iterations may get "trapped"
in a suboptimal subspace spanned by $$f_1$$, $$f_2$$, etc.
Think of it like this:
we would be varying up to $$n$$ coefficients $$\alpha_m$$
to try to optimize a $$D$$-dimensional $$f_{n+1}$$, where usually $$D \gg n$$.
It is almost impossible to find a decent optimum in this way!

This problem is easy to fix,
by mixing in a small amount of the raw residuals $$R_m$$
to "nudge" $$f_{n+1}$$ towards the right subspace,
where $$\beta \in [0,1]$$ is a tunable parameter:

$$\begin{aligned}
	\boxed{
		f_{n+1}
		= \sum_{m = N}^n \alpha_m (f_m + \beta R_m)
	}
\end{aligned}$$

In this way, the raw outputs $$f_m^\mathrm{new}$$ are (rightfully) included via $$R_m$$,
but we still give more weight to iterations with smaller residuals.

Pulay mixing is very effective for certain types of problems,
e.g. density functional theory,
where it can accelerate convergence by up to two orders 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.