The Mathematics of DIIS

Suppose that we have a set of trial vectors $\left\{ {\mathbf p}^i \right\}$ which have been generated during the iterative solution of a problem. Now let us form a set of ``residual'' vectors defined as

$$\Delta {\mathbf p}^i = {\mathbf p}^{i+1} - {\mathbf p}^i.$$ (1)

The DIIS method assumes that a good approximation to the final solution p^f can be obtained as a linear combination of the previous guess vectors

$${\mathbf p} = \sum_i^m c_i {\mathbf p}^i,$$ (2)

where m is the number of previous vectors (in practice, only the most recent few vectors are used). The coefficients c_i are obtained by requiring that the associated residual vector

$$\Delta {\mathbf p} = \sum_i^m c_i \left( \Delta {\mathbf p}^i \right)$$ (3)

approximates the zero vector in a least-squares sense. Furthermore, the coefficients are required to add to one,

$$\sum_i^m c_i = 1.$$ (4)

The motivation for the latter requirement can be seen as follows. Each of our trial solutions pⁱ can be written as the exact solution plus an error term, p^f + eⁱ. Then, the DIIS approximate solution is given by

$$\sum_i^m c_i \left( {\mathbf p}^f + {\mathbf e}^i \right)$$ p = (5)

$${\mathbf p}^f \sum_i^m c_i + \sum_i^m c_i {\mathbf e}^i.$$ =

Hence, we wish to minimize the actual error, which is the second term in the equation above (of course, in practice, we don't know eⁱ, only $\Delta {\mathbf p}^i$); doing so would make the second term vanish, leaving only the first term. For p = p^f, we must have $\sum_i^m c_i = 1$.

Thus, we wish to minimize the norm of the residuum vector

$$\langle \Delta {\mathbf p} \vert \Delta {\mathbf p} \rangle =... ...angle \Delta {\mathbf p}^i \vert \Delta {\mathbf p}^j \rangle,$$ (6)

subject to the constraint (4). These requirements can be satisfied by minimizing the following function with Lagrangian multiplier $\lambda$

$${\cal L} = {\mathbf c}^{\dag } {\mathbf B} {\mathbf c} - \lambda \left( 1 - \sum_i^m c_i \right),$$ (7)

where B is the matrix of overlaps

$$B_{ij} = \langle \Delta {\mathbf p}^i \vert {\Delta \mathbf p}^j \rangle.$$ (8)

We can minimize ${\cal L}$ with respect to a coefficient c_k to obtain (assuming real quantities)

$$\frac{\partial {\cal L}}{\partial c_k} = 0 \sum_j^m c_j B_{kj} + \sum_i^m c_i B_{ik} - \lambda$$ = (9)

$$2 \sum_i^m c_i B_{ki} - \lambda.$$ =

We can absorb the factor of 2 into $\lambda$ to obtain the following matrix equation, which is eq. (6) of Pulay [3]:

$$\left( \begin{array}{ccccc} B_{11} & B_{12} & \cdots & B_{1... ...} 0 \\ 0 \\ \cdots \\ 0 \\ -1 \\ \end{array}\right)$$ (10)