A density matrix approach to the convergence of the self-consistent field iteration

Abstract

In this paper, we present a local convergence analysis of the self-consistent field (SCF) iteration using the density matrix as the state of a fixed-point iteration. Sufficient and almost necessary conditions for local convergence are formulated in terms of the spectral radius of the Jacobian of a fixed-point map. The relationship between convergence and certain properties of the problem is explored by deriving upper bounds expressed in terms of higher gaps. This gives more information regarding how the gaps between eigenvalues of the problem affect the convergence, and hence these bounds are more insightful on the convergence behaviour than standard convergence results. We also provide a detailed analysis to describe the difference between the bounds and the exact convergence factor for an illustrative example. Finally we present numerical examples and compare the exact value of the convergence factor with the observed behaviour of SCF, along with our new bounds and the characterization using the higher gaps. We provide heuristic convergence factor estimates in situations where the bounds fail to well capture the convergence.

keywords

self-consistent field iteration, convergence analysis, nonlinear eigenvalue problems, eigenvector nonlinearity, electronic structure calculations, iterative methods

Parikshit Upadhyaya

Lindstedtsvägen 25,

Department of Mathematics,

SeRC - Swedish e-Science research center,

Royal Institute of Technology, SE-11428 Stockholm, Sweden

Elias Jarlebring

Lindstedtsvägen 25,

Department of Mathematics,

SeRC - Swedish e-Science research center,

Royal Institute of Technology, SE-11428 Stockholm, Sweden

Emanuel H. Rubensson

Division of Scientific Computing, Department of Information Technology,

Uppsala University, Box 337, SE-75105 Uppsala, Sweden

1 Introduction

Let $A:M\to M$ , where $M\subset\mathbb{C}^{n\times n}$ denotes the set of Hermitian matrices. In this work we consider the associated nonlinear eigenvalue problem consisting of determining $(X_{1},\Lambda_{1})\in\mathbb{C}^{n\times p}\times\mathbb{R}^{p\times p}$ such that $(X_{1},\Lambda_{1})$ is an invariant pair of $A(X_{1}X_{1}^{H})$ , i.e.,


$\displaystyle A(X_{1}X_{1}^{H})X_{1}$	$\displaystyle=$	$\displaystyle X_{1}\Lambda_{1},$	(1a)
$\displaystyle X_{1}^{H}X_{1}$	$\displaystyle=$	$\displaystyle I,$	(1b)

where $\Lambda_{1}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{p})$ and $\lambda_{1},\ldots,\lambda_{n}$ are the eigenvalues of $A(X_{1}X_{1}^{H})$ , numbered in ascending order. This is one of the fundamental computational challenges in quantum chemistry and related fields (Hartree-Fock and Kohn-Sham density functional theory).¹¹ 1 In these settings, the columns of $X_{1}$ contain the basis-expansion coefficients of the molecular orbitals. However, this problem also arises as a trace ratio maximization problem in linear discriminant analysis for dimension reduction. See [14],[24] and [25] for more on this application. The self-consistent field (SCF) iteration consists of computing iterates satisfying the linear eigenvalue problem

A(V_{k}V_{k}^{H})V_{k+1}=V_{k+1}S_{k+1},

(2)

where $V_{k+1}^{H}V_{k+1}=I$ and $S_{k+1}\in\mathbb{R}^{p\times p}$ is diagonal. In the case of convergence, $V_{k}\to X_{1}$ and $S_{k}\to\Lambda_{1}$ . In this paper we provide a local convergence analysis of this algorithm. SCF is rarely used on its own as a solution method and most of the state-of-the-art procedures are based on its enhancements and improvements, for example, Pulay’s DIIS (Direct Inversion in the Iterative Subspace) acceleration in [15]. See also standard references [22],[5] and further literature discussion below.

There is an extensive amount of literature on the convergence of the SCF iteration and its variants. We mention some main approaches to the convergence theory, without an ambition of a complete description. A number of recent works are based on the optimization viewpoint, e.g., [3, 9, 10, 11]. This is natural, since the problem in (1) often stems from the first order optimality condition of an energy minimization problem, as in [20, Section 2.1]. In particular, the Roothaan algorithm with level-shifting and damping are studied in [3]. This analysis was used as a basis for the gradient analysis in [9], which provided explicit estimates of the convergence rate for the algorithms applied to the Hartree-Fock equations. The convergence of the DIIS acceleration scheme has been studied separately in [16].

Various approaches are based on measuring the subspace angle and other using chordal norms, e.g., [10, 2], leading to local convergence as well as a global convergence analysis. In contrast to these approaches, we use a density matrix based analysis and derive bounds involving higher gaps (as we explain below). The analysis in [11] provides precise conditions for local convergence (and some global convergence conditions), under the assumption that $A$ only depends on the diagonal of the density matrix $X_{1}X_{1}^{H}$ , which in many discretization settings corresponds to the charge density. Not all problems are nonlinear only in the diagonal of the density matrix, as e.g., the example in Section 4.2. The work in [23] also provides a convergence analysis, mostly based on a non-zero temperature filter function approach. We note also that a precise local convergence criterion for the classical version was presented in [21], not involving a density matrix analysis.

A model of interacting bosons which has received considerable attention is the Gross-Pitaevskii equation. This corresponds to (1) with $p=1$ . Convergence results for the SCF iteration for this case can be found for example in [1].

Our convergence analysis is focused on establishing a precise characterization of the convergence factor as well as natural upper bounds. We provide an exact formula for the convergence factor, which turns out to be the spectral radius of a matrix (the Jacobian of the fixed-point map). Using this exact formula, we derive upper bounds which can be phrased in terms of higher gaps (as we define later in Definition 1.2) and the action of a linear operator on the outer products of eigenvectors. We also provide an example where the convergence cannot be characterized based on the first gap alone, which illustrates the importance of taking into account the higher gaps in the convergence analysis. This should be viewed in contrast to the analysis in [10, 2, 23], which is primarily focused on the first gap.

We will use a formulation of the SCF iteration in terms of the density matrix. In our context, a density matrix is defined by

P_{k}:=V_{k}V_{k}^{H}\in M.

Given $P_{k}$ we can compute $A(P_{k})=A(V_{k}V_{k}^{H})$ from which we can compute $V_{k+1}$ and in principle construct $P_{k+1}=V_{k+1}V_{k+1}^{H}$ . Hence, the iteration (2) is equivalent to a fixed point iteration in $P_{k}$ . We will refer to this fixed point map as $\Psi$ , i.e.,

P_{k+1}=\Psi(P_{k}).

(3)

Although our conclusions hold for general problems, we restrict our analysis to the case where the operator $A$ has the form

A(P)=A_{0}+\mathcal{L}(P),

(4)

where $A_{0}\in M$ and $\mathcal{L}:\mathbb{C}^{n\times n}\to\mathbb{C}^{n\times n}$ is a complex linear operator, i.e., $\mathcal{L}(zA)=z\mathcal{L}(A)$ for all $z\in\mathbb{C}$ . The density matrix formulation in (3) is the starting point of several linear scaling variants of the SCF-algorithm that avoid explicit construction of $V_{k+1}$ [4]. Before we proceed to the next section, we will introduce some necessary notation and definitions.

Let $X,\Lambda$ correspond to a complete eigenvalue decomposition of $A$ evaluated in a solution to (1), i.e.,

A(X_{1}X_{1}^{H})X=X\Lambda,

where

\Lambda=\begin{bmatrix}\Lambda_{1}&\\ &\Lambda_{2}\end{bmatrix}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{p},\lambda_{p+1},\ldots,\lambda_{n})

(5)

and

X=[X_{1}\;\;X_{2}]=[x_{1}\;\;\ldots\;\;x_{n}].

Definition 1.1 (Gap).

The smallest distance between the diagonal elements of $\Lambda_{1}$ and the diagonal elements of $\Lambda_{2}$ is defined as as the gap and denoted $\delta$ .

Due to the numbering of eigenvalues, the gap is given by

\delta=\min_{i\leq p,\;j\geq p+1}|\lambda_{i}-\lambda_{j}|=\lambda_{p+1}-\lambda_{p}.

Throughout this paper, we assume that $\delta\neq 0$ . Otherwise the decomposition (5) is not unique and violates the common uniform well-posedness hypothesis of [3].

Definition 1.2 (Higher gap).

The $j$ -th smallest distance between the diagonal elements of $\Lambda_{1}$ and $\Lambda_{2}$ is denoted by $\delta_{j}$ .

Note that $\delta_{1}=\delta$ . As an example, the second gap is given by

\delta_{2}=\underset{(i,j)\neq(p,p+1)}{\min_{i\leq p,\;j\geq p+1}}|\lambda_{i}-\lambda_{j}|=\min(\lambda_{p+2}-\lambda_{p},\lambda_{p+1}-\lambda_{p-1}).

Definition 1.3 (Lower triangular vectorization).

Let $m=n(n+1)/2$ . The operator $\mathit{{\operatorname{vech}}}\colon M\to\mathbb{C}^{m}$ is defined as

{\operatorname{vech}}(W)=\begin{bmatrix}w_{1,1}&\cdots&w_{n,1}&w_{2,2}&\cdots&w_{n,2}&w_{3,3}&\cdots&w_{n,n}\end{bmatrix}^{T}.

This is the vectorization operator adapted for hermitian matrices, and returns the vectorizaton of the lower triangular part. Similarly, we define the inverse operator ${\operatorname{vech}}^{-1}\colon\mathbb{C}^{m}\to M$ which maps any vector $v\in\mathbb{C}^{m}$ to a corresponding $W\in M$ . The relation between ${\operatorname{vec}}$ and ${\operatorname{vech}}$ is given by

{\operatorname{vech}}(W)=T{\operatorname{vec}}(W),

(6)

where $T\in\mathbb{R}^{m\times n^{2}}$ . The matrix $T$ is in general non-unique (as discussed in [12] and [6]). In this paper, we will specifically use (6) with

T=\operatorname{diag}(I_{n},\begin{bmatrix}0&I_{n-1}\end{bmatrix},\begin{bmatrix}0&0&I_{n-2}\end{bmatrix},\cdots,1).

2 Convergence characterization

2.1 Main theory

The following characterization involves the matrix consisting of reciprocal gaps, which we denote $R\in\mathbb{R}^{n\times n}$ , and is given by

R_{i,j}=\begin{cases}\frac{1}{\lambda_{j}(0)-\lambda_{i}(0)},\;\;&\textrm{ if }i\leq p\textrm{ and }j>p\\ \frac{1}{\lambda_{i}(0)-\lambda_{j}(0)},\;\;&\textrm{ if }i>p\textrm{ and }j\leq p\\ 0&\textrm{ otherwise}.\end{cases}

(7)

where $\lambda_{1}(t),\ldots,\lambda_{n}(t)$ are eigenvalues of a parameter dependent matrix $B(t)$ .²² 2 Note that the matrix $R$ has appeared with different names in other papers before. For example, in [21], it is referred to as the ”density perturbation”. In [11], it is called ”first divided difference matrix”. The matrix $R$ is symmetric with the following structure:

R=\begin{bmatrix}0&R_{p}^{T}\\ R_{p}&0\end{bmatrix},

where $R_{p}\in\mathbb{R}^{(n-p)\times p}$ . We need the following perturbation result whose variants exist in quantum mechanical perturbation theory, for example in chapter 15.III of [13].

Lemma 2.1 (Density matrix derivatives).

Consider a matrix-valued function $B$ depending on a complex parameter such that $B(t)=B_{0}+B_{1}t$ , where $B_{0},B_{1}$ are Hermitian. Let $X,\Lambda$ correspond to a parameter dependent diagonalization for a sufficiently small neighborhood of $t_{0}=0$ , i.e.,

X(t)\Lambda(t)X(t)^{H}=B(t),

with $X(t)^{H}X(t)=I$ and $\Lambda(0)=\operatorname{diag}(\lambda_{1}(0),\ldots,\lambda_{n}(0))$ , where $\lambda_{1}(0)<\cdots<\lambda_{p}(0)<\lambda_{p+1}(0)<\cdots<\lambda_{n}(0)$ . Let $X(t)\in\mathbb{C}^{n\times n}$ be decomposed as $X(t)=[X_{1}(t),X_{2}(t)]$ , where $X_{1}(t)\in\mathbb{C}^{n\times p}$ and $P(t):=X_{1}(t)X_{1}(t)^{H}$ . Then,

{\operatorname{vec}}(P^{\prime}(0))=-(\overline{X(0)}\otimes X(0))D(X(0)^{T}\otimes X(0)^{H}){\operatorname{vec}}(B_{1})

where

D:=\operatorname{diag}({\operatorname{vec}}(R)).

(8)

Proof.

Let us consider the domain $\mathcal{D}=\left(-\infty,\lambda_{mid}\right)\cup\left(\lambda_{mid},\infty\right)$ , where $\lambda_{mid}=\frac{\lambda_{p}(0)+\lambda_{p+1}(0)}{2}$ and define the step function $h:\mathcal{D}\to\mathbb{R}$ , which acts as a filter for the $p$ smallest eigenvalues of $B(0)=B_{0}$ :

h(x)=\begin{cases}1,\quad x<\lambda_{mid}\\ 0,\quad x>\lambda_{mid}\end{cases}.

(9)

We can generalize $h:M\to M$ as a matrix function (in the sense of [7]) and rewrite the density matrix function as

P(t)=X_{1}(t)X_{1}(t)^{H}=X(t)\begin{bmatrix}I_{p}&0\\ 0&0\\ \end{bmatrix}X(t)^{H}=X(t)h(\Lambda(t))X(t)^{H}=h(B(t)),

(10)

where we assume that $t$ lies in a sufficiently small neighbourhood of $t_{0}=0$ such that $\lambda_{p}(t)<\lambda_{mid}<\lambda_{p+1}(t)$ . Then,

	$\displaystyle P^{\prime}(0)=$	$\displaystyle\lim\limits_{\epsilon\to 0}$	$\displaystyle\frac{h(B(\epsilon))-h(B(0))}{\epsilon}$
	$\displaystyle=$	$\displaystyle\lim\limits_{\epsilon\to 0,\|\|E\|\|\to 0}$	$\displaystyle\frac{X(0)\Big(h(\Lambda(0)+X(0)^{H}EX(0))-h(\Lambda(0))\Big)X(0)^{H}}{\epsilon},$		(11)

where

E=\epsilon B_{1}.

(12)

If we denote by $L(F,G)$ the Fréchet derivative of $h$ evaluated at $F$ and applied to $G$ , then from (2.1), we get

$\displaystyle P^{\prime}(0)$	$\displaystyle=$	$\displaystyle\lim\limits_{\epsilon\to 0,\|\|E\|\|\to 0}$	$\displaystyle\frac{X(0)\Big(L\big(\Lambda(0),X(0)^{H}EX(0)\big)+O(\|\|E\|\|^{2})\Big)X(0)^{H}}{\epsilon}$
	$\displaystyle=$	$\displaystyle\lim\limits_{\epsilon\to 0,\|\|E\|\|\to 0}$	$\displaystyle\frac{X(0)L\big(\Lambda(0),X(0)^{H}EX(0)\big)X(0)^{H}}{\epsilon}$
	$\displaystyle=$		$\displaystyle X(0)L\big(\Lambda(0),X(0)^{H}B_{1}X(0)\big)X(0)^{H}$	(13)

The last equation is a consequence of (12) and the fact that $O(||E||^{2})/\epsilon$ goes to zero as $\epsilon$ goes to zero and $L$ being linear in the second argument. From the Daleckiĭ-Kreĭn theorem[7, Theorem 3.11],[8], we have

L(\Lambda(0),X(0)^{H}B_{1}X(0))=U\circ(X(0)^{H}B_{1}X(0)),

(14)

where $U$ denotes the matrix of divided differences

U:=\begin{cases}\frac{h(\lambda_{i}(0))-h(\lambda_{j}(0))}{\lambda_{i}(0)-\lambda_{j}(0)}&i\neq j\\ h^{\prime}(\lambda_{i}(0))&i=j.\end{cases}

By definition of $h$ in (9) and from (7), we see that $U$ reduces to $-R$ . Vectorizing (2.1) and using (14) leads us to

	$\displaystyle{\operatorname{vec}}(P^{\prime}(0))$	$\displaystyle=-{\operatorname{vec}}\Big(X(0)\big(R\circ(X(0)^{H}B_{1}X(0))\big)X(0)^{H}\Big)$		(15)
		$\displaystyle=-(\overline{X(0)}\otimes X(0))D(X(0)^{T}\otimes X(0)^{H}){\operatorname{vec}}(B_{1}),$

which is a result of repeated application of the matrix product vectorization identity ${\operatorname{vec}}(LMN)=(N^{T}\otimes L){\operatorname{vec}}(M)$ . ∎

In practice, the computation of the density matrix as in (3) can be done by using (10). Moreover, in SCF iterations, variations of the problem can be solved, where an approximation of the step function is used. A common choice for this approximation is the Fermi-Dirac distribution: $f_{\mu,\beta}(t)=\frac{1}{1+e^{\beta(t-\mu)}}$ , where the parameter $\mu$ is usually selected such that ${\mathrm{Tr}}(P(0))=p$ . The function $f_{\mu,\beta}$ tends to the step function in the limit $\beta\to 0$ . Note that Lemma 2.1 can be generalized for such an approximation of the density matrix $P_{f}$ as follows:

{\operatorname{vec}}(P^{\prime}_{f}(0))=(\overline{X(0)}\otimes X(0))D_{f}(X(0)^{T}\otimes X(0)^{H}){\operatorname{vec}}(B^{\prime}(0)),

where $D_{f}=\operatorname{diag}({\operatorname{vec}}(R_{f}))$ and

R_{f}=\begin{cases}\frac{f(\lambda_{i}(0))-f(\lambda_{j}(0))}{\lambda_{i}(0)-\lambda_{j}(0)}&i\neq j\\ f^{\prime}(\lambda_{i}(0))&i=j\end{cases}.

Theorem 1 (Density matrix local convergence).

Let $P_{*}=X_{1*}X_{1*}^{H}\in\mathbb{C}^{n\times n}$ be a fixed point of $\Psi$ , i.e., $P_{*}=\Psi(P_{*})$ . Then, the SCF iteration satisfies

{\operatorname{vech}}(P_{k+1}-P_{*})={\operatorname{vech}}(\Psi(P_{k})-P_{*})=J_{P}{\operatorname{vech}}(P_{k}-P_{*})+O(\|{\operatorname{vech}}(P_{k}-P_{*})\|^{2})

where

J_{P}=-T(\overline{X_{*}}\otimes X_{*})D({X_{*}}^{T}\otimes{X_{*}}^{H})L^{\prime}

(16)

and $L^{\prime}\in\mathbb{C}^{n^{2}\times m}$ is defined by

L^{\prime}=\left({\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{1})),\ldots,{\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{m})))\right),

(17)

$X_{*}=[X_{1*},X_{2*}]$ and $e_{1},e_{2},\ldots,e_{m}$ are the first $m$ columns of the identity matrix $I_{n}$ .

Proof.

Applying the operator ${\operatorname{vech}}(\cdot)$ to our iteration, we get

	$\displaystyle{\operatorname{vech}}(P_{k+1})$	$\displaystyle={\operatorname{vech}}(\Psi(P_{k}))$
		$\displaystyle={\operatorname{vech}}({\operatorname{vec}}^{-1}({\operatorname{vec}}(\Psi({\operatorname{vech}}^{-1}({\operatorname{vech}}(P_{k})))))).$

Hence, the fixed point iteration can be re-written as

{\operatorname{vech}}(P_{k+1})=f({\operatorname{vech}}(P_{k})),

where $f:\mathbb{R}^{m}\to\mathbb{R}^{m}$ ,

f(v)={\operatorname{vech}}({\operatorname{vec}}^{-1}({\operatorname{vec}}(\Psi({\operatorname{vech}}^{-1}(v)))))\quad\forall v\in\mathbb{R}^{m}.

A Taylor expansion around the fixed-point ${\operatorname{vech}}(P_{*})$ gives us,

{\operatorname{vech}}(P_{k+1}-P_{*})=J_{P}{\operatorname{vech}}(P_{k}-P_{*})+O(\|{\operatorname{vech}}(P_{k}-P_{*})\|^{2}),

where $J_{P}$ is the Jacobian of $f$ evaluated in ${\operatorname{vech}}(P_{*})$ . The $j$ -th column of $J_{P}$ is given by

	$\displaystyle J_{P}(:,j)$	$\displaystyle=\lim_{\epsilon\to 0}\frac{f({\operatorname{vech}}(P_{})+\epsilon e_{j})-f({\operatorname{vech}}(P_{}))}{\epsilon}$		(18)
		$\displaystyle=\lim_{\epsilon\to 0}\frac{{\operatorname{vech}}({\operatorname{vec}}^{-1}({\operatorname{vec}}(\Psi(P_{}+\epsilon{\operatorname{vech}}^{-1}(e_{j}))-\Psi(P_{}))))}{\epsilon}.$		(18)

By using linearity of the vectorization operators and of $\mathcal{L}$ , we can now invoke Lemma 2.1 with $B(\alpha)=A+\mathcal{L}(P_{*}+\alpha{\operatorname{vech}}^{-1}(e_{j}))$ to get

\lim_{\epsilon\to 0}\frac{{\operatorname{vec}}(\Psi(P_{*}+\epsilon{\operatorname{vech}}^{-1}(e_{j}))-\Psi(P_{*}))}{\epsilon}=\\ (\overline{X_{*}}\otimes X_{*})D({X_{*}}^{T}\otimes{X_{*}}^{H}){\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{j}))).

Using this and equation (18), we have

	$\displaystyle J_{P}(:,j)$	$\displaystyle={\operatorname{vech}}({\operatorname{vec}}^{-1}((\overline{X_{}}\otimes{X_{}})D({X_{}}^{T}\otimes{X_{}}^{H}){\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{j})))))$		(19)
		$\displaystyle=T(\overline{X_{}}\otimes{X_{}})D({X_{}}^{T}\otimes{X_{}}^{H}){\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{j}))).$		(19)

Due to the fact that ${\operatorname{vec}}(\mathcal{L}({\operatorname{vech}}^{-1}(e_{j})))$ is the only component in (19) depending on $j$ , we obtain (16) by factorizing the matrix $T(\overline{X_{*}}\otimes{X_{*}})D({X_{*}}^{T}\otimes{X_{*}}^{H})$ . ∎

3 Convergence factor bounds and their interpretation

3.1 Spectral-norm bounds

Since (3) is a nonlinear fixed-point map and the Jacobian is evaluated at a fixed point in Theorem 1, the convergence factor is

c=\rho(J_{P})

where $\rho(J_{P})$ denotes the spectral radius of $J_{P}$ . Moreover, a sufficient and almost necessary condition for local convergence is $c<1$ . Due to the fact that the spectral radius is smaller than any operator norm, we have in particular for the spectral norm:

c\leq||J_{P}||_{2}=:c_{2}.

3.1.1 Naive bounds

Now note that $\overline{X_{*}}\otimes X_{*}$ and $X_{*}^{T}\otimes X_{*}^{H}$ are orthogonal matrices, and that $D$ defined by (8) is a diagonal matrix whose largest element is the reciprocal gap, such that

\|D\|_{2}=\max_{i}|d_{i,i}|=\frac{1}{\delta}.

By using this and the Cauchy-Schwartz inequality we obtain a straight-forward upper bound for $c$ ,

c\leq||T||_{2}||\overline{X}\otimes{X}||_{2}||D||_{2}||({X}^{T}\otimes{X}^{H})||_{2}||L^{\prime}||_{2}\leq\frac{||L^{\prime}||_{2}}{\delta}:=c_{\rm naive},

(20)

where we dropped subscript $*$ in the eigenvector matrix $X_{*}$ for notational convenience. We can conclude from (20) that a small gap implies a larger value of the upper bound $c_{naive}$ , indicating slow convergence. This is consistent with the well-known fact that problems with a small gap are more difficult to solve using the SCF iteration which is concluded in several convergence analysis works, e.g. [23]. Note that the bound (20) does not depend on the gap alone but also on $||L^{\prime}||_{2}$ , which can be large and difficult to analyze. The matrix $L^{\prime}$ depends on the action of the operator $\mathcal{L}$ , which leads us to the pursuit of other bounds which may quantify this dependence in a way that is easier to interpret.

3.1.2 Cycled permutation

Different bounds can be derived by using the fact that the spectral radius does not change when we reverse the order of multiplication of matrices, i.e., $\rho(AB)=\rho(BA)$ . Therefore, from the definition of $J_{P}$ and

c=\rho(J_{P})=\rho(T(\overline{X}\otimes X)D({X}^{T}\otimes{X}^{H})L^{\prime}),

we obtain variants based on cyclic permutation


$\displaystyle c$	$\displaystyle=$	$\displaystyle\rho((\overline{X}\otimes X)D({X}^{T}\otimes{X}^{H})L^{\prime}T)$	(21a)
	$\displaystyle=$	$\displaystyle\rho(D({X}^{T}\otimes{X}^{H})L^{\prime}T(\overline{X}\otimes X)).$	(21b)

Both equations in (21) lead to the bound

c\leq\|D({X}^{T}\otimes{X}^{H})L^{\prime}T\|_{2}=:c_{2,a}.

The cyclic permutation can be continued such that


$\displaystyle c$	$\displaystyle=$	$\displaystyle\rho(({X}^{T}\otimes{X}^{H})L^{\prime}T(\overline{X}\otimes X)D)$	(22a)
	$\displaystyle=$	$\displaystyle\rho(L^{\prime}T(\overline{X}\otimes X)D({X}^{T}\otimes{X}^{H})).$	(22b)

Equation (22) leads to the bound

c\leq\|L^{\prime}T(\overline{X}\otimes X)D\|_{2}=:c_{2,b}.

(23)

In the following we need the symmetrization operator formally defined as

S(X):=\sum_{j=1}^{m}{\operatorname{vech}}(X)_{j}{\operatorname{vech}^{-1}}(e_{j})

(24)

or equivalently $S(L+D+R)=L+D+L^{T}$ , where $L+D+R$ is the decomposition into the lower triangular, diagonal and upper triangular matrices. Using (24), we have the identity

	$\displaystyle L^{\prime}{\operatorname{vech}}(X)$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{m}{\operatorname{vec}}(\mathcal{L}({\operatorname{vech}^{-1}}(e_{j}))){\operatorname{vech}}(X)_{j}$		(25)
		$\displaystyle=$	$\displaystyle{\operatorname{vec}}(\mathcal{L}(S(X))),$		(25)

due to the definition of $L^{\prime}$ in (17) and the linearity of $\mathcal{L}$ . The columns of the matrix in (23) can be expressed as

	$\displaystyle(L^{\prime}T(\overline{X}\otimes X)D)_{:,j}$	$\displaystyle=$	$\displaystyle L^{\prime}T{\operatorname{vec}}(x_{\ell}x_{m}^{H})d_{j,j}$
		$\displaystyle=$	$\displaystyle L^{\prime}{\operatorname{vech}}(x_{\ell}x_{m}^{H})d_{j,j}$

where $j=n(m-1)+\ell$ and where we used (6) in the last step. Hence, from (25) we have

(L^{\prime}T(\overline{X}\otimes X)D)_{:,j}={\operatorname{vec}}(\mathcal{L}(S(x_{\ell}x_{m}^{H})))d_{j,j}.

(26)

The columns of the matrix inside the norm in (23) is given by (26). This formula can be interpreted as follows. We clearly see that the action of the linear operator applied to the outer products of eigenvectors $\mathcal{L}(S(x_{\ell}x_{m}^{H}))$ has significance. The weighting with $d_{j,j}$ implies that only pairs of eigenvectors of different occupancy are relevant (in this bound). This quantity is further described in the following section.

3.2 Higher gaps

We begin by decomposing the matrix obtained by cycled permutation in (22) as follows,

L^{\prime}T(\overline{X}\otimes X)D({X}^{T}\otimes{X}^{H})=\\ L^{\prime}T(\overline{X}\otimes X)(D-D_{1}-\cdots-D_{2q})({X}^{T}\otimes{X}^{H})+L^{\prime}T(\overline{X}\otimes X)(D_{1}+\cdots+D_{2q})({X}^{T}\otimes{X}^{H})

(27)

where $D_{j}$ are rank one diagonal matrices (and we take $2q$ terms for symmetry reasons). We have $q=p(n-p)$ unique gaps and they occur twice each on the diagonal of $D$ . Hence, this decomposition reveals the dependence of the convergence on eigenvalue gaps other than just the smallest gap $\delta$ . In the following theorem we quantify this dependence along with the dependence on the outer products of the eigenvectors $\mathcal{L}(S(x_{\ell}x_{m}^{H}))$ as introduced in the previous subsection.

In the formulation of the theorem, we use the set $\Omega_{q}\subset[1,n]\times[1,n]$ which contains the indices of $R$ that have entries corresponding to the $q$ smallest gaps. As a result, $\Omega_{q}$ contains $2q$ elements.

Refer to caption — Figure 1: Schematic illustration of elements of $\Omega_{3}$ as indices of $R$ for the real-valued problem in subsection 4.1 with $n=7,p=3,\alpha=10.0$ .

In figure 1, we visualize the elements of $\Omega_{3}$ . The set $\Omega_{3}$ comprises the indices of the reciprocal gap matrix $R$ which correspond to these gaps, that is, $\Omega_{3}=\{(4,3),(3,4),(4,2),(2,4),(5,3),(3,5)\}$ .

Theorem 3.1 (Higher gaps).

The convergence factor of the SCF-iteration is bounded by

\rho(J_{p})\leq\frac{\|L^{\prime}\|_{2}}{\delta_{q+1}}+\sum_{(\ell,m)\in\Omega_{q}}\frac{1}{|\lambda_{\ell}-\lambda_{m}|}\|\mathcal{L}(S(x_{\ell}x_{m}^{H}))\|_{F}:=c_{gap,q}

(28)

for any $q\in[0,p(n-p)]$ , where $\delta_{p(n-p)+1}:=\infty$ and $\delta_{1}=\delta$ .

Proof.

For notational convenience, we express $D$ in terms of $R$ (defined in (7)), i.e.,

$\displaystyle D$	$\displaystyle=$	$\displaystyle\operatorname{diag}({\operatorname{vec}}(R))=\sum_{\ell,m}^{n}r_{\ell,m}{\operatorname{vec}}(e_{\ell}e_{m}^{T}){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}$
	$\displaystyle=$	$\displaystyle\sum_{(\ell,m)\not\in\Omega_{q}}r_{\ell,m}{\operatorname{vec}}(e_{\ell}e_{m}^{T}){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}+$
		$\displaystyle\sum_{(\ell,m)\in\Omega_{q}}r_{\ell,m}{\operatorname{vec}}(e_{\ell}e_{m}^{T}){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}.$

The idea of the proof is to take the last sum in this equation as $\sum_{j=1}^{2q}D_{j}$ with the decomposition in (27). By the triangle inequality, we have

c\leq\|L^{\prime}T(\overline{X}\otimes X)(D-D_{1}-\cdots-D_{q})\|+\sum_{j=1}^{2q}\|L^{\prime}T(\overline{X}\otimes X)D_{j}\|.

(29)

The first term in (29) is of the form used in the naive bound (20), except that the diagonal matrix is modified by setting the contribution corresponding to the $q$ first gaps to zero. We obtain directly the first term in (28),

\|L^{\prime}T(\overline{X}\otimes X)(D-D_{1}-\cdots-D_{2q})\|\leq\frac{\|L^{\prime}\|_{2}}{\delta_{q+1}}.

The second term in (29) is

\sum_{j=1}^{2q}\|L^{\prime}T(\overline{X}\otimes X)D_{j}\|=\sum_{(\ell,m)\in\Omega_{q}}r_{\ell,m}\|L^{\prime}T(\overline{X}\otimes X){\operatorname{vec}}(e_{\ell}e_{m}^{T}){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}\|.

This can be simplified with the identity (26) which implies that


$\displaystyle\\|L^{\prime}T(\overline{X}\otimes X){\operatorname{vec}}(e_{\ell}e_{m}^{T}){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}\\|$	$\displaystyle=$	$\displaystyle\\|{\operatorname{vec}}(\mathcal{L}(S(x_{\ell}x_{m}^{H}))){\operatorname{vec}}(e_{\ell}e_{m}^{T})^{T}\\|$	(30a)
	$\displaystyle=$	$\displaystyle\\|{\operatorname{vec}}(\mathcal{L}(S(x_{\ell}x_{m}^{H})))\\|$	(30b)
	$\displaystyle=$	$\displaystyle\\|\mathcal{L}(S(x_{\ell}x_{m}^{H}))\\|_{F}.$	(30c)

The last two equalities follow from the fact that the spectral norm of a matrix with one non-zero column, is the two-norm of that column vector, which is the Frobenius norm of the corresponding matrix. The proof is concluded by combining (29) with (30) and noting that $r_{\ell,m}=\frac{1}{|\lambda_{\ell}-\lambda_{m}|}$ .

∎

Theorem 3.1 should be further interpreted as follows. The parameter $q$ is free and the theorem therefore provides us with a family of bounds parameterized by $q$ . For example, $q=0$ gives us $\rho(J_{P})\leq\frac{||L^{\prime}||_{2}}{\delta_{1}}$ , which is the naive bound from (20). For $q=1$ , we have

\rho(J_{p})\leq\frac{||L^{\prime}||_{2}}{\delta_{2}}+\frac{||\mathcal{L}(S(x_{p}x_{p+1}^{H}))||_{F}+||\mathcal{L}(S(x_{p+1}x_{p}^{H}))||_{F}}{\delta_{1}}.

By induction, $q=k$ gives us a bound that is a function of the $k+1$ smallest gaps and the norm of the action of $\mathcal{L}$ on the outer products of eigenvector-pairs corresponding to the gap indices of the $k$ smallest gaps.

3.3 Illustrative example

In order to illustrate the insight provided by Theorem 3.1, we provide an example showing how in certain situations, using a bound with a higher value of $q$ provides us tighter upper bounds. Consider the following problem parameterized by $\epsilon$ .

\left(\begin{bmatrix}0&\epsilon&0\\ \epsilon&1+d&\epsilon\\ 0&\epsilon&10\\ \end{bmatrix}+\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&100\\ \end{bmatrix}\circ x_{1}x_{1}^{H}\right)X=X\Lambda.

(31)

Let $d=0.16$ . In the notation in equation (4), we have,

A_{0}=\begin{bmatrix}0&\epsilon&0\\ \epsilon&1+d&\epsilon\\ 0&\epsilon&10\\ \end{bmatrix}\quad\mathcal{L}(P)=L\circ P=\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&100\\ \end{bmatrix}\circ P.

First note that for $\epsilon=0$ , the solution is given by

\Lambda(0)=\begin{bmatrix}\lambda_{1}&0&0\\ 0&\lambda_{2}&0\\ 0&0&\lambda_{3}\\ \end{bmatrix}=\begin{bmatrix}1&0&0\\ 0&1+d&0\\ 0&0&10\\ \end{bmatrix},\quad X(0)=\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&1\\ \end{bmatrix}.

Varying $\epsilon$ , solving the resulting problem instances and plotting the convergence factor and its bounds gives us figure 2. In figure 3, we plot the norm of $L^{\prime}$ along with the norm of action of $\mathcal{L}$ on the outer products of the eigenvectors ( $x_{1}x_{2}^{H}$ and $x_{1}x_{3}^{H}$ ) with $\epsilon$ . In figure 4, we visualize how the gaps $\delta_{1}$ and $\delta_{2}$ vary with $\epsilon$ . Note that we have $||L^{\prime}||_{2}=100$ , which is independent of $\epsilon$ . As can be seen from figures 2 and 4, $c_{naive}$ shows a direct inverse dependence on $\delta_{1}$ , which is expected from (20) because $||L^{\prime}||$ is constant.

We can see from figures 2 and 3 that $c_{naive}$ is not a good approximation to the convergence factor $c$ , in comparison to our bounds $c_{gap,2}$ , $c_{2}$ and $c_{2,a}$ . The increase in $c$ as we increase $\epsilon$ (which means slower convergence for larger values of $\epsilon$ ) is not captured by $c_{naive}$ or $c_{gap,1}$ as they are essentially constant for small $\epsilon$ . However, from figure 2, we see that the increase in $c$ coincides with the increase of the norm of action of $\mathcal{L}$ on the outer products of the eigenvectors ( $x_{1}x_{2}^{H}$ and $x_{1}x_{3}^{H}$ ), as seen in figure 3. This behaviour is better captured in the formula for the upper bound $c_{gap,2}$ ,

	$\displaystyle c_{gap,2}$	$\displaystyle=\frac{\|\|\mathcal{L}(S(x_{1}x_{2}^{H}))\|\|_{F}+\|\|\mathcal{L}(S(x_{2}x_{1}^{H}))\|\|_{F}}{\delta_{1}}$
		$\displaystyle+\frac{\|\|\mathcal{L}(S(x_{1}x_{3}^{H}))\|\|_{F}+\|\|\mathcal{L}(S(x_{3}x_{1}^{H}))\|\|_{F}}{\delta_{2}}.$

Although the bounds $c_{gap,2}$ , $c_{2}$ and $c_{2,a}$ are better approximations of $c$ as compared to $c_{naive}$ , there is still a discrepancy in the slopes in figure 2, and the rate of increase of $c$ is faster than that of the bounds. We now provide a more detailed analysis. First note that differentiating (31) with respect to $\epsilon$ , and setting $\epsilon=0$ , we obtain

X^{\prime}(0)=\begin{bmatrix}0&\frac{1}{\lambda_{2}-\lambda_{1}}&0\\ \frac{1}{\lambda_{1}-\lambda_{2}}&0&\frac{1}{\lambda_{3}-\lambda_{2}}\\ 0&\frac{1}{\lambda_{2}-\lambda_{3}}&0\end{bmatrix},\quad\Lambda^{\prime}(0)=0,\quad D^{\prime}(0)=0.

(32)

Let $J(\epsilon)$ denote the parameter dependent Jacobian evaluated at the solution,

J(\epsilon)=-T\left(X(\epsilon)\otimes X(\epsilon)\right)D(\epsilon)\left(X(\epsilon)^{H}\otimes X(\epsilon)^{H}\right)L^{\prime}.

Differentiating w.r.t $\epsilon$ , setting $\epsilon=0$ , and using $X(0)=I$ ,

\begin{split}J^{\prime}(0)=-T\bigg(&\left(X^{\prime}(0)\otimes I+I\otimes X^{\prime}(0)\right)D(0)+D^{\prime}(0)\\ &-D(0)\left(X^{\prime}(0)\otimes I+I\otimes X^{\prime}(0)\right)\bigg)L^{\prime}.\\ \end{split}

(33)

Using the formulae from (32) and substituting into (33), we get,

J^{\prime}(0)=\frac{1}{{(\lambda_{2}-\lambda_{1})}^{2}}\begin{bmatrix}-e_{2}&0&0&e_{2}&0&0\end{bmatrix}.

From the structure of $J^{\prime}(0)$ and the fact that all eigenvalues of $J^{\prime}(0)$ are zero, we have,

	$\displaystyle\rho(J^{\prime}(0))$	$\displaystyle=$	$\displaystyle 0,$		(34)
	$\displaystyle\|\|J^{\prime}(0)\|\|_{2}$	$\displaystyle=$	$\displaystyle\frac{1}{{(\lambda_{2}-\lambda_{1})}^{2}}\quad.$

This allows us to carry out a Taylor series analysis for $\rho\left(J(\epsilon)\right)$ and $||J(\epsilon)||_{2}$ around 0,

c=\rho(J(\epsilon))=\rho(J(0)+\epsilon J^{\prime}(0)+\mathcal{O}(\epsilon^{2}))\approx\epsilon\rho(J^{\prime}(0))+\mathcal{O}(\epsilon^{2})=\mathcal{O}(\epsilon^{2}).

(35)

Similarly,

c_{2}=||J(\epsilon)||_{2}=\epsilon||J^{\prime}(0)||_{2}+\mathcal{O}(\epsilon^{2})=\frac{\epsilon}{{(\lambda_{2}-\lambda_{1})}^{2}}+\mathcal{O}(\epsilon^{2})=\mathcal{O}(\epsilon).

(36)

Hence, from (35) and (36), we expect $c$ to vary at a rate that is an order of magnitude faster than $c_{2}$ for very small values of $\epsilon$ , which is exactly what we observe in figure 2. We clearly see from (34),(35) and (36) that this is because $J^{\prime}(0)$ has zero eigenvalues (and hence zero spectral radius) but non-zero norm. This illustrates how the two-norm based bounds can overestimate $c$ .

4 Numerical examples

4.1 Discrete Laplacian example

In this subsection, we apply our theory to a minor variation of the problem type discussed in [10, Section 5]. In the context of this paper, this translates to

A_{0}=\begin{bmatrix}2/h^{2}&-1/h^{2}+i/2h&0&\ldots&0\\ -1/h^{2}-i/2h&2/h^{2}&-1/h^{2}+i/2h&\ddots&\vdots\\ 0&-1/h^{2}-i/2h&2/h^{2}&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&-1+i/2h\\ 0&\ldots&0&-1/h^{2}-i/2h&2/h^{2}\\ \end{bmatrix},

which is the discretized 1D differential operator: $\dfrac{\partial^{2}}{\partial x^{2}}+i\dfrac{\partial}{\partial x}$ . In a PDE setting, this would correspond to a diffusion term added with a complex convection term discretized with a central difference scheme with grid spacing $h$ . We also have

\mathcal{L}(P)=\alpha Diag\left(Re(A_{0})^{-1}diag\left(P\right)\right).

Note that $\mathcal{L}(\cdot)$ depends only on the diagonal of $P$ .

Figure 5(a) shows us that the predicted convergence rate $c$ agrees perfectly with SCF convergence history. The norm based bounds $c_{2}$ and $c_{naive}$ are slightly worse than the exact rate $c$ . As expected, $c_{naive}$ is the least accurate upper bound, but $c_{gap,2}$ is only slightly better. As seen from figure 5(b) the gaps between the eigenvalues are not very well separated, that is, the higher gaps are not much larger than $\delta$ . More precisely, $\delta_{3}$ is not much larger than $\delta_{1}$ and $\delta_{2}$ . This explains why $c_{gap,2}$ is not a good approximation to the norm based bounds for this problem. Figure 5(d) shows a linear increase in the value of the convergence factor and the upper bounds with $\alpha$ , which is expected from the linear dependence of $L^{\prime}$ on $\alpha$ and equation (21). From figure 5(c), we also see that convergence becomes faster(that is $c$ decreases) with increase in problem size for a constant value of $\alpha$ and $p$ . To make a comparison of our upper bounds with the convergence factor derived from [10, Theorem 4.2], we change the problem by setting

A_{0}=\begin{bmatrix}2/h^{2}&-1/h^{2}&0&\ldots&0\\ -1/h^{2}&2/h^{2}&-1/h^{2}&\ddots&\vdots\\ 0&-1/h^{2}&2/h^{2}&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&-1/h^{2}\\ 0&\ldots&0&-1/h^{2}&2/h^{2}\\ \end{bmatrix}

since the analysis in that paper is presented for real-valued problems. The operator $\mathcal{L}(\cdot)$ is the same as before. We plot the upper bound $c_{Liu}=\frac{2\alpha\sqrt{n}||A_{0}^{-1}||_{2}}{\delta_{1}}$ along with the other upper bounds for this modified real-valued problem in figure 6.

Figure 6 suggests to us that the upper bounds discussed in this paper are improvements over the upper bound in [10].

4.2 Water molecule example

In this example, we apply the SCF iteration to a problem that originates from the modelling of a water molecule system. The discretization involves a restricted Hartree-Fock approximation with a set of $n=13$ basis functions and $p=5$ . For any nonlinear eigenvalue problem that results from a Hartree-Fock approximation, we have

A(X_{1}X_{1}^{H})=H_{core}+2G(R^{-1}X_{1}X_{1}^{H}{R^{-1}}^{H})

where $G(\cdot)$ is a linear operator and $H_{core}$ is a sum of two matrices that correspond to terms for kinetic energy and the nuclear-electron interaction energy. In our context, $\mathcal{L}(P)=2G(R^{-1}P{R^{-1}}^{H})$ . Here, $R$ is a lower triangular matrix that results from a cholesky decomposition of the ”overlap matrix”. The overlap matrix is hermitian and obtained by computing integrals of products of basis functions, as explained in [17, Section 2.4]. For the purpose of reproducibility, we provide the coordinates of the nuclei of the Oxygen and Hydrogen atoms in the following table. Note that all data is in atomic units.

Atom	Charge( $\mathbf{e}$ )	x( $\mathbf{a_{0}}$ )	y( $\mathbf{a_{0}}$ )
O	8.0	0.0	0.0
H	1.0	-1.809	0.0
H	1.0	0.453549	1.751221

The computation was performed using Ergo[18, 19], which is a software package for large-scale SCF calculations. The standard Gaussian basis set 3-21G was used and the starting guess for the density matrix was projected from the calculation using a smaller STO-3G basis set. We plot the SCF convergence history and the exact convergence factor $c$ . We have not plotted the other upper bounds that we derived because they overestimate $c$ by a large margin. Instead, based on the theory in Theorem 3.1 we investigate the bound that neglects certain terms such that

\widetilde{c_{2}}=||T(\overline{X}\otimes X)D_{1}({X}^{T}\otimes{X}^{H})L^{\prime}||_{2}

which is the spectral norm of a 2-rank approximation of the Jacobian $J_{P}$ (taking into account only the entries that contain $\delta_{1}$ in $D$ ). As we can see from figure 7, the observed behaviour of SCF convergence agrees with that predicted by the exact value of the spectral radius, $c$ .

5 Conclusions and outlook

The SCF algorithm is an important algorithm in many fields. We have provided a new convergence characterization for the algorithm using a density matrix based analysis of a fixed point map. The upper bounds derived for the spectral radius of the Jacobian of the fixed point map illustrate how the convergence depends on the different problem parameters and physical properties. In particular, Theorem 3.1 provides a mathematical footing for studying how the gaps interact with the outer products of eigenvectors to affect the convergence properties. This is a quantification of Stanton’s observation in [21, Section IV], where he points out that typically, divergence in SCF calculations is not due to a single very low energy excitation, but to the interaction of several moderately low excitations. The discussion of the illustrative example in section 3.3 explains how when the Hessian has zero spectral radius but non-zero norm, an upper bound based on the interaction of higher gaps is needed to give a more accurate picture of the convergence behaviour. Finally, the application of our upper bounds to practical problems in subsections 4.1 and 4.2 reveal that our bounds are slightly better approximations to the convergence factor than the bounds that exist in previous literature.

References

[1] Z. Bai, D. Lu and B. Vandereycken, Robust Rayleigh quotient minimization and nonlinear eigenvalue problems, preprint, URL http://web.cs.ucdavis.edu/~bai/publications/bailuvandereycken18.pdf.
[2] Y. Cai, L.-H. Zhang, Z. Bai and R.-C. Li, On an eigenvector-dependent nonlinear eigenvalue problem, preprint, URL http://web.cs.ucdavis.edu/~bai/publications/caizhangbaili17.pdf.
[3] E. Cancès and C. L. Bris, On the convergence of SCF algorithms for the Hartree-Fock equations, M2AN, Math. Model. Numer. Anal., 34 (2000), 749–774, URL https://doi.org/10.1051/m2an:2000102.
[4] D.R.Bowler and T.Mizayaki, $\mathcal{O}$ ( $n$ ) methods in electronic structure calculations, Rep. Prog. Phys., 75 (2012), 036503, URL http://iopscience.iop.org/article/10.1088/0034-4885/75/3/036503.
[5] T. Helgaker, P. Jorgensen and J. Olsen, Molecular Electronic-Structure Theory, John Wiley and Sons, 2000.
[6] H. V. Henderson and S. R. Searle, Vec and vech operators for matrices, with some uses in jacobians and multivariate statistics, Can J. Stat., 7 (1979), 65–81, URL http://dx.doi.org/10.2307/3315017.
[7] N. Higham, Functions of Matrices, Society for Industrial and Applied Mathematics, 2008, URL https://epubs.siam.org/doi/abs/10.1137/1.9780898717778.
[8] J.L.Daleckiĭ and G.Kreĭn, Integration and differentiation of functions of hermitian operators and applications to the theory of perturbations, Amer. Math. Soc. Trans.: Series 2, 47 (1965), 1–30.
[9] A. Levitt, Convergence of gradient-based algorithms for the Hartree-Fock equations, ESAIM: Math. Model. Numer. Anal., 46 (2012), 1321–1336, URL https://doi.org/10.1051/m2an/2012008.
[10] X. Liu, X. Wang, Z. Wen and Y. Yuan, On the convergence of the self-consistent field iteration in Kohn-Sham density functional theory, SIAM J. Matrix Anal. Appl., 35 (2014), 546–558, URL https://doi.org/10.1137/130911032.
[11] X. Liu, Z. Wen, X. Wang, M. Ulbrich and Y. Yuan, On the analysis of the discretized Kohn-Sham density functional theory, SIAM J. Numer. Anal., 53 (2015), 1758–1785, URL https://doi.org/10.1137/140957962.
[12] C. E. McCulloch, Symmetric matrix derivatives with applications, J. Amer. Stat. Assoc., 77 (1982), 679–682, URL https://doi.org/10.2307/2287736.
[13] A. Messiah, Quantum Mechanics, Dover Publications, 1999.
[14] T. T. Ngo, M. Bellalij and Y. Saad, The trace ratio optimization problem for dimensionality reduction, SIAM J. Matrix Anal. Appl., 31 (2010), 2950–2971, URL https://doi.org/10.1137/090776603.
[15] P. Pulay, Convergence acceleration of iterative sequences. the case of scf iteration, Chem. Phys. Lett., 73 (1979), 393–398, URL https://doi.org/10.1016/0009-2614(80)80396-4.
[16] T. Rohwedder and R. Schneider, An analysis for the DIIS acceleration method used in quantum chemistry calculations, J. Math. Chem., 49 (2011), 1889–1914.
[17] E. Rudberg, Quantum Chemistry for Large Scale Systems, PhD thesis, Royal Institute of Technology, 2007, URL http://urn.kb.se/resolve?urn=urn:nbn:se:kth:diva-4561.
[18] E. Rudberg, E. H. Rubensson and P. Sałek, Kohn-Sham density functional theory electronic structure calculations with linearly scaling computational time and memory usage, J. Chem. Theory Comput., 7 (2011), 340–350, URL https://doi.org/10.1021/ct100611z.
[19] E. Rudberg, E. H. Rubensson, P. Sałek and A. Kruchinina, Ergo: An open-source program for linear-scaling electronic structure calculations, SoftwareX, 7 (2018), 107 – 111, URL https://doi.org/10.1016/j.softx.2018.03.005.
[20] Y. Saad, J. T. Chelikowsky and S. M. Shontz, Numerical methods for electronic structure calculations of materials, SIAM Rev., 52 (2010), 3–54, URL "https://epubs.siam.org/doi/pdf/10.1137/060651653".
[21] R. E. Stanton, Intrinsic convergence in closed-shell SCF calculations. a general criterion, J. Chem. Phys., 75 (1981), 5416–5422, URL https://doi.org/10.1063/1.441942.
[22] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications, 1996.
[23] C. Yang, W. Gao and J. C. Meza, On the convergence of the self-consistent field iteration for a class of nonlinear eigenvalue problems, SIAM J. Matrix Anal. Appl., 30 (2009), 1773–1788.
[24] L.-H. Zhang, L.-Z. Liao and M. K. Ng, Fast algorithms for the generalized foley–sammon discriminant analysis, SIAM J. Matrix Anal. Appl., 31 (2010), 1584–1605, URL https://doi.org/10.1137/080720863.
[25] L.-H. Zhang, W. H. Yang and L.-Z. Liao, A note on the trace quotient problem, Opt. Lett., 8 (2014), 1637–1645, URL https://doi.org/10.1007/s11590-013-0680-z.