
==== Front
J Sci Comput
J Sci Comput
Journal of Scientific Computing
0885-7474
1573-7691
Springer US New York

2612
10.1007/s10915-024-02612-3
Article
Riemannian Newton Methods for Energy Minimization Problems of Kohn–Sham Type
http://orcid.org/0000-0002-4161-6704
Altmann R. robert.altmann@ovgu.de

1
Peterseim D. 2
Stykel T. 2
1 https://ror.org/00ggpsq73 grid.5807.a 0000 0001 1018 4307 Institute of Analysis and Numerics, Otto von Guericke University Magdeburg, Universitätsplatz 2, 39106 Magdeburg, Germany
2 https://ror.org/03p14d497 grid.7307.3 0000 0001 2108 9006 Institute of Mathematics and Centre for Advanced Analytics and Predictive Sciences (CAAPS), University of Augsburg, Universitätsstr. 12a, 86159 Augsburg, Germany
13 8 2024
13 8 2024
2024
101 1 67 8 2023
12 6 2024
21 6 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
This paper is devoted to the numerical solution of constrained energy minimization problems arising in computational physics and chemistry such as the Gross–Pitaevskii and Kohn–Sham models. In particular, we introduce Riemannian Newton methods on the infinite-dimensional Stiefel and Grassmann manifolds. We study the geometry of these two manifolds, its impact on the Newton algorithms, and present expressions of the Riemannian Hessians in the infinite-dimensional setting, which are suitable for variational spatial discretizations. A series of numerical experiments illustrates the performance of the methods and demonstrates their supremacy compared to other well-established schemes such as the self-consistent field iteration and gradient descent schemes.

Keywords

Riemannian optimization
Stiefel manifold
Grassmann manifold
Newton method
Kohn–Sham model
Gross–Pitaevskii eigenvalue problem
Mathematics Subject Classification

65K10
65N25
81Q10
http://dx.doi.org/10.13039/100019180 HORIZON EUROPE European Research Council 865751 Peterseim D. issue-copyright-statement© Springer Science+Business Media, LLC, part of Springer Nature 2024
==== Body
pmcIntroduction

The Kohn–Sham model [31, 37, 38] is a prototypical example of a constrained energy minimization problem stated on the infinite-dimensional Stiefel manifold. This means that the sought-after minimizer is a p-frame of L2-orthonormal functions. Another well-known example is the Gross–Pitaevskii model for Bose–Einstein condensates of ultracold bosonic gases [39, 43]. Here, the special case p=1 is of interest, where we seek a single (minimizing) function on the unit sphere in L2, representing a unit mass constraint. Since these two applications are relevant for different communities, numerical methods are mostly considered separately. One aim of this paper is to give a unified approach to solving energy minimization problems. More precisely, we introduce Riemannian Newton methods for minimizing energy functionals of Kohn–Sham type, which also includes the Gross–Pitaevskii model. In general, the here considered PDE problems require a special treatment in terms of sparsity and dimension-independent methods, which is not part of existing general purpose optimization packages.

The numerical solution of the Gross–Pitaevskii model has been studied extensively in recent years. The most common numerical techniques are iterative methods based on Riemannian (conjugate) gradient descent methods or discretized Riemannian gradient flows in various metrics [10, 11, 20, 22, 26, 29, 36, 44, 53]. A conceptually different approach is the J-method [5, 33] with its inimitable sensitivity with regard to spectral shifts, allowing remarkable speed-ups in a Rayleigh quotient iteration manner. Reformulating the minimization problem as an eigenvalue problem with eigenvector nonlinearity—also known as nonlinear eigenvector problem—the self-consistent field iteration (SCF) can be employed; see [14, 23, 45]. This method involves the solution of a linear eigenvalue problem in each step and is strongly connected to the Newton method [28, 34]. Considering the extended nonlinear system including the normalization constraint also allows a direct application of Newton or Newton-type methods [12, 13, 24]. For an extended review on numerical methods for the Gross–Pitaevskii model, we refer to [28].

Most of the above approaches (with appropriate adjustments) have been applied to the Kohn–Sham model as well. This includes the direct constrained minimization algorithm [3, 47, 52] and the energy-adaptive gradient descent method [8]—both based on Riemannian optimization—as well as gradient flow schemes [21, 32]. Moreover, the SCF algorithm with different types of mixing is very popular in the computational chemistry community; see, e.g., [9, 15, 18, 19, 40]. For a discretized and simplified Kohn–Sham model (without the external potential and the exchange-correlation energy), global convergence and local second-order convergence of an inexact Riemannian Newton method on the Grassmann manifold has been shown in [54]. An overview of existing software packages for density functional theory problems can be found in [35].

In this paper, the point of origin is an energy functional defined on the infinite-dimensional Stiefel manifold, which we introduce in Sect. 2. For a better understanding, we recall definitions and properties of the Stiefel manifold and corresponding retractions in Sect. 3. Moreover, we provide formulae for the Riemannian gradient and the Riemannian Hessian which are needed for the Newton iteration. Since the considered energy functional is invariant under orthogonal matrices, we also discuss the infinite-dimensional Grassmann manifold and examine a connection of its tangent space to a certain subspace of the tangent space of the Stiefel manifold. The resulting Newton algorithms are then subject of Sect. 4. In particular, we present an inexact Riemannian Newton method on the Grassmann manifold. In Sect. 5, we consider the two mentioned examples of the Gross–Pitaevskii and the Kohn–Sham model in more detail. For both applications, we derive the formulae including a spatial discretization and illustrate the supremacy of the inexact Newton approach compared to well-established methods such as the SCF iteration and gradient descent schemes.

Notation

The sets of p×p real symmetric and skew-symmetric matrices are denoted by Ssym(p) and Sskew(p), respectively. For M∈Rp×p, we write symM=12(M+MT) for the symmetric part, and trM denotes the trace of M. Further, Ip and 0p denote the p×p identity and zero matrices, respectively. The expression diag(M) defines the column vector consisting of the diagonal elements of M∈Rn×n and Diag(v) denotes the diagonal matrix with components of the vector v∈Rn on the diagonal.

The Energy Functional and Nonlinear Eigenvector Problems

For a given spatial domain Ω⊆Rd, d≤3, we consider the Hilbert spaces L2(Ω) and V~⊆H1(Ω). For p≥1, we further define the Hilbert spaces V=V~p and H=[L2(Ω)]p of p-frames. Throughout this paper, we assume that V is dense in H and that V⊆H⊆V∗ form a Gelfand triple, where V∗ denotes the dual space of V.

For v=(v1,⋯,vp),w=(w1,⋯,wp)∈H, we define the dot productv·w=∑j=1pvjwj.

On the pivot space H, we further introduce an outer product1 〚v,w〛H=(v1,w1)L2(Ω)⋯(v1,wp)L2(Ω)⋮⋱⋮(vp,w1)L2(Ω)⋯(vp,wp)L2(Ω)∈Rp×p

and an inner product2 (v,w)H=∑j=1p(vj,wj)L2(Ω)=tr〚v,w〛H.

The inner product (2) induces the norm ‖v‖H=(v,v)H on H. The canonical identification I:V→V∗ is defined by⟨Iv,w⟩=(v,w)Hfor allv,w∈V,

where ⟨·,·⟩ denotes the duality pairing on V∗×V. This identification operator can also be written as I=j∗∘iH∘j with the trivial embedding j:V→H (the injective identity operator), the Riesz isomorphism iH:H→H∗, which reads iH(u)=(u,·)H, and the adjoint operator j∗:H∗→V∗ satisfying j∗(f)=f∘j for all f∈H∗. Since all these operators act componentwisely, we have I(vΛ)=I(v)Λ for all v∈V and Λ∈Rp×p. Moreover, since V is a dense subspace of H, so is j(V). Hence, j∗ is injective, and as the composition of injective operators, I is also injective. As a result, I has a left inverse J:V∗→V such that JIv=v for all v∈V.

Energy and Applications

For a p-frame ϕ∈V, we consider the energy functional3 E(ϕ)=12∫Ωtr((∇ϕ(x))T∇ϕ(x))dx+∫Ωϑ(x)ρ(ϕ(x))dx+12∫ΩΓ(ρ(ϕ(x)))dx

with an external potential ϑ, the density function ρ(ϕ)=ϕ·ϕ, and a smooth nonlinearity Γ(ρ). Our aim is to minimize this energy functional on the infinite-dimensional Stiefel manifold of index p given by4 St(p,V)={ϕ∈V:〚ϕ,ϕ〛H=Ip}.

In other words, we are interested in solving the constrained minimization problem5 minϕ∈St(p,V)E(ϕ).

A state of lowest energy is called the ground state. Such states play an important role in quantum-mechanical models as they represent a most stable configuration of atoms and molecules. These models include two famous applications in computational physics and chemistry.

Example 1

(Gross–Pitaevskii model) For p=1 and Γ(ρ)=12κρ2 with κ∈R, the energy functional takes the form6 EGP(ϕ)=12∫Ω‖∇ϕ(x)‖2dx+∫Ωϑ(x)ϕ(x)2dx+κ4∫Ωϕ(x)4dx.

This is the well-known Gross–Pitaevskii energy used in the modeling of Bose–Einstein condensates of ultracold bosonic gases [39, 43]. Here, ϑ∈L∞(Ω) is the magnetic trapping potential, ϕ∈H01(Ω) is the quantum state of the Bose–Einstein condensate, and κ characterizes the strength and the direction of particle interactions.

Example 2

(Kohn–Sham model) The (non-local) nonlinearityΓ(ρ)=ρ∫Ωρ(ϕ(y))‖x-y‖dy+2ρϵxc(ρ)

yields the Kohn–Sham energy functional7 EKS(ϕ)=12∑j=1p∫Ω‖∇ϕj(x)‖2dx+∫Ωϑion(x)ρ(ϕ(x))dx+12∫Ω∫Ωρ(ϕ(x))ρ(ϕ(y))‖x-y‖dydx+∫Ωρ(ϕ(x))ϵxc(ρ(ϕ(x)))dx,

where ϕ denotes a wave function with p components called single-particle orbitals and ρ(ϕ) is the electronic charge density. Moreover, ϑion is the ionic potential, and ϵxc(ρ) is the exchange-correlation energy per particle in a homogeneous electron gas of density ρ. This model is based on the so-called density functional theory [31], which allows a significant reduction of the degrees of freedom [17, 37, 38]. The last integral in (7) is a local density approximation to the exchange-correlation energy obtained by using semi-empirically knowledge of the model [42]. In the Kohn–Sham model, a ground state corresponds to the low-energy wave function of the considered molecule and the orthogonality condition 〚ϕ,ϕ〛H=Ip means that there is no interaction between the electrons in different orbitals.

At this point, it should be emphasized that, since the energy functional E in (3) is invariant under orthogonal transformations, i.e. E(ϕ)=E(ϕQ) for all orthogonal matrices Q∈Rp×p, the optimal solution to the minimization problem (5) is not unique. To overcome this difficulty, we will transfer this problem to the infinite-dimensional Grassmann manifold defined in Sect. 3.2.

Connection to Nonlinear Eigenvector Problems

We observe that the directional derivative of E from (3) at ϕ∈V along w∈V has the formDE(ϕ)[w]=aϕ(ϕ,w),

where8 aϕ(v,w)=∫Ωtr((∇v)T∇w)dx+2∫Ωϑv·wdx+∫Ωγ(ρ(ϕ))v·wdx

with γ(ρ)=ddρΓ(ρ). One can see that for fixed ϕ∈V, aϕ is a symmetric bilinear form on V×V. Further note that aϕ exhibits a special structure, namely9 aϕ(v,w)=∑j=1pa~ϕ(vj,wj)

with a symmetric bilinear form a~ϕ:V~×V~→R given bya~ϕ(v,w)=∫Ω(∇v)T∇wdx+2∫Ωϑvwdx+∫Ωγ(ρ(ϕ))vwdx.

Within this paper, we assume that a~ϕ is bounded and coercive on V~×V~. Obviously, the bilinear form aϕ inherits these properties such that aϕ is also bounded and coercive on V×V.

Introducing the Lagrangian L(ϕ,Λ)=E(ϕ)-12tr(ΛT(〚ϕ,ϕ〛H-Ip)) with a Lagrange multiplier Λ∈Ssym(p), the first-order necessary optimality conditions for the minimization problem (5) yield the nonlinear eigenvector problem (NLEVP) 10a aϕ∗(ϕ∗,w)-(ϕ∗Λ∗,w)H=0for allw∈V,

10b 〚ϕ∗,ϕ∗〛H-Ip=0p

with unknown ϕ∗∈V, which is referred to as the eigenvector, and Λ∗∈Ssym(p), whose eigenvalues are the lowest p eigenenergies of the system. Yet another formulation of the NLEVP () follows from the special structure of the bilinear form aϕ given in (9): seek ϕ∗=(ϕ∗,1,⋯,ϕ∗,p)∈St(p,V) and p eigenvalues λ1,⋯,λp∈R such that11 a~ϕ∗(ϕ∗,j,v)=λj(ϕ∗,j,v)L2(Ω)for allv∈V~.

For fixed ϕ∈V, we introduce the operator Aϕ:V→V∗, called the Hamiltonian, which is defined by⟨Aϕv,w⟩=aϕ(v,w)for allv,w∈V.

Then the NLEVP () can be written as 12a Aϕ∗ϕ∗-I(ϕ∗Λ∗)=0∗,

12b 〚ϕ∗,ϕ∗〛H-Ip=0p,

where 0∗∈V∗ is the zero functional. Using the left inverse J of I, we find that13 Λ∗=〚ϕ∗,ϕ∗〛HΛ∗=〚ϕ∗,ϕ∗Λ∗〛H=〚ϕ∗,JAϕ∗ϕ∗〛H.

Remark 1

Due to the symmetry of the bilinear form a~ϕ, we conclude that(ϕi,(JAϕϕ)j)L2(Ω)=a~ϕ(ϕi,ϕj)=a~ϕ(ϕj,ϕi)=(ϕj,(JAϕϕ)i)L2(Ω)

for i,j=1,…,p. This means that 〚ϕ,JAϕϕ〛H is symmetric for any ϕ∈V.

The Infinite-Dimensional Stiefel and Grassmann Manifolds

In this section, we summarize definitions and properties of the infinite-dimensional Stiefel and Grassmann manifolds and their tangent spaces, which lay the foundation of the Riemannian optimization schemes in the upcoming section.

The Stiefel Manifold

We consider the infinite-dimensional Stiefel manifold St(p,V) defined in (4). It is an embedded submanifold of the Hilbert space V and has co-dimension p(p+1)/2; see [8]. The tangent space of St(p,V) at ϕ∈St(p,V) is given by14 TϕSt(p,V)={η∈V:〚η,ϕ〛H+〚ϕ,η〛H=0p}.

The Riemannian structure of the Stiefel manifold St(p,V) strongly depends on an underlying metric. Within this paper, we equip St(p,V) with the metric given by15 g(η,ζ)=(η,ζ)H=tr〚η,ζ〛H,η,ζ∈TϕSt(p,V).

The normal space with respect to g is then defined asTϕ⊥St(p,V)={x∈V:g(x,η)=0for allη∈TϕSt(p,V)}.

It can also be represented as16 Tϕ⊥St(p,V)={ϕS∈V:S∈Ssym(p)}.

Further, any y∈V can be decomposed as y=Pϕ(y)+Pϕ⊥(y), where17 Pϕ(y)=y-ϕsym〚ϕ,y〛HandPϕ⊥(y)=ϕsym〚ϕ,y〛H

are the orthogonal projections onto the tangent and normal spaces, respectively.

The Riemannian gradient of a smooth function E:St(p,V)→R with respect to the metric g is the unique element gradE(ϕ)∈TϕSt(p,V) satisfying the conditiong(gradE(ϕ),η)=DE¯(ϕ)[η]for allη∈TϕSt(p,V),

where E¯ denotes a smooth extension of E around ϕ in V and DE¯(ϕ) is the Fréchet derivative of E¯ in V.

For the energy functional E in (3), the Riemannian gradient at ϕ∈St(p,V) with respect to the metric g can be determined by using the L2-Sobolev gradient ∇E¯(ϕ)∈V which is defined as the Riesz representation of DE¯(ϕ) in the Hilbert space V with respect to the inner product (·,·)H. Then, for all w∈V, we have⟨Aϕϕ,w⟩=aϕ(ϕ,w)=DE¯(ϕ)[w]=(∇E¯(ϕ),w)H=〈I∇E¯(ϕ),w〉

and, hence, ∇E¯(ϕ)=JAϕϕ. Furthermore, for all η∈TϕSt(p,V), we obtain(gradE(ϕ),η)H=DE¯(ϕ)[η]=(∇E¯(ϕ),η)H.

This implies that18 gradE(ϕ)=Pϕ(∇E¯(ϕ))=Pϕ(JAϕϕ)=JAϕϕ-ϕ〚ϕ,JAϕϕ〛H.

The Riemannian Hessian of E at ϕ∈St(p,V) with respect to the metric g, denoted by HessE(ϕ), is a linear mapping on the tangent space TϕSt(p,V) into itself which is defined byHessE(ϕ)[η]=∇ηgradE¯(ϕ)for allη∈TϕSt(p,V),

where ∇η denotes the covariant derivative along η with respect to the connection ∇, cf. [1, Sect. 5.3] for the finite-dimensional case.

The following theorem provides two expressions for the Riemannian Hessian of E in terms of the directional derivative of gradE(ϕ) and the L2-Sobolev Hessian ∇2E¯(ϕ) of E¯, which is a linear operator mapping v∈V onto the Riesz representation of D2E¯(ϕ)[v,·] with respect to the inner product (·,·)H.

Theorem 1

Let ϕ∈St(p,V) and η∈TϕSt(p,V). Then the Riemannian Hessian of a smooth function E:St(p,V)→R admits the expressions19 HessE(ϕ)[η]=Pϕ(DgradE(ϕ)[η])

20 =Pϕ(∇2E¯(ϕ)[η]-ηsym〚ϕ,∇E¯(ϕ)〛H),

where ∇E¯(ϕ) and ∇2E¯(ϕ) denote, respectively, the L2-Sobolev gradient and the L2-Sobolev Hessian of a smooth extension E¯ of E around ϕ in V.

Proof

Since St(p,V) is an embedded submanifold of the Hilbert space V, the expression (19) can be shown similarly to the finite-dimensional case [1, Prop. 5.3.2].

In order to prove (20), we first compute the directional derivative21 DgradE(ϕ)[η]=D(Pϕ(∇E¯(ϕ))[η]=Pϕ(∇2E¯(ϕ)[η])+DPϕ[η]∇E¯(ϕ).

Let c(t)⊂St(p,V) be a smooth curve defined on a neighborhood of t=0 such that c(0)=ϕ and ddtc(0)=η. Then for all y∈V, we haveDPϕ[η]y=limt→01t(Pc(t)(y)-Pϕ(y))=limt→01t(y-c(t)sym〚c(t),y〛H-y+c(0)sym〚c(0),y〛H)=-limt→01t(c(t)sym〚c(t)-c(0),y〛H+(c(t)-c(0))sym〚c(0),y〛H)=-ϕsym〚η,y〛H-ηsym〚ϕ,y〛H.

Inserting (21) into (19) and taking into account thatPϕ(DPϕ[η]∇E¯(ϕ))=-Pϕ(ϕsym〚η,∇E¯(ϕ)〛H+ηsym〚ϕ,∇E¯(ϕ)〛H)=-Pϕ(ηsym〚ϕ,∇E¯(ϕ)〛H),

we obtain (20). □

In order to derive a formula for the Riemannian Hessian of the energy functional E in (3), we first compute the second-order derivativeD2E¯(ϕ)[v,w]=limt→01t〈Aϕ+tv(ϕ+tv)-Aϕϕ,w〉=limt→01t(∫Ω(tr((∇(ϕ+tv))T∇w)-tr((∇ϕ)T∇w))dx+2∫Ωϑ((ϕ+tv)·w-ϕ·w)dx+∫Ω(γ(ρ(ϕ+tv))(ϕ+tv)·w-γ(ρ(ϕ))ϕ·w)dx)=∫Ωtr((∇v)T∇w)dx+2∫Ωϑv·wdx+∫Ωγ(ρ(ϕ))v·wdx+2∫Ωβ(ρ(ϕ))(ϕ·v)(ϕ·w)dx=⟨Aϕv+Bϕv,w⟩,

where β(ρ)=ddργ(ρ) and the operator Bϕ:V→V∗ has the form22 ⟨Bϕv,w⟩=2∫Ωβ(ρ(ϕ))(ϕ·v)(ϕ·w)dx.

Hence, the L2-Sobolev Hessian of E¯ is given by ∇2E¯(ϕ)[v]=JAϕv+JBϕv for all v∈V. By the definition of the orthogonal projection onto TϕSt(p,V) in (17), we conclude that for η∈TϕSt(p,V), the Riemannian Hessian of E is given by23 HessE(ϕ)[η]=Pϕ(JAϕη+JBϕη-η〚ϕ,JAϕϕ〛H)=JAϕη+JBϕη-η〚ϕ,JAϕϕ〛H-ϕsym〚ϕ,JAϕη〛H-ϕsym〚ϕ,JBϕη〛H+ϕsym(〚ϕ,η〛H〚ϕ,JAϕϕ)〛H).

Within optimization methods, we need to transfer data from the tangent space to the manifold to keep the iterations on the search space. For this purpose, we can use retractions defined as follows. Let TSt(p,V) be the tangent bundle to St(p,V). A smooth mapping R:TSt(p,V)→St(p,V) is called a retraction if for all ϕ∈St(p,V), the restriction of R to TϕSt(p,V), denoted by Rϕ, satisfies the following properties: Rϕ(0ϕ)=ϕ, where 0ϕ denotes the origin of TϕSt(p,V), and

ddtRϕ(tη)|t=0=η for all η∈TϕSt(p,V).

Retractions provide first-order approximations to the exponential mapping on a Riemannian manifold and are often much easier to compute. A retraction R on St(p,V) is of second-order, if it satisfies d2dt2Rϕ(tη)|t=0∈Tϕ⊥St(p,V) for all (ϕ,η)∈TSt(p,V).

In [8], several retractions on the Stiefel manifold St(p,V) have been introduced. They can be considered as an extension of the corresponding concepts on the matrix Stiefel manifold (see, e.g., [2, 46]) to the infinite-dimensional case.

For v∈V with linearly independent components, we consider the qR decomposition v=qR, where q∈St(p,V) and R∈Rp×p is upper triangular. Such a decomposition exists and is unique if we additionally require that R has positive diagonal elements. Then the qR decomposition based retraction is defined as RqR(ϕ,η)=qf(ϕ+η), where qf(ϕ+η) denotes the factor from St(p,V) in the qR decomposition of ϕ+η. Such a factor can be computed, e.g., by the modified Gram-Schmidt orthonormalization procedure on V presented in [8].

An alternative retraction can be defined by using the polar decomposition v=uS, where u∈St(p,V) and S∈Rp×p is symmetric and positive definite. Assuming that the components of v are linearly independent, S=〚v,v〛H1/2 and u=v〚v,v〛H-1/2 are uniquely defined. This leads to the polar decomposition based retractionRpol(ϕ,η)=(ϕ+η)〚ϕ+η,ϕ+η〛H-1/2,

which is of second-order. Indeed, computing the second-order derivative of Rϕpol(tη) at t=0 and exploiting (16), we obtain thatd2dt2Rϕpol(tη)|t=0=-ϕ〚η,η〛H∈Tϕ⊥St(p,V).

Note that second-order retractions are advantageous for second-order Riemannian optimization methods; see, e.g. [1, Sect. 6.3].

The Grassmann Manifold

Let O(p) be the orthogonal group of Rp×p. Following [47], we define the infinite-dimensional Grassmann manifold as the quotientGr(p,V)=St(p,V)/O(p)

of the Stiefel manifold St(p,V) with respect to the equivalence relationϕ∼ϕ^⟺ϕ^=ϕQforsomeQ∈O(p).

The Grassmann manifold Gr(p,V) can be interpreted as the set of the equivalence classes given by

for ϕ∈St(p,V). Similarly to the Grassmann matrix manifold [1, Prop. 3.4.6], one can show that Gr(p,V) admits a unique structure of quotient manifold. A canonical projection from the Stiefel manifold into the Grassmann manifold is defined by

and is a smooth submersion. This means that Dπ(ϕ) is surjective, and, hence, the equivalence class is an embedded submanifold of St(p,V); see [1, Prop. 3.4.4.].

In the following, we examine a useful connection of the Stiefel manifold and the Grassmann manifold. More precisely, we show that there is a one-to-one relation between the tangent space of the Grassmann manifold and the so-called horizontal space, a subspace of the tangent space of the Stiefel manifold. The tangent space TϕSt(p,V) at ϕ∈St(p,V) defined in (14) can be splitted with respect to the projection π and the metric g as TϕSt(p,V)=Vϕ⊕Hϕ, where 24

is the vertical space at ϕ and Hϕ=Vϕ⊥={x∈TϕSt(p,V):g(x,v)=0for allv∈Vϕ}={x∈TϕSt(p,V):〚ϕ,x〛H=0p}

is the horizontal space at ϕ; see [47, Lem. 2]. The orthogonal projection of a tangent vector η∈TϕSt(p,V) onto Hϕ is given by 25 Pϕh(η)=η-ϕ〚ϕ,η〛H.

One can see that, moving on a curve in the Stiefel manifold St(p,V) with direction in the vertical space Vϕ, we stay in the equivalence class . The tangent space of the Grassmann manifold Gr(p,V) can then be identified with the horizontal space Hϕ in the sense that for any , there exists a unique ψϕh∈Hϕ such that Dπ(ϕ)[ψϕh]=ψ. The unique element ψϕh is called the horizontal lift of ψ at ϕ. This relation allows us to introduce a metric on the Grassmann manifold Gr(p,V), namely

where ψϕh,ζϕh∈Hϕ are the horizontal lifts of ψ and ζ at ϕ, respectively. Due to ψϕQh=ψϕhQ for all Q∈O(p), one can show that this metric does not depend on the choice of the representative  ϕ of the equivalence class  .

The connection of and Hϕ makes it possible to introduce optimization methods on the Grassmann manifold, while still working on the tangent space of the corresponding Stiefel manifold. Using the canonical projection π, the minimization problem (5) on the Stiefel manifold St(p,V) can be written as the minimization problem 26

on the Grassmann manifold Gr(p,V), where the cost functional F:Gr(p,V)→R is induced by E as E(ϕ)=F(π(ϕ)) and . Note that this definition is justified by the fact that E(ϕ)=E(ϕQ) for all  Q∈O(p). The horizontal lift of the Riemannian gradient with respect to the metric gGr is given by 27

To obtain the horizontal lift of the Riemannian Hessian , we proceed as before but replace the projection  Pϕ by the orthogonal projection  Pϕh onto the horizontal space; see Eq. (25). This leads to 28

Retractions on the Grassmann manifold are inherited from that on the Stiefel manifold applied to the horizontal lift; see [1, Prop. 4.1.3]. For all and , we have

Note that these retractions are independent of the chosen point ϕ, providing the same equivalence class on Gr(p,V).

Similar to the matrix case [1, 25], we can also derive an explicit expression for the Grassmann exponential Exp:TGr(p,V)→Gr(p,V), which maps to the end point of the unique geodesic starting at and going in the direction  ψ. Let ψϕh=uΣWT be a singular value decomposition of the horizontal lift ψϕh of  ψ, where u∈St(p,V), W∈O(p), and Σ∈Rp×p is diagonal with nonnegative diagonal elements. Then the Grassmann exponential is given by

Using ψϕh∈Hϕ, one can verify that ϕWcosΣ+usinΣ∈St(p,V). Therefore, it can be considered as a representative of the resulting equivalence class.

Remark 2

For p=1, the Stiefel manifold coincides with the Grassmann manifold and equals the unit sphere S={ϕ∈V~:‖ϕ‖L2(Ω)=1}.

Its tangent space is given by TϕS={η∈V~:(η,ϕ)L2(Ω)=0} and the orthogonal projection onto this space takes the form Pϕ(y)=y-(ϕ,y)L2(Ω)ϕ for y∈V~. Furthermore, for (ϕ,η)∈TS, the second-order retraction and the exponential mapping on S are given by R(ϕ,η)=ϕ+η‖ϕ+η‖L2(Ω),Exp(ϕ,η)=cos(‖η‖L2(Ω))ϕ+sin(‖η‖L2(Ω))η‖η‖L2(Ω),

respectively.

Riemannian Newton Methods

In this section, we present Riemannian Newton methods on the Stiefel manifold as well as on the Grassmann manifold and discuss the inexact version of the latter.

Within the Riemannian Newton method on the Stiefel manifold St(p,V), for given iterate ϕk∈St(p,V), we first compute the Newton search direction ηk∈TϕkSt(p,V) by solving the Newton equation 29 HessE(ϕk)[ηk]=-gradE(ϕk).

The iterate is then updated by  ϕk+1=R(ϕk,ηk) for any retraction on St(p,V). It should, however, be noted that due to the non-uniqueness of the minimizer of (5) caused by the invariance of E under orthogonal transformations, we cannot expect that HessE(ϕk) is invertible on TϕkSt(p,V). This is, indeed, vindicated by the following theorem.

Theorem 2

Let ϕ∈St(p,V) and let Vϕ be the vertical space given in (24). Then the Riemannian Hessian of the energy functional  E from (3) in  ϕ is non-invertible on  Vϕ, i.e., (ξ,HessE(ϕ)[η])H=0 for all η,ξ∈Vϕ.

Proof

Let η,ξ∈Vϕ be arbitrary. Then there exist the matrices Θη,Θξ∈Sskew(p) such that η=ϕΘη and ξ=ϕΘξ. Using the definition of HessE(ϕ) in (23) and the symmetry of the matrix  〚ϕ,JAϕϕ〛H shown in Remark 1, we have HessE(ϕ)[η]=JAϕϕΘη+JBϕϕΘη-ϕΘη〚ϕ,JAϕϕ〛H-ϕsym(〚ϕ,JAϕϕ〛HΘη)-ϕsym(〚ϕ,JBϕϕ〛HΘη)+ϕsym(Θη〚ϕ,JAϕϕ〛H)=JAϕϕΘη-ϕ〚ϕ,JAϕϕ〛HΘη+JBϕϕΘη-ϕsym(〚ϕ,JBϕϕ〛HΘη)

and, hence, (ξ,HessE(ϕ)[η])H=trΘξT〚ϕ,JBϕϕ〛HΘη-ΘξTsym(〚ϕ,JBϕϕ〛HΘη).

We now show that 〚ϕ,JBϕϕ〛HΘη=0. With the skew-symmetric matrix Θη=[θij]i,j=1p, we first observe that ϕ·(ϕΘη)=∑j=1pϕj∑i=1pϕiθij=∑i=1pϕi∑j=1pϕjθij=-∑i=1pϕi∑j=1pϕjθji=-ϕ·(ϕΘη).

This implies ϕ·(ϕΘη)=0 and, hence, 〚ϕ,JBϕϕ〛HΘη=0. As a result, we conclude that (ξ,HessE(ϕ)[η])H=0 for all η,ξ∈Vϕ. □

It follows from Theorem 2 that if Eq. (29) is solvable, its solution is not unique. To overcome this difficulty, we pass on to the Grassmann manifold Gr(p,V). Given , the Newton direction is computed by solving the Newton equation

By applying the horizontal lift expressions (27) and (28), this equation leads to 30 Pϕkh(DgradE(ϕk)[(ψk)ϕkh])=-gradE(ϕk)

with unknown (ψk)ϕkh∈Hϕk being the horizontal lift of  ψk at  ϕk. Note that this equation is well-defined, since  -gradE(ϕk) is an element of the horizontal space Hϕk.

For solving the Newton equation (30), we can employ any matrix-free iterative linear solver which does not require the storage of the coefficient matrix explicitly but accesses it by computing the matrix–vector product or—as in our case—by evaluating the linear operator given in (28). The resulting Riemannian Newton method on the Grassmann manifold is presented in Algorithm 1. Algorithm 1 Riemannian Newton method on the Grassmann manifold

The Newton equation (30) can also be formulated as a saddle point problem. To this end, we introduce the bilinear form a^ϕ(ψ,w)=⟨Aϕψ,w⟩+⟨Bϕψ,w⟩-(ψ〚ϕ,JAϕϕ〛H,w)H.

Then, the equivalent problem to (30) reads: find (ψk)ϕkh∈V and a Lagrange multiplier Mk∈Rp×p such that 31a a^ϕk((ψk)ϕkh,w)+tr(MkT〚ϕk,w〛H)=-aϕk(ϕk,w)for allw∈V,

31b 〚ϕk,(ψk)ϕkh〛H=0p.

The constraint (31b) implies that (ψk)ϕkh∈Hϕk. Further note that any function w∈Hϕk satisfies  〚ϕk,w〛H=0p. Hence, for all w∈Hϕk, Eq. (31a) reads a^ϕ((ψk)ϕkh,w)=-aϕk(ϕk,w)=-(gradE(ϕk),w)H,

which is equivalent to the Newton equation (30).

One important property guaranteeing an isolated local minimum of the energy is that the Hessian is positive at a stationary point. For a global minimizer of (5), denoted by  ϕ∗, we consider the following linear eigenvalue problem: seek ϕ∈V~ and λ∈R such that 32 a~ϕ∗(ϕ,v)=λ(ϕ,v)L2(Ω)for allv∈V~.

Then, due to (11), we know that the components of  ϕ∗=(ϕ∗,1,⋯,ϕ∗,p) satisfy Eq. (32) together with the smallest p eigenvalues denoted by 0<λ1≤⋯≤λp. In the following, we will assume that these eigenfunctions can be extended to a basis of  V~.

Assumption 3

(Basis and spectral gap) The eigenfunctions  ϕ∗,1,ϕ∗,2,…∈V~ of the eigenvalue problem (32) form an L2-orthonormal basis of V~. The corresponding eigenvalues  λ1≤λ2≤⋯ are ordered by size with a spectral gap  λp<λp+1.

Theorem 4

(Positive Hessian) Let  ϕ∗ be a global minimal solution of (5) and let the corresponding eigenvalue problem (32) satisfy Assumption 3. Furthermore, assume that the operator Bϕ∗ fulfills  (JBϕ∗ψϕ∗h,ψϕ∗h)H≥0 for all ψϕ∗h∈Hϕ∗. Then the Riemannian Hessian of  F at is positive, i.e.,

for all nonzero  .

Proof

We extend the proof of [54, Th. 5.1], which considers the finite-dimensional case for the simplified Kohn–Sham problem, to the infinite-dimensional case in a more general setting. We know from (28) that the horizontal lift of the Riemannian Hessian of F at ϕ∗ takes the form

For the first term, we make the following considerations. Due to the orthogonality, each component of ψϕ∗h satisfies (ϕ∗,l,(ψϕ∗h)j)L2(Ω)=0 for l,j=1,⋯,p. Hence, the jth component of ψϕ∗h takes the form ψj=(ψϕ∗h)j=∑l>pαljϕ∗,l for some coefficients  αlj∈R and  ϕ∗,l denoting the basis from Assumption 3. As a consequence, we get (〚ϕ∗,JAϕ∗ψϕ∗h〛H)ij=(ϕ∗,i,(JAϕ∗ψϕ∗h)j)L2(Ω)=a~ϕ∗(ϕ∗,i,ψj)=∑l>pαlja~ϕ∗(ϕ∗,i,ϕ∗,l)=∑l>pαljλi(ϕ∗,i,ϕ∗,l)L2(Ω)=0

for all i,j=1,…,p and, hence, T1(ψϕ∗h)=Pϕ∗h(JAϕ∗ψϕ∗h)=JAϕ∗ψϕ∗h-ϕ∗〚ϕ∗,JAϕ∗ψϕ∗h〛H=JAϕ∗ψϕ∗h.

For the second term, we get with the assumption on  Bϕ∗ that g(T2(ψϕ∗h),ψϕ∗h)=(T2(ψϕ∗h),ψϕ∗h)H=(JBϕ∗ψϕ∗h-ϕ∗〚ϕ∗,JBϕ∗ψϕ∗h〛H,ψϕ∗h)H=(JBϕ∗ψϕ∗h,ψϕ∗h)H-〚ϕ∗,JBϕ∗ψϕ∗h〛HT(ϕ∗,ψϕ∗h)H=(JBϕ∗ψϕ∗h,ψϕ∗h)H≥0.

Finally, by using (13), the third term takes the form T3(ψϕ∗h)=-Pϕ∗h(ψϕ∗h〚ϕ∗,JAϕ∗ϕ∗〛H)=-ψϕ∗h〚ϕ∗,JAϕ∗ϕ∗〛H+ϕ∗〚ϕ∗,ψϕ∗h〛H〚ϕ∗,JAϕ∗ϕ∗〛H=-ψϕ∗hΛ∗.

Let the columns of U∈O(p) form a basis of eigenvectors corresponding to the eigenvalues λ1,…,λp of Λ∗ and let ψϕ∗hU=(ψ~1,…,ψp~). Due to the assumed spectral gap, this yields all together

which completes the proof. □

Remark 3

(Connection to the Lagrange–Newton method) The optimal solution of the constrained minimization problem (5) can also be determined by the Lagrange–Newton method. Based on the first-order optimality conditions () with a symmetric Lagrange multiplier, we aim to solve the nonlinear system of equations f(ϕ,Λ)=JAϕϕ-ϕΛ〚ϕ,ϕ〛H-IpΛ-ΛT=00p0p.

Computing the Jacobian of f, the Lagrange-Newton iteration is given as follows: for given ϕk∈V and Λk∈Rp×p, solve the equations 33a JAϕkηk+JBϕkηk-ηkΛk-ϕkΞk=-(JAϕkϕk-ϕkΛk),

33b 〚ϕk,ηk〛H+〚ηk,ϕk〛H=-(〚ϕk,ϕk〛H-Ip),

33c Ξk-ΞkT=-(Λk-ΛkT)

for ηk∈V, Ξk∈Rp×p and update ϕk+1=ϕk+ηk, Λk+1=Λk+Ξk. Note that ϕk+1 does not necessarily belong to St(p,V). Assuming ϕk∈St(p,V) and Λk∈Ssym(p), however, Eqs. (33b) and (33c) imply that ηk∈TϕkSt(p,V) and Ξk∈Ssym(p), respectively. Resolving Eq. (33a) for symmetric Ξk, we find that Ξk=sym(〚ϕk,JAϕkηk〛H+〚ϕk,JBϕkηk〛H-〚ϕk,ηk〛HΛk)+〚ϕk,JAϕkϕk〛H-Λk.

Inserting this matrix into (33a) yields the Newton equation (29). This shows that the Lagrange–Newton method with the modified update ϕk+1=R(ϕk,ηk),Λk+1=〚ϕk+1,JAϕk+1ϕk+1〛H

is equivalent to the Riemannian Newton method on the Stiefel manifold.

Examples and Numerical Experiments

This section is devoted to the numerical investigation of the Riemannian Newton methods. To this end, we consider the Gross–Pitaevskii eigenvalue problem from Example 1 and the Kohn–Sham model from Example 2.

Gross–Pitaevskii Eigenvalue Problem

The minimization of the Gross–Pitaevskii energy functional EGP in (6) leads to the following nonlinear eigenvector problem: find ϕ∈V~=H01(Ω) with ‖ϕ‖L2(Ω)=1 and λ∈R such that 34 -Δϕ+2ϑϕ+κ|ϕ|2ϕ=λϕ

for some space-dependent external potential  ϑ≥0 and an interaction constant κ>0. The latter means that the particle interactions are repulsive, i.e., we consider the so-called defocussing regime. In this case, we get the operators 35a ⟨Aϕv,w⟩=∫Ω(∇v)T∇wdx+2∫Ωϑvwdx+κ∫Ωϕ2vwdx,

35b ⟨Bϕv,w⟩=2κ∫Ωϕ2vwdx

for v,w∈V~. One can see that the bilinear form defined through (35a) corresponds to the Laplacian with the L2-shift  2ϑ+κϕ2. Assuming this shift to be constant and Ω=(0,1)d as the spatial domain, Assumption 3 is satisfied; see [50, Ch. 12]. Moreover, the nonlinear operator from (35b) fulfills (JBϕψ,ψ)L2(Ω)=〈Bϕψ,ψ〉=2κ∫Ωϕ2ψ2dx≥0,

such that Theorem 4 is applicable.

For the spatial discretization of the Gross–Pitaevskii problem (34), we use a biquadratic finite element method on a Cartesian mesh of width h; see [16] for the corresponding error analysis. The resulting discrete eigenvalue problem reads Aφ+2Mϑφ+κMφ2φ=λMφ,φTMφ=1

with φ∈Rn, where n denotes the number of degrees of freedom. Here, A is the stiffness matrix, M is the mass matrix, and Mϑ and Mφ2 are the weighted mass matrices, respectively, where  φ2 should be understood as the elementwise product. Then the discrete version of  〚ϕ,JAϕϕ〛H=(ϕ,JAϕϕ)L2(Ω) equals λφ=φT(A+2Mϑ+κMφ2)φ,

and the Newton equation takes the form (I-MφφT)((A+2Mϑ+3κMφ2)ψ-λφMψ)=-(I-MφφT)(A+2Mϑ+κMφ2)φ

with unknown ψ∈{ξ∈Rn:ξTMφ=0}=im(I-φφTM).

We demonstrate the performance of the resulting Riemannian Newton method in comparison with the SCF iteration combined with the optimal damping algorithm (ODA) proposed in [23] and the energy-adaptive Riemannian gradient descent method (RGD) of [29] with a non-monotone step size control as outlined in [8]. The numerical experiments are performed on a sufficiently large bounded domain Ω=(-L,L)2, L=8, for two types of trapping potentials. In Sect. 5.1.1, we consider a simple harmonic trap, whereas in Sect. 5.1.2, we add an additional disorder potential. The interaction parameter  κ as well as the spatial resolution h will be specified separately for each case.

Ground State in a Harmonic Trap

For the first numerical experiment, we consider the harmonic trapping potential 36 ϑharm(x)=12‖x‖2

and the interaction parameters κ=10,100,1000. The resulting ground states computed on a Cartesian mesh of width h/(2L)=2-10 are depicted in Fig. 1.Fig. 1 Ground state in the harmonic trap (potential in gray, properly rescaled) for κ=10,100,1000 (from left to right)

To generate a joint and sufficiently accurate initial value for the three solvers of the discretized nonlinear eigenvector problem, we run the energy-adaptive RGD method starting from the biquadratic finite element interpolation of the constant 1 (respecting the homogeneous Dirichlet boundary condition). We stopped this iteration once the residual fell below the 10-2 tolerance and used the approximated ground state as the initial state to compare the asymptotic behavior of the three different solvers. The corresponding convergence histories are presented in Fig. 2 showing the evolution of the residuals during the iteration processes. It can be observed that the Riemannian Newton method (with sparse direct solution of the Newton equation using the Sherman–Morrison formula [48]) reaches the tolerance of 10-8 in only three steps. While the performances of the SCF iteration and the energy-adaptive RGD method abate with increasing κ, the Riemannian Newton scheme appears to be extremely robust. We would like to emphasize that, although one Newton step is slightly more expensive than one step of any other competing method, the overall costs are much smaller, especially for increasing  κ.Fig. 2 Convergence history of the residuals for the ground state in the harmonic trap for κ=10,100,1000 (from left to right)

The convergence behavior of the Riemannian Newton scheme is also robust to the underlying mesh size h and, hence, independent of the dimension of the discretization space. This is demonstrated in Fig. 3 with a fixed choice of κ=1000. We consider a sequence of meshes with h/(2L)=2-1,…,2-10 and use the same procedure as above to generate initial guesses with residuals of order 10-2. The left graph shows the number of (outer) iterations of the Riemannian Newton method to fall below the tolerance of 10-10 for each of these mesh sizes. An increase in the number of iterations with smaller mesh size is not observed. In our experience, this mesh independence of the Riemannian Newton optimization scheme is representative for many other choices of potentials and interaction parameters. Note, however, that this does not mean that the costs of a Newton step are independent of h. As for every Laplace-type problem, methods such as multigrid need to be implemented in order to obtain a mesh-independence also for the inner iteration. This holds for all competing methods in the same way. For completeness, Fig. 3 also shows the corresponding errors in the minimum energy approximation as a function of the mesh size, demonstrating the optimal fourth-order convergence rate of the biquadratic finite element implementation [30].Fig. 3 Computing the ground state in the harmonic trap for κ=1000: iteration count of the Riemannian Newton method to fall below the tolerance 10-10 (left) and the error in the minimal energy (right) versus mesh size h. The dashed line indicates order h4

Localized Ground State in a Disorder Potential

The second experiment considers the computationally more difficult case where the external potential is the sum of the harmonic potential ϑharm defined in (36) and a potential ϑrand reflecting a high degree of disorder. The disorder part ϑrand is chosen as a piecewise constant function on the Cartesian mesh of width 2Lε, ε=2-6, taking values 0 or ε-2 as depicted in Fig. 4. For a potential in such a scaling regime, the low-energy eigenstates essentially localize in terms of an exponential decay of their moduli relative to the small parameter ε. For the linear case, i.e., for κ=0, this has been analyzed in [4]. For growing κ, the ground state consists of a growing number of localized peaks; see Fig. 4. Further details on the phenomenon of localization in the Gross–Pitaevskii equation and the onset of delocalization can be found in [6, 7]. As in the previous experiments, we use biquadratic finite elements on a Cartesian mesh of width h/(2L)=2-10. To illustrate the localization behavior that occurs with the current parameter scaling for κ≲1, we consider the interaction parameters κ=0.1,1,10. The ground states for κ=1,10 are shown in Fig. 4. The ground state for κ=0.1 is hardly distinguishable from the one for κ=1 and, therefore, it is not shown in a separate figure.Fig. 4 Piecewise constant disorder potential ϑrand (left, black elements refer to the value ε-2, white elements refer to the value 0, ε=2-6) and the corresponding ground states for κ=1 (middle) and κ=10 (right)

Figure 5 displays the convergence history of the residuals for κ=0.1,1,10. We employed the same strategy as in Sect. 5.1.1 to generate suitable initial guesses with the residuals of order 10-2 used for all methods. The results clearly indicate that the ground state computations with the disorder potential are already challenging for smaller values of κ. Particularly, the energy-adaptive RGD method needs much larger iteration counts, which according to [27], may be related to smaller spectral gaps between the first and second eigenvalue. The Riemannian Newton method, on the other hand, still performs well and reaches the prescribed tolerance 10-8 for the residual in only a few steps in all three examples. For comparison, the SCF iteration converges very fast in the almost linear case but suffers from larger values of κ as the energy-adaptive RGD method. Here, again, the higher costs per Newton step are compensated by far by the very small number of needed iteration steps.Fig. 5 Convergence history of the residuals for the ground state in a disorder potential for κ=0.1,1,10 (from left to right)

Kohn–Sham Model

For the Kohn–Sham energy functional EKS introduced in (7), we have ⟨Aϕv,w⟩=∫Ωtr((∇v)T∇w)dx+2∫Ωϑionv·wdx+2∫Ω(∫Ωρ(ϕ(y))‖x-y‖dy)v·wdx+2∫Ωμxc(ρ(ϕ))v·wdx

with μxc(ρ)=ddρ(ρϵxc(ρ)). Moreover, the operator Bϕ has the form ⟨Bϕv,w⟩=4∫Ω(∫Ωϕ·v‖x-y‖dy)ϕ·wdx+4∫Ωζxc(ρ(ϕ))(ϕ·v)(ϕ·w)dx,

where ζxc(ρ)=ddρμxc(ρ). The exchange-correlation function ϵxc(ρ) can additively be decomposed as ϵxc(ρ)=ϵx(ρ)+ϵc(ρ), where the exchange component ϵx(ρ) has the particular analytical expression ϵx(ρ)=-34(3πρ)1/3 and the correlation component ϵc(ρ) is usually unknown, but can be fitted by using quantum Monte-Carlo data [41]. For the numerical experiments, we use the MATLAB toolbox KSSOLV [35, 51], in which the correlation component is implemented as ϵc(ρ)=a1+a2r(ρ)+(a3+a4r(ρ))ln(r(ρ)),ifr(ρ)<1,[0.3em](1+b2r(ρ)+b3r(ρ))/b1,ifr(ρ)≥1,

where r(ρ)=(4π3ρ)-1/3 is the Wigner-Seitz radius, and aj,bj∈R are fitted constants; see [42, App. C].

For the spatial discretization, we employ the planewave discretization method as implemented in KSSOLV. With n denoting the number of degrees of freedom, the matrix Φ∈Cn×p contains the coefficients of the approximation of the wave function ϕ. Then the discretized Kohn–Sham energy functional is given by E(Φ)=12tr(Φ∗(L+2Dion)Φ)+12ρh(Φ)TL+ρh(Φ)+ρh(Φ)Tϵxc(ρh(Φ)),

where Φ∗ denotes the complex conjugate transpose of Φ, L∈Cn×n is the discrete Laplace matrix, L+∈Cn×n is its pseudoinverse, Dion∈Rn×n is the discretized ionic potential, and ρh(Φ)=diag(ΦΦ∗)∈Rn is the discretized electronic charge density. Note that the matrix L is Hermitian and Dion is diagonal. In this setting, the minimization problem minΦ∈St(p,n)E(Φ)

on the (compact) Stiefel manifold St(p,n)={Φ∈Cn×p:Φ∗Φ=Ip} leads to the finite-dimensional nonlinear eigenvector problem A(Φ)Φ-ΦΛ=0,Φ∗Φ-Ip=0,

where the discrete Kohn–Sham Hamiltonian is given by A(Φ)=L+2Dion+2Diag(L+ρh(Φ)+μxc(ρh(Φ))).

Further, the Riemannian gradient of E(Φ) becomes 37 gradE(Φ)=(I-ΦΦ∗)A(Φ)Φ=A(Φ)Φ-Φ(Φ∗A(Φ)Φ).

In the Riemannian Newton method on the Stiefel manifold St(p,n), we need to solve the equation 38 PΦ(A(Φ)Ψ+B(Φ,Ψ)-ΨΦ∗A(Φ)Φ)=-(I-ΦΦ∗)A(Φ)Φ

for Ψ belonging to the tangent space TΦSt(p,n). Therein, PΦ(Y)=Y-12Φ(Φ∗Y+Y∗Φ)

is the orthogonal projector onto TΦSt(p,n) and B(Φ,Ψ)=2Diag((L++Diag(ζxc(ρh(Φ))))diag(ΦΨ∗+ΨΦ∗))Φ

is the discretization of Bϕ. On the Grassmann manifold Gr(p,n)=St(p,n)/U(p) with the unitary group U(p), the Newton equation takes the form 39 (I-ΦΦ∗)(A(Φ)ΨΦh+B(Φ,ΨΦh)-ΨΦhΦ∗A(Φ)Φ)=-(I-ΦΦ∗)A(Φ)Φ

for ΨΦh∈HΦ={Ψ∈TΦSt(p,n):Φ∗Ψ=0}.

In our experiments, we compare the calculation of the ground state by using the SCF iteration with the Anderson charge mixing scheme [51], the energy-adaptive RDG with non-monotone step size control [8], and the Riemannian Newton methods on the Stiefel manifold (RNS) and on the Grassmann manifold (RNG). In all these methods, we choose the same initial guess for the wave function by performing one SCF step with a randomly generated starting point and stop the iterations once the Frobenius norm of the residual R(Φk)=A(Φk)Φk-ΦkΛk with Λk=Φk∗A(Φk)Φk is smaller than the tolerance 10-8. Note that due to (37), ‖R(Φk)‖F=‖gradE(Φk)‖F, i.e., the norms of the residuals provide the information on the size of the Riemannian gradients. In both Newton methods and the energy-adaptive RGD method, we use the qR decomposition based retractions. The reference minimal energy Emin is computed by the RNG method with the tolerance 10-10.

All algorithms are performed in an inexact manner, i.e., the occurring linear systems are only solved up to a certain tolerance. In RNG, for instance, we follow Algorithm 1 using the MATLAB built-in function minres as a linear system solver with the adaptive tolerance min(1/k,10-3‖gradE(Φk-1)‖F) and the maximal number of inner iterations ℓmax=15. The remaining parameters are chosen as η=10-8, δ=0.5, and σ=10-4. In RNS, we proceed similarly, with the only difference that instead of (39) we solve the Newton equation of the form (38). For solving the linear eigenvalue problems in SCF, we employ the KSSOLV built-in LOBPCG algorithm for the pentacene model in Sect. 5.2.1 and the MATLAB built-in function eigs for the graphene model in Sect. 5.2.2. Switching to another eigenvalue solver is necessary due to the ill-conditioning in LOBPCG for the latter example. In both cases, the tolerance for the inner iterations is set to be min(10-3,10-3‖gradE(Φk-1)‖F). For the linear system solvers, we use the kinetic energy preconditioner, which provides, especially for the energy-adaptive RGD method, better numerical results than the KSSOLV built-in Teter–Payne–Allan preconditioner [49].

Finally, we would like to mention that similarly to the Gross–Pitaevskii example in Sect. 5.1, the Riemannian Newton schemes for the Kohn–Sham model are again robust in terms of the dimension of the discretization space used in KSSOLV. This means that the number of Newton iterations needed to fall below a certain tolerance is not effected by finer discretizations (as long as the number of inner iterations is sufficiently large).

Pentacene Molecule

In the first numerical experiment, we calculate the ground state for the pentacene molecule C22H14 with p=51 electron orbitals. A spatial planewave discretization on a 80×55×160 sampling grid gives the discrete model of dimension n=44791. In Fig. 6, we present the convergence history of the residuals and the energy reduction during the iterations. One can see that the RNS and RNG methods have very similar behavior and converge within 8 and 10 iterations, respectively. In comparison, the SCF method and the energy-adaptive RGD method require 18 and 54 iterations to converge, respectively. In terms of computing time, all methods perform quite similarly in this experiment. This can also be seen in Table 1, which shows the values of the energy functional, the reached residuals, the number of (outer) iterations, the total number of Hamiltonian evaluations, and the CPU time.Fig. 6 Convergence history for the pentacene molecule: residuals (left) and energy reduction (right)

Fig. 7 Convergence history for the graphene lattice: residuals (left) and energy reduction (right)

Table 1 Numerical results for the pentacene and the graphene models

Method	Energy	Residual	# iter	# Ham. eval.	CPU time [s]	
	Pentacene molecule	n=44791, p=51	
RNS	-1.3189e+2	3.2731e-9	10	208	1880.64	
RNG	-1.3189e+2	6.4308e-9	8	172	1572.17	
SCF	-1.3189e+2	7.4357e-9	18	276	1709.78	
eaRGD	-1.3189e+2	6.9843e-9	54	369	1739.08	
	Graphene lattice	n=12279, p=67	
RNS	-1.7360e+2	5.5035e-9	8	657	568.69	
RNG	-1.7360e+2	1.0569e-9	10	693	741.55	
SCF	-1.7312e+2	4.8657e-4	100	69616	3422.67	
eaRGD	-1.7360e+2	9.4645e-9	64	914	991.22	

Graphene Lattice

As the second model, we consider a graphene lattice consisting of carbon atoms arranged in 9 hexagons with p=67 electron orbitals. We use a 32×55×160 sampling grid for the wave function and get a discrete model of dimension n=12279. Figure 7 presents the evolution of the residuals and errors in the energy. We observe again that both Newton methods converge very fast compared to the energy-adaptive RGD method which needs 64 iterations to achieve the tolerance 10-8 for the residual. In contrary, the SCF iteration has difficulties to converge. Also other mixing strategies implemented in KSSOLV do not improve the convergence property of SCF for the graphene model. This behaviour may be explained by a missing spectral gap between the excited and non-excited states. A detailed comparison, including the overall CPU time is part of Table 1. In this experiment (with the particular implementation and used hardware), one can say that the computational complexity of the methods follows the rule of thumb 1 step Newton≈2 steps SCF≈4 steps eaRGD.

Overall, this example clearly shows the supremacy of the Newton approach for more challenging examples.

Conclusion

In this paper, we have derived Riemannian Newton methods on the infinite-dimensional Stiefel and Grassmann manifolds for Kohn–Sham type energy minimization problems. Starting from an energy functional, we present a unified approach for applications in computational physics (e.g., the Gross–Pitaevskii eigenvalue problem) and computational chemistry (e.g., the Kohn–Sham model). The remarkable gain in computational efficiency of the Riemannian Newton methods compared to the so far more popular methods such as SCF and gradient descent methods is demonstrated by a series of numerical experiments.

Funding

Open Access funding enabled and organized by Projekt DEAL. The work of D. Peterseim is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 865751 – RandomMultiScales).

Data Availability

Data is available from the corresponding author on reasonable request.

Declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Absil PA Mahony R Sepulchre R Optimization Algorithms on Matrix Manifolds 2008 Princeton, NJ Princeton University Press
Absil, P.A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ (2008)
2. Absil PA Malick J Projection-like retractions on matrix manifolds SIAM J. Optim. 2012 22 1 135 158 10.1137/100802529
Absil, P.A., Malick, J.: Projection-like retractions on matrix manifolds. SIAM J. Optim. 22(1), 135–158 (2012). 10.1137/100802529
3. Alouges F Audouze C Preconditioned gradient flows for nonlinear eigenvalue problems and application to the Hartree-Fock functional Numer. Methods Partial Differ. Equ. 2009 25 2 380 400 10.1002/num.20347
Alouges, F., Audouze, C.: Preconditioned gradient flows for nonlinear eigenvalue problems and application to the Hartree-Fock functional. Numer. Methods Partial Differ. Equ. 25(2), 380–400 (2009). 10.1002/num.20347
4. Altmann R Henning P Peterseim D Quantitative Anderson localization of Schrödinger eigenstates under disorder potentials Math. Models Methods Appl. Sci. 2020 30 5 917 955 10.1142/S0218202520500190
Altmann, R., Henning, P., Peterseim, D.: Quantitative Anderson localization of Schrödinger eigenstates under disorder potentials. Math. Models Methods Appl. Sci. 30(5), 917–955 (2020). 10.1142/S0218202520500190
5. Altmann R Henning P Peterseim D The J-method for the Gross–Pitaevskii eigenvalue problem Numer. Math. 2021 148 575 610 10.1007/s00211-021-01216-5
Altmann, R., Henning, P., Peterseim, D.: The -method for the Gross–Pitaevskii eigenvalue problem. Numer. Math. 148, 575–610 (2021). 10.1007/s00211-021-01216-5
6. Altmann R Henning P Peterseim D Localization and delocalization of ground states of Bose–Einstein condensates under disorder SIAM J. Appl. Math. 2022 82 1 330 358 10.1137/20M1342434
Altmann, R., Henning, P., Peterseim, D.: Localization and delocalization of ground states of Bose–Einstein condensates under disorder. SIAM J. Appl. Math. 82(1), 330–358 (2022). 10.1137/20M1342434
7. Altmann R Peterseim D Localized computation of eigenstates of random Schrödinger operators SIAM J. Sci. Comput. 2019 41 B1211 B1227 10.1137/19M1252594
Altmann, R., Peterseim, D.: Localized computation of eigenstates of random Schrödinger operators. SIAM J. Sci. Comput. 41, B1211–B1227 (2019). 10.1137/19M1252594
8. Altmann R Peterseim D Stykel T Energy-adaptive Riemannian optimization on the Stiefel manifold ESAIM Math. Model. Numer. Anal. 2022 56 5 1629 1653 10.1051/m2an/2022036
Altmann, R., Peterseim, D., Stykel, T.: Energy-adaptive Riemannian optimization on the Stiefel manifold. ESAIM Math. Model. Numer. Anal. 56(5), 1629–1653 (2022). 10.1051/m2an/2022036
9. Bai Z Li RC Lu D Sharp estimation of convergence rate for self-consistent field iteration to solve eigenvector-dependent nonlinear eigenvalue problems SIAM J. Matrix Anal. Appl. 2022 43 1 301 327 10.1137/20M136606X
Bai, Z., Li, R.C., Lu, D.: Sharp estimation of convergence rate for self-consistent field iteration to solve eigenvector-dependent nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl. 43(1), 301–327 (2022). 10.1137/20M136606X
10. Bao W Chern IL Lim FY Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose–Einstein condensates J. Comput. Phys. 2006 219 2 836 854 10.1016/j.jcp.2006.04.019
Bao, W., Chern, I.L., Lim, F.Y.: Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose–Einstein condensates. J. Comput. Phys. 219(2), 836–854 (2006). 10.1016/j.jcp.2006.04.019
11. Bao W Du Q Computing the ground state solution of Bose–Einstein condensates by a normalized gradient flow SIAM J. Sci. Comput. 2004 25 5 1674 1697 10.1137/S1064827503422956
Bao, W., Du, Q.: Computing the ground state solution of Bose–Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput. 25(5), 1674–1697 (2004). 10.1137/S1064827503422956
12. Bao W Tang W Ground-state solution of Bose–Einstein condensate by directly minimizing the energy functional J. Comput. Phys. 2003 187 1 230 254 10.1016/S0021-9991(03)00097-4
Bao, W., Tang, W.: Ground-state solution of Bose–Einstein condensate by directly minimizing the energy functional. J. Comput. Phys. 187(1), 230–254 (2003). 10.1016/S0021-9991(03)00097-4
13. Caliari M Ostermann A Rainer S Thalhammer M A minimisation approach for computing the ground state of Gross–Pitaevskii systems J. Comput. Phys. 2009 228 2 349 360 10.1016/j.jcp.2008.09.018
Caliari, M., Ostermann, A., Rainer, S., Thalhammer, M.: A minimisation approach for computing the ground state of Gross–Pitaevskii systems. J. Comput. Phys. 228(2), 349–360 (2009). 10.1016/j.jcp.2008.09.018
14. Cancès, E.: SCF algorithms for HF electronic calculations. In: Mathematical Models and Methods for Ab Initio Quantum Chemistry, Lecture Notes in Chemistry, vol. 174, pp. 17–43. Springer, Berlin, Heidelberg (2000). 10.1007/978-3-642-57237-1_2
15. Cancès E Self-consistent field algorithms for Kohn–Sham models with fractional occupation numbers J. Chem. Phys. 2001 114 24 10616 10622 10.1063/1.1373430
Cancès, E.: Self-consistent field algorithms for Kohn–Sham models with fractional occupation numbers. J. Chem. Phys. 114(24), 10616–10622 (2001). 10.1063/1.1373430
16. Cancès E Chakir R Maday Y Numerical analysis of nonlinear eigenvalue problems J. Sci. Comput. 2010 45 90 117 10.1007/s10915-010-9358-1
Cancès, E., Chakir, R., Maday, Y.: Numerical analysis of nonlinear eigenvalue problems. J. Sci. Comput. 45, 90–117 (2010). 10.1007/s10915-010-9358-1
17. Cancès E Chakir R Maday Y Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models ESAIM Math. Model. Numer. Anal. 2012 46 2 341 388 10.1051/m2an/2011038
Cancès, E., Chakir, R., Maday, Y.: Numerical analysis of the planewave discretization of some orbital-free and Kohn-Sham models. ESAIM Math. Model. Numer. Anal. 46(2), 341–388 (2012). 10.1051/m2an/2011038
18. Cancès E Kemlin G Levitt A Convergence analysis of direct minimization and self-consistent iterations SIAM J. Matrix Anal. Appl. 2021 42 1 243 274 10.1137/20M1332864
Cancès, E., Kemlin, G., Levitt, A.: Convergence analysis of direct minimization and self-consistent iterations. SIAM J. Matrix Anal. Appl. 42(1), 243–274 (2021). 10.1137/20M1332864
19. Cancès E Le Bris C On the convergence of SCF algorithms for the Hartree–Fock equations ESAIM Math. Model. Numer. Anal. 2000 34 4 749 774 10.1051/m2an:2000102
Cancès, E., Le Bris, C.: On the convergence of SCF algorithms for the Hartree–Fock equations. ESAIM Math. Model. Numer. Anal. 34(4), 749–774 (2000). 10.1051/m2an:2000102
20. Chen Z Lu J Lu Y Zhang X On the convergence of Sobolev gradient flow for the Gross–Pitaevskii eigenvalue problem SIAM J. Numer. Anal. 2024 62 2 667 691 10.1137/23M1552553
Chen, Z., Lu, J., Lu, Y., Zhang, X.: On the convergence of Sobolev gradient flow for the Gross–Pitaevskii eigenvalue problem. SIAM J. Numer. Anal. 62(2), 667–691 (2024). 10.1137/23M1552553
21. Dai X Wang Q Zhou A Gradient flow based Kohn–Sham density functional theory model Multiscale Model. Simul. 2020 18 4 1621 1663 10.1137/19M1276170
Dai, X., Wang, Q., Zhou, A.: Gradient flow based Kohn–Sham density functional theory model. Multiscale Model. Simul. 18(4), 1621–1663 (2020). 10.1137/19M1276170
22. Danaila I Protas B Computation of ground states of the Gross–Pitaevskii functional via Riemannian optimization SIAM J. Sci. Comput. 2017 39 6 B1102 B1129 10.1137/17M1121974
Danaila, I., Protas, B.: Computation of ground states of the Gross–Pitaevskii functional via Riemannian optimization. SIAM J. Sci. Comput. 39(6), B1102–B1129 (2017). 10.1137/17M1121974
23. Dion CM Cancès E Ground state of the time-independent Gross–Pitaevskii equation Comput. Phys. Comm. 2007 177 10 787 798 10.1016/j.cpc.2007.04.007
Dion, C.M., Cancès, E.: Ground state of the time-independent Gross–Pitaevskii equation. Comput. Phys. Comm. 177(10), 787–798 (2007). 10.1016/j.cpc.2007.04.007
24. Du CE Liu CS Newton-Noda iteration for computing the ground states of nonlinear Schrödinger equations SIAM J. Sci. Comput. 2022 44 4 A2370 A2385 10.1137/21M1435793
Du, C.E., Liu, C.S.: Newton-Noda iteration for computing the ground states of nonlinear Schrödinger equations. SIAM J. Sci. Comput. 44(4), A2370–A2385 (2022). 10.1137/21M1435793
25. Edelman A Arias TA Smith ST The geometry of algorithms with orthogonality constraints SIAM J. Matrix Anal. Appl. 1998 20 2 303 353 10.1137/S0895479895290954
Edelman, A., Arias, T.A., Smith, S.T.: The geometry of algorithms with orthogonality constraints. SIAM J. Matrix Anal. Appl. 20(2), 303–353 (1998). 10.1137/S0895479895290954
26. García-Ripoll JJ Pérez-García VM Optimizing Schrödinger functionals using Sobolev gradients: applications to quantum mechanics and nonlinear optics SIAM J. Sci. Comput. 2001 23 4 1316 1334 10.1137/S1064827500377721
García-Ripoll, J.J., Pérez-García, V.M.: Optimizing Schrödinger functionals using Sobolev gradients: applications to quantum mechanics and nonlinear optics. SIAM J. Sci. Comput. 23(4), 1316–1334 (2001). 10.1137/S1064827500377721
27. Henning P The dependency of spectral gaps on the convergence of the inverse iteration for a nonlinear eigenvector problem Math. Models Methods Appl. Sci. 2023 33 07 1517 1544 10.1142/S0218202523500343
Henning, P.: The dependency of spectral gaps on the convergence of the inverse iteration for a nonlinear eigenvector problem. Math. Models Methods Appl. Sci. 33(07), 1517–1544 (2023). 10.1142/S0218202523500343
28. Henning, P., Jarlebring, E.: The Gross–Pitaevskii equation and eigenvector nonlinearities: numerical methods and algorithms. Preprint (2022)
29. Henning P Peterseim D Sobolev gradient flow for the Gross–itaevskii eigenvalue problem: global convergence and computational efficiency SIAM J. Numer. Anal. 2020 58 3 1744 1772 10.1137/18M1230463
Henning, P., Peterseim, D.: Sobolev gradient flow for the Gross–itaevskii eigenvalue problem: global convergence and computational efficiency. SIAM J. Numer. Anal. 58(3), 1744–1772 (2020). 10.1137/18M1230463
30. Henning P Yadav M On discrete ground states of rotating Bose–Einstein condensates Math. Comput. 2024 10.1090/mcom/3962
Henning, P., Yadav, M.: On discrete ground states of rotating Bose–Einstein condensates. Math. Comput. (2024). 10.1090/mcom/3962
31. Hohenberg P Kohn W Inhomogeneous electron gas Phys. Rev. 1964 136 B864 B871 10.1103/PhysRev.136.B864
Hohenberg, P., Kohn, W.: Inhomogeneous electron gas. Phys. Rev. 136, B864–B871 (1964). 10.1103/PhysRev.136.B864
32. Hu G Wang T Zhou J A linearized structure-preserving numerical scheme for a gradient flow model of the Kohn–Sham density functional theory East Asian J. Appl. Math. 2023 13 2 299 319 10.4208/eajam.2022-134.081022
Hu, G., Wang, T., Zhou, J.: A linearized structure-preserving numerical scheme for a gradient flow model of the Kohn–Sham density functional theory. East Asian J. Appl. Math. 13(2), 299–319 (2023). 10.4208/eajam.2022-134.081022
33. Jarlebring E Kvaal S Michiels W An inverse iteration method for eigenvalue problems with eigenvector nonlinearities SIAM J. Sci. Comput. 2014 36 4 A1978 A2001 10.1137/S1064827500366124
Jarlebring, E., Kvaal, S., Michiels, W.: An inverse iteration method for eigenvalue problems with eigenvector nonlinearities. SIAM J. Sci. Comput. 36(4), A1978–A2001 (2014). 10.1137/S1064827500366124
34. Jarlebring E Upadhyaya P Implicit algorithms for eigenvector nonlinearities Numer. Algorithms 2022 90 301 321 10.1007/s11075-021-01189-4
Jarlebring, E., Upadhyaya, P.: Implicit algorithms for eigenvector nonlinearities. Numer. Algorithms 90, 301–321 (2022). 10.1007/s11075-021-01189-4
35. Jiao S Zhang Z Wu K Wan L Ma H Li J Chen S Qin X Liu J Ding Z Yang J Li Y Hu W Lin L Yang C KSSOLV 2.0: An efficient MATLAB toolbox for solving the Kohn–Sham equations with plane-wave basis set Comput. Phys. Comm. 2022 279 108424 10.1016/j.cpc.2022.108424
Jiao, S., Zhang, Z., Wu, K., Wan, L., Ma, H., Li, J., Chen, S., Qin, X., Liu, J., Ding, Z., Yang, J., Li, Y., Hu, W., Lin, L., Yang, C.: KSSOLV 2.0: An efficient MATLAB toolbox for solving the Kohn–Sham equations with plane-wave basis set. Comput. Phys. Comm. 279, 108424 (2022). 10.1016/j.cpc.2022.108424
36. Kazemi P Eckart M Minimizing the Gross–Pitaevskii energy functional with the Sobolev gradient—Analytical and numerical results Int. J. Comput. Methods 2010 7 3 453 475 10.1142/S0219876210002301
Kazemi, P., Eckart, M.: Minimizing the Gross–Pitaevskii energy functional with the Sobolev gradient—Analytical and numerical results. Int. J. Comput. Methods 7(3), 453–475 (2010). 10.1142/S0219876210002301
37. Kohn W Sham LJ Self-consistent equations including exchange and correlation effects Phys. Rev. 1965 140 A1133 A1138 10.1103/PhysRev.140.A1133
Kohn, W., Sham, L.J.: Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133–A1138 (1965). 10.1103/PhysRev.140.A1133
38. Le Bris C Computational chemistry from the perspective of numerical analysis Acta Numer. 2005 14 363 444 10.1017/S096249290400025X
Le Bris, C.: Computational chemistry from the perspective of numerical analysis. Acta Numer. 14, 363–444 (2005). 10.1017/S096249290400025X
39. Lieb EH Seiringer R Yngvason J A rigorous derivation of the Gross–Pitaevskii energy functional for a two-dimensional Bose gas Commun. Math. Phys. 2001 224 1 17 31 10.1007/s002200100533
Lieb, E.H., Seiringer, R., Yngvason, J.: A rigorous derivation of the Gross–Pitaevskii energy functional for a two-dimensional Bose gas. Commun. Math. Phys. 224(1), 17–31 (2001). 10.1007/s002200100533
40. Liu X Wen Z Wang X Ulbrich M Yuan Y On the analysis of the discretized Kohn–Sham density functional theory SIAM J. Numer. Anal. 2015 53 4 1758 1785 10.1137/140957962
Liu, X., Wen, Z., Wang, X., Ulbrich, M., Yuan, Y.: On the analysis of the discretized Kohn–Sham density functional theory. SIAM J. Numer. Anal. 53(4), 1758–1785 (2015). 10.1137/140957962
41. Perdew J Wang Y Accurate and simple analytic representation of the electron-gas correlation energy Phys. Rev. B 1992 45 23 13244 13249 10.1103/PhysRevB.45.13244
Perdew, J., Wang, Y.: Accurate and simple analytic representation of the electron-gas correlation energy. Phys. Rev. B 45(23), 13244–13249 (1992). 10.1103/PhysRevB.45.13244
42. Perdew JP Zunger A Self-interaction correction to density-functional approximations for many-electron systems Phys. Rev. B 1981 23 5048 5079 10.1103/PhysRevB.23.5048
Perdew, J.P., Zunger, A.: Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 23, 5048–5079 (1981). 10.1103/PhysRevB.23.5048
43. Pitaevskii LP Stringari S Bose–Einstein Condensation 2003 Oxford Oxford University Press
Pitaevskii, L.P., Stringari, S.: Bose–Einstein Condensation. Oxford University Press, Oxford (2003)
44. Raza N Sial S Siddiqi SS Lookman T Energy minimization related to the nonlinear Schrödinger equation J. Comput. Phys. 2009 228 7 2572 2577 10.1016/j.jcp.2008.12.016
Raza, N., Sial, S., Siddiqi, S.S., Lookman, T.: Energy minimization related to the nonlinear Schrödinger equation. J. Comput. Phys. 228(7), 2572–2577 (2009). 10.1016/j.jcp.2008.12.016
45. Roothaan CCJ New developments in molecular orbital theory Rev. Mod. Phys. 1951 23 69 89 10.1103/RevModPhys.23.69
Roothaan, C.C.J.: New developments in molecular orbital theory. Rev. Mod. Phys. 23, 69–89 (1951). 10.1103/RevModPhys.23.69
46. Sato H Aihara K Cholesky QR-based retraction on the generalized Stiefel manifold Comput. Optim. Appl. 2019 72 2 293 308 10.1007/s10589-018-0046-7
Sato, H., Aihara, K.: Cholesky QR-based retraction on the generalized Stiefel manifold. Comput. Optim. Appl. 72(2), 293–308 (2019). 10.1007/s10589-018-0046-7
47. Schneider R Rohwedder T Neelov A Blauert J Direct minimization for calculating invariant subspaces in density functional computations of the electronic structure J. Comput. Math. 2009 27 2–3 360 387
Schneider, R., Rohwedder, T., Neelov, A., Blauert, J.: Direct minimization for calculating invariant subspaces in density functional computations of the electronic structure. J. Comput. Math. 27(2–3), 360–387 (2009)
48. Sherman J Morrison WJ Adjustment of an inverse matrix corresponding to a change in one element of a given matrix Ann. Math. Statist. 1950 21 1 124 127 10.1214/aoms/1177729893
Sherman, J., Morrison, W.J.: Adjustment of an inverse matrix corresponding to a change in one element of a given matrix. Ann. Math. Statist. 21(1), 124–127 (1950). 10.1214/aoms/1177729893
49. Teter MP Payne MC Allan DC Solution of Schrödinger’s equation for large systems Phys. Rev. B 1989 40 12255 12263 10.1103/PhysRevB.40.12255
Teter, M.P., Payne, M.C., Allan, D.C.: Solution of Schrödinger’s equation for large systems. Phys. Rev. B 40, 12255–12263 (1989). 10.1103/PhysRevB.40.12255
50. Wloka J Partial Differential Equations 1987 Cambridge Cambridge University Press
Wloka, J.: Partial Differential Equations. Cambridge University Press, Cambridge (1987)
51. Yang C Meza JC Lee B Wang LW KSSOLV—a MATLAB toolbox for solving the Kohn–Sham equations ACM Trans. Math. Softw. 2009 36 2 1 35 10.1145/1499096.1499099
Yang, C., Meza, J.C., Lee, B., Wang, L.W.: KSSOLV—a MATLAB toolbox for solving the Kohn–Sham equations. ACM Trans. Math. Softw. 36(2), 1–35 (2009). 10.1145/1499096.1499099
52. Yang C Meza JC Wang LW A constrained optimization algorithm for total energy minimization in electronic structure calculation J. Comput. Phys. 2006 217 2 709 721 10.1016/j.jcp.2006.01.030
Yang, C., Meza, J.C., Wang, L.W.: A constrained optimization algorithm for total energy minimization in electronic structure calculation. J. Comput. Phys. 217(2), 709–721 (2006). 10.1016/j.jcp.2006.01.030
53. Zhang Z Exponential convergence of Sobolev gradient descent for a class of nonlinear eigenproblems Commun. Math. Sci. 2022 20 377 403 10.4310/CMS.2022.v20.n2.a4
Zhang, Z.: Exponential convergence of Sobolev gradient descent for a class of nonlinear eigenproblems. Commun. Math. Sci. 20, 377–403 (2022). 10.4310/CMS.2022.v20.n2.a4
54. Zhao Z Bai ZJ Jin XQ A Riemannian Newton algorithm for nonlinear eigenvalue problems SIAM J. Matrix Anal. Appl. 2015 36 2 752 774 10.1137/140967994
Zhao, Z., Bai, Z.J., Jin, X.Q.: A Riemannian Newton algorithm for nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl. 36(2), 752–774 (2015). 10.1137/140967994
