# SELF-CONSISTENT GRADIENT-LIKE EIGEN DECOM## POSITION IN SOLVING SCHR ¨ODINGER EQUATIONS

**Anonymous authors**
Paper under double-blind review

ABSTRACT

The Schr¨odinger equation is at the heart of modern quantum mechanics. Since
exact solutions of the ground state are typically intractable, standard approaches
approximate Schr¨odinger’s equation as forms of nonlinear generalized eigenvalue problems F (V )V = SV Λ in which F (V ), the matrix to be decomposed, is a function of its own top-k smallest eigenvectors V, leading to a “selfconsistency problem”. Traditional iterative methods heavily rely on high-quality
initial guesses of V generated via domain-specific heuristics methods based on
quantum mechanics. In this work, we eliminate such a need for domain-specific
heuristics by presenting a novel framework, Self-consistent Gradient-like Eigen
Decomposition (SCGLED) that regards F (V ) as a special “online data generator”, thus allows gradient-like eigendecomposition methods in streaming k-PCA
to approach the self-consistency of the equation from scratch in an iterative way
similar to online learning. With several critical numerical improvements, SCGLED is robust to initial guesses, free of quantum-mechanism-based heuristics
designs, and neat in implementation. Our experiments show that it not only can
simply replace traditional heuristics-based initial guess methods with large performance advantage (achieved averagely 25x more precise than the best baseline in
similar wall time), but also is capable of finding highly precise solutions independently without any traditional iterative methods.

1 INTRODUCTION

While the many-body Schr¨odinger equation governs all quantum-mechanical systems, it can be
solved analytically only for extremely simple ones like an isolated hydrogen atom. Nonetheless, a
series of important approximation theories are developed so that it can be solved numerically for
many relevant molecules and materials including tens to thousands of electrons, making it possible
for physicists and chemists to explore the properties of matters in an ab-initio or first-principle way,
serving as a foundation of computational physics, chemistry and material science research.

Among those approximation theories, two mainstream directions involve the Hartree-Fock (Hartree,
1928; Hartree & Hartree, 1935) and density functional theories (Hohenberg & Kohn, 1964) that
either approximate the Schr¨odinger equation as a “Hartree-Fock” equation or as a “Kohn-Sham
equation” (Kohn & Sham, 1965). While very different in nature, both approximated equations are
finally formed as a nonlinear, generalized eigenvalue problem F (V )V = SV Λ in which F (V ) and
_S are N_ _N real symmetric matrices, Λ = diag(λ1,_ _, λk) is a k_ _k matrix containing the_
_×_ _· · ·_ _×_
top-k smallest eigenvalues, and V = [v1, _, vk] is an N_ _k matrix containing the corresponding_
_· · ·_ _×_
top-k eigenvectors. Interestingly, the matrix F which we aim to decompose is explicitly defined as
a given function of the eigenvectors V (denoted as F (V ) : R[N] _[×][k]_ _→_ R[N] _[×][N]_ ), leading to a “causeand-effect dilemma” or “self-consistency problem”. That is, the form of the equation is defined by
its final solution, but how does one know the solution if the form of the equation is not determined?

The dominant method in solve such equations is generally referred to as Self-Consistent Fields
(SCFs). SCFs are the main component of nearly all mainstream quantum chemistry and solid-state
physics software. The general idea is to perform a fixed-point iteration, that is, guess an initial
solution V0 from scratch, then solve F (Vk 1)Vk = SVkΛk repeatedly to obtain eigenvectors Vk at
_−_
the kth iteration, k = 1, 2, · · ·, until ∥F (Vk) _−_ _F_ (Vk−1)∥ is less than a given convergence threshold.
Unfortunately, there is no guarantee that such an iteration can always lead to a converged solution


-----

(Froese Fischer, 1987). Actually, SCFs are extremely sensitive to the selection of initial guess V0.
Without a high-quality initialisation, SCFs easily fail and oscillate between two or more different
results.

Till now, all popular choices for the initial guess are heuristics-based, highly reliant on specialized
prior knowledge of quantum mechanics. For example, the SAD method (Van Lenthe et al., 2006)
applies the physical intuition that molecular electron densities can be approximated via the sum of
atomic electron densities. Here, the core Hamiltonian method ignores all interactions between electrons, which simplifies the equation to a standard eigen-decomposition problem. More methods and
details can be found in (Lehtola, 2019). Besides the prerequisite of domain-specific theories, such
methods are also sophisticated in implementation to achieve satisfactory convergence performance,
which is a main barrier for practitioners to build high-performance solvers from scratch.

Meanwhile, in the machine learning community, researchers are facing the challenge of analyzing
real-time, online data represented as streams. In contrast to the normal case that data is fully available offline, online data is received sequentially and may even change over time, thus specialized
methods are developed to dynamically adapt to new patterns in the data stream, which is called
stochastic, streaming or online algorithms. Specially, a series of works (Oja, 1982; Tang, 2019;
Allen-Zhu & Li, 2017; Gemp et al., 2021) focus on stochastic k-PCA, which aim to estimate the
top-k principal components of a data stream in a real-time manner. To handle the stochasticity
of data, such stochastic k-PCA methods usually contain gradient-like eigen decomposition. When
a new data sample is available, the eigenvectors will be updated towards a gradient-like direction
computed by the new sample, with a small learning rate.

In this work, we noticed that F (V ) in the aforementioned eigen-decomposition problem is sub
ject to concurrent update during the sequential update of V towards self-consistency: V1 _F =F (V1)_

_V2_ _F =F (V2)_ . In this case F (V ) can be regarded as a special “online data generator”, which−−−−−−→
_−−−−−−→· · ·_
shares similarities with the online learning setting that online data is concurrently updated during
the sequential learning process. In this way, we explore the possibility of applying gradient-like
eigen-decomposition in stochastic k-PCA methods to handle the self-consistency of approximated
Schr¨odinger equations. Together with several numerical improvements to enhance the smoothness
of optimization, we developed Self-consistent Gradient-like Eigen Decomposition (SCGLED), an
efficient solver for the self-consistent approximated Schr¨odinger equations, SCGLED is robust to
initial guesses V0, free of quantum-mechanism-based heuristics design, and neat in implementation.
While it can simply replace traditional heuristics-based initial guess methods with performance advantage, it is also capable of finding highly precise solutions without any traditional SCF iteration.

2 RELATED WORK

The computational theories for quantum many-body systems, especially for determining the wave
function of Schr¨odinger equations, has a long history starting from 1920s. There are three mainstream theories: Hartree-Fock theory (Hartree, 1928; Hartree & Hartree, 1935), density functional
theory (Hohenberg & Kohn, 1964) and quantum Monte Carlo. With the prosperity of deep learning
and differentiable optimization in recent years, there are a series of works focusing on deep-learningaided wave function representation of quantum Monte Carlo methods such as FermiNet (Pfau et al.,
2020) and PauliNet (Hermann et al., 2020). For density functional theory, there are also some works
using neural networks to learn the exchange-correlation functional (Li et al., 2021; Kasim & Vinko,
2021). However, these works focus more on improving simulation accuracy towards physical reality by using neural networks as a better functional approximator, while our work’s focus is very
different, stressing on the optimization efficiency while the equation is completely given. There
are also works on direct optimization for Kohn-Sham equations by total energy minimization (Yang
et al., 2006), which model the problem as a constrained optimization problem instead, and more
domain-specific knowledge is involved, while our work tries to solve the equation from a purely
optimization-based aspect.

For gradient-like eigen-decomposition, most of the works focus on stochastic k-PCA, which estimates the top-k principal components of a data stream. Let x ∈ R[d] denote a random data sample at
time step t, vi ∈ R[d], i = 1, 2, · · ·, k denote the estimate of ith principal component (eigenvector of


-----

𝐹(𝑉)


Generalized Eigendecomposition

“Self-consistent”

Definition of 𝐹(𝑉)

|𝜆1|Col2|
|---|---|
||𝜆2|


𝜆1 < 𝜆2 < ⋯

|Col1|Col2|
|---|---|
|𝑣|𝑣|
|1|2|
|||

|Col1|Col2|
|---|---|
|𝑣|𝑣|
|1|2|
|||

|Col1|Col2|Col3|Col4|
|---|---|---|---|
|||||
||𝐹(|𝑉|)|
|||||

|Col1|Col2|Col3|Col4|
|---|---|---|---|
|||||
|||𝐻||
|||||

|Col1|Col2|Col3|Col4|
|---|---|---|---|
||2𝑉|𝑉|⊤|
|||||
|||||


Figure 1: An visualization of the problem. S, H ∈ R[N] _[×][N]_ and Ueff : R[N] _[×][N]_ _→_ R[N] _[×][N]_ are given
input. The solution V _[∗]_ should obey both the generalized eigen-decomposition and the definition of
_F_ (V ), which leads to a “self-consistent” problem.

E[XX _[⊤]]) at time t, corresponding to the ith largest eigenvalue. η denotes the learning rate._

Oja’s algorithm: vi[′] [+][ η][(][xx][⊤][v][i][)][,][ ∀][i][ = 1][,][ · · ·][, k][ and][ (][v][1][,][ · · ·][, v][k][)][ ←] _[QR][(][v]1[′]_ _[,][ · · ·][, v]k[′]_ [)]

_[←]_ _[v][i]_

Here QR(v1, · · ·, vk) is the Gram-Schmidt decomposition that orthonormalizes v1[t] _[,][ · · ·][, v]k[t]_ [. There]
are also several follow-up works (Sanger, 1989; Tang, 2019) that improve the performance, in which
a most recent one is EigenGame (Gemp et al., 2021)


_x[⊤]vi, x[⊤]vj_

EigenGame: _vi_ 2x _x[⊤]vi_ _⟨_ _⟩_ _,_ _i = 1,_ _, k_
_∇_ _←_ h _−_ Xj<i _⟨x[⊤]vj, x[⊤]vj⟩_ _[x][⊤][v][j]i_ _∀_ _· · ·_

_vi_ _vi_ _vi, vi_ _vi, vi[′]_ [+][ η][∇][R][v][i] [and][ v][i] _vi[′]_
_∇[R]_ _←∇_ _−⟨∇_ _⟩_ _[←]_ _[v][i]_ _[←]_ _∥vi∥_ _[,][ ∀][i][ = 1][,][ · · ·][, k.]_


EigenGame: _vi_ 2x _x[⊤]vi_
_∇_ _←_ _−_
h


_j<i_


3 PROBLEM DESCRIPTION

In this paper, while what we want to solve is the approximated Schr¨odinger Equations, whose physical background is elaborated in Appendix A, we abstract it mathematically as the following type of
nonlinear generalized eigenvalue problem[1]

_F_ (V )V = SV Λ, (1)

where

-  F (V ): an N × N real symmetric matrix to be decomposed, which is defined in equation 2
as a function of V .

-  S: an N × N positive semi-definite matrix, which is a constant input in the problem.

-  Λ: Λ = diag(λ1, _, λk) is a k_ _k diagonal matrix containing the top-k smallest eigen-_
_· · ·_ _×_
values.

-  V : V = [v1, _, vk] is an N_ _k matrix containing k column eigenvectors corresponding_
_· · ·_ _×_
to the top-k smallest eigenvalues.

The definition of F (V ) : R[N] _[×][k]_ _→_ R[N] _[×][N]_ is as follows:

def
_F_ (V ) = H + Ueff(2V V ), (2)

_[⊤]_

in which H is an N × N real symmetric matrix, which is given in this problem. Ueff : R[N] _[×][N]_ _→_
R[N] _[×][N]_ is a given function. We also define P (V ) = 2V V _[⊤]_ for convenience. To conclude, the input

1While the regular eigenvalue problem can be described as finding V that obeys AV = V Λ in which A is
an N × N matrix to be decomposed, Λ = diag(λ1, · · ·, λN ) contains all eigenvalues and V = (v1, · · ·, vN )
contains the corresponding eigenvectors, a generalized eigenvalue problem adds an additional matrix B with
the form AV = BV Λ. When B = I, it degenerates to a regular eigenvalue problem.


-----

of the problem is S, H, Ueff( ) and k, and the output of the problem is the eigenvectors V that

_·_ _[∗]_
obeys both equation 1 and equation 2. The top-k smallest eigenvalues of F (V _[∗]) are guaranteed to_
be negative.

To make it more clear, a toy example is provided as follows:


0.6953 1.1204 0.9584

_, H =_ _−_ _−_

1.0  −0.9584 _−1.1204_

_PλσEuvλσ_ _PλσEuλσv,_
_λ,σ_ _−_ 2[1] _λ,σ_

X X


1.0 0.6953
_S =_
0.6953 1.0



_, k = 1_



[Ueff(P )]uv =


0.7746 0.4441
_E11 =_
0.4441 0.5697



0.4441 0.2970
_, E12 = E21 =_
0.2970 0.4441



0.5697 0.4441
_, E22 =_
0.4441 0.7746



in which P = 2V V _[⊤], Ueff(P_ ) ∈ R[2][×][2], E is a 2 × 2 × 2 × 2 tensor. The solution of the toy example

0.3655 0.5939
should be V _[∗]_ = (0.5489, 0.5489)[⊤], in this case F (V _[∗]) =_ −−0.5939 _−−0.3655, and the result of_

top-1 eigen-decomposition

0.3655 0.5939 1.0 0.6953
_−_ _−_ _V =_ _V Λ_
−0.5939 _−0.3655_ 0.6953 1.0 

will happen to be exactly V _[∗]_ = (0.5489, 0.5489)[⊤] with the smallest eigenvalue λ1 = −0.5782,
while the other one is [1.2115, −1.2115][⊤] with eigenvalue λ2 = 0.6703.

4 SOLVING APPROXIMATED SCHR ¨ODINGER EQUATIONS USING
SELF-CONSISTENT GRADIENT-LIKE EIGEN DECOMPOSITION

First, notice that the stochastic k-PCA methods usually can be generalized to decompose a real
symmetric matrix M . For Oja’s algorithm (Oja & Karhunen, 1985), that is

_vi[′]_ [+][ ηMv][i] _i = 1,_ _, k and (v1,_ _, vk)_ _QR(v1[′]_ _[,][ · · ·][, v]k[′]_ [)][,] (3)

_[←]_ _[v][i]_ _∀_ _· · ·_ _· · ·_ _←_

in which vi[t] [is the eigenvector corresponding to the][ i][th largest eigenvalue.]

Due to the lack of efficiency, such decomposition methods are not favorable for traditional eigendecomposition (M is a given constant input) compared with classical algorithms such as QR iteration. However, as stochastic algorithms, they have a unique advantage that they can adapt to dynamic
change of M, which suits the sticking point in equation 1 that F (V ) is subject to concurrent update
when V changes. This motivates us to apply them to tackle the approximated Schr¨odinger equations.

Together with the orthogonalization technique introduced in Appendix B which transforms the generalized eigenvalue problem into a standard one, we replace M in equation 3 with F _[′]_ = X _[⊤]F_ (V )X,
and update F _[′]_ in each time step immediately after each update of V to maintain the self-consistency,
which forms the initial version of our proposed algorithm, self-consistent gradient-like eigen decomposition (SCGLED), shown in Algorithm 1. Note that we need the eigenvectors corresponding to
top-k smallest eigenvalues rather than the largest ones, so we decompose −F _[′]_ instead of F _[′]. Dif-_
ferent from traditional SCF method shown in Appendix B where Vt[′] [is purely a temporary variable]
that should be discarded after each iteration, Vt[′] [in our proposed algorithm is more likely a “training]
variable” that is first randomly initialized and then trained during the whole iterative process.

While the convergence analysis of Algorithm 1 is heavily blocked by the intrinsic nonlinearity of
approximated Schrodinger equations, which is fundamentally complex in computational physics, we
instead use heuristics to boost the empirical convergence performance in later paragraphs, and prove
the correctness of our algorithm under some robustness assumption of F (V ) on the convergence
point, with the following proposition
**Proposition 1. If Algorithm 1 converges to a stable convergence point V** _[∗]_ _and for a small pertur-_
_bation ϵ towards V_ _[∗], [X_ _[⊤]F_ (V _[∗])X −_ _X_ _[⊤]F_ (V _[∗]_ + ϵ)X](V _[∗]_ + ϵ) = O(ϵ[2]), then V _[∗]_ _is the solution_
_of equation 1._

The proof is deferred to Appendix C.


-----

**Algorithm 1 The vanilla version of self-consistent gradient-like eigen decomposition (SCGLED)**

**Input: H, S, Ueff(·) in equation 1 and equation 2, learning rate η > 0**
**Output: V** _[∗], the solution of equation 1_

1: Initialize V _[′]_ _∈_ R[N] _[×][k]_ randomly.
2: Find X satisfying X _[⊤]SX = I._
3: while V _[′]_ is not converged do
4: _V ←_ _XV_ _[′]_

5: _F_ _[′]_ _←_ _X_ _[⊤]F_ (V )X following equation 2

6: Update V _[′]_ for one-step decomposition of −F _[′]_ using equation 3.

7: end while
8: V _[∗]_ _←_ _XV_ _[′]_


However, this initial version is of poor efficiency and stability. The computational efficiency is
blocked by the Gram-Schmidt process which is highly order-preserving thus cannot be decentralized/vectorized. Also, due to the existence of F (V ) = H +Ueff(P (V )) that is sensitive to the change
of V, the iteration process is highly nonlinear. In this case oscillation will easily occur which hinders
the iteration process toward convergence. To tackle these problems, we made several improvements
to Algorithm 1 as follows.

First, we replace the classical Oja’s algorithm with the decentralized version of EigenGame (Gemp
et al., 2021). EigenGame replaces the time-consuming QR decomposition by a generalized GramSchmidt step in its gradient _vi, which makes it theoretically harder to analyse (since the orthonor-_
_∇_
mality of eigenvectors is not guaranteed during the iteration), but with better empirical performance
due to its decentralized nature. Its procedure is as follows:


_vi[⊤][Mv][j]_

_vj_ _,_ _vi_ _vi_ _vi, vi_ _vi_ _i = 1,_ _, k_ (4)
_vj[⊤][Mv][j]_ _∇[R]_ _←∇_ _−⟨∇_ _⟩_ _∀_ _· · ·_
i


_vi_ 2M _vi_
_∇_ _←_ _−_
h


_j<i_


_vi[′]_ _[←]_ _[v][i]_ [+][ η][∇][R][v][i][,][ and][ v][i] _[←]_ _∥vvii[′][′][∥]_ _∀i = 1, · · ·, k_ (5)

Second, to reduce oscillation during the iteration process, we introduce the damping technique for
the update of F . That is,
_Ft = (1_ _α)Ft_ 1 + αF (Vt), (6)
_−_ _−_
in which α is the mixing hyperparameter between 0 and 1. Here we set it as 0.2.

Third, the selection of learning rate η is highly tricky for different molecule inputs. To enhance
the robustness of the algorithm towards learning rate, we introduce the momentum method for the
update of V _[′]_ as follows:

_mt = βmt−1 + η∇V_ _[′],_ (7)

_V_ = V + mt, (8)

_[′]_ _[′]_

in which the momentum term β is set to 0.9.

Fourth, instead of updating F in every iteration, we control the update interval of F via a parameter IF . The first reason is for efficiency, since V _[′]_ only update slightly in each time step with a
small learning rate η, it may not be necessary to do a fresh computation of F (V ) in each time step,
especially considering that the computation of Ueff(P (V )) is relatively time-consuming. The second reason is for stability, since the damping technique in equation 6 will degenerate if previously
updated Ft 1 and current F (V ) are too close.
_−_

Summarizing all considerations above, our proposed algorithm is shown in Algorithm 2

5 EXPERIMENTS

In this section, we perform extensive performance benchmarks on W4-17 dataset(Karton et al.,
2017), which is also applied in prior benchmark work (Lehtola, 2019). All the 160 singlet molecules
in the W4-17 dataset are used to evaluate the performance of our proposed algorithm. Our algorithm


-----

**Algorithm 2 Self-consistent gradient-like eigen decomposition (SCGLED)**

**Input: H, S, Ueff(·) in equation 1 and equation 2, total number of iterations T**, learning rate η, F ’s
update ratio IF

**Output: V** _[∗], the solution of equation 1_

1: Initialize V _[′]_ _∈_ R[N] _[×][k]_ randomly.
2: F ← _H_
3: m ← **0[N]** _[×][k]_

4: Find X satisfying X _[⊤]SX = I._
5: for t = 0, 1, 2, · · ·, T do
6: **if t mod IF = 0 then**

7: _V ←_ _XV_ _[′]_

8: _F ←_ (1 − _α)F + αF_ (V )

9: _F_ _[′]_ _←_ _X_ _[⊤]FX_

10: **end if**

11: Obtain ∇[R]V _[′]_ for one-step decomposition of −F _[′]_ using equation 4.

12: _m ←_ _βm + η∇[R]V_ _[′]_

13: _V_ _V_ + mt

_[′]_ _←_ _[′]_

14: Normalize all column vectors in V _[′]_ following equation 5.

15: end for
16: V _[∗]_ _←_ _XV_ _[′]_


is implemented in Python with NumPy, and optimized via Numba (Lam et al., 2015). The learning
rate η is set to be 10[−][2]. The effective potential matrix function Ueff(·) in Algorithm 2 is based on
Hartree-Fock theory and provided by PySCF (Sun et al., 2018). We use the standard 6-31G basis
set (Ditchfield et al., 1971) for the computation of all molecules. All experiments are run on a AMD
Ryzen 7 5800H CPU (8 cores, 16 threads, 3.2-4.4GHz) with 16GB memory. For the reproducibility


of the result, in all experiments we initialize V _[′]_ in Algorithm 2 by

matrix.

5.1 SCGLED AS AN INITIAL GUESS METHOD


where I is a k × k identity


A direct application of our proposed method is to provide an initial guess for the traditional SCF
method shown in Algorithm 3. We compare with the following popular initial guess methods that
are widely applied in quantum chemistry and solid-state physics software. All of them highly rely
on domain-specific heuristics.

-  Hcore: Core Hamiltonian method, which obtains the guess orbitals from the diagonalization
of the core Hamiltonian H, ignoring all interelectronic interactions.

-  atom (Van Lenthe et al., 2006): Superposition of atomic HF density matrices.

-  minao: Superposition of atomic densities projected in a minimal basis obtained from the
first contracted functions in the cc-pVTZ or cc-pVTZ-PP basis set.

-  huckel (Karton et al., 2017): Parameter-free H¨uckel guess.

In the experiment, we use PySCF’s integrated implementation for all the initial guess baselines.
PySCF (Python-based Simulations of Chemistry Framework, (Sun et al., 2018)) is an open-source,
highly popular quantum chemistry software in Python, whose computationally critical parts are
implemented and optimized in C to guarantee efficiency.

We evaluate the performance of an initial guess method from two aspects

-  Precision: the distance between the initial guess and the final converged solution of equation 1, measured by the averaged energy error over 160 molecules. (expected to be as small
as possible)

-  Convergence: the ratio that the initial guess can successfully lead to a converged solution
via traditional SCF iteration, measured by the number of molecules that failed to converge
in SCF iteration with a certain initial guess method. (expected to be as few as possible)


-----

T = 200 100

10[0]

T = 500

10 2 T = 1000 80 T = 200

60

10 4 T = 2000

SCGLED 40 SCGLED

10 6 Hcore Hcore

Energy error (Hartree) 10 8 10 3 minaoatomhuckel 10 2 T = 5000T = 1000010 1 #(molecules that are not converged) 20 10 3 minaoatomhuckel T = 50010T = 10002 T = 2000T = 5000T = 1000010 1

Wall Time (s) Wall Time (s)


(a) Time-precision evaluation


(b) Time-convergence evaluation


Figure 2: Performance evaluation of initial guess methods

Initial guess methods are expected to be very rapid compared with main SCF iterations, so we restrict
the range of T from 200 to 10,000, corresponding to 5-250 milliseconds, while traditional methods
take around 30 milliseconds. We also set F ’s update interval IF = 50 in this experiment.

The result of performance evaluation is shown in Figure 2. While traditional initial guess methods are parameter-free heuristics, SCGLED is iterative thus the performance and time cost are directly influenced by the number of iterations T . As a result, SCGLED is shown as a curve in timeperformance space while traditional methods are shown as points. In both Figure 2a and Figure 2b,
SCGLED’s curve lies at the lower-left direction of traditional methods, indicating that our method
achieves better results for both precision and convergence ability. Comparing the result of the best
baseline “minao” and SCGLED at T = 1000 in Figure 2a, we can find that averagely SCGLED
is 25x more precise than minao (0.013379 Hartree vs. 0.343194 Hartree in energy error) while the
wall time are very close (26.3 ms vs. 28.6 ms). Our proposed method not only performs better or
close to traditional methods given the same time, but also provides more flexible options, such as
paying more time to achieve better performance in time-rich tasks, or sacrificing performance to satisfy rigorous time limitation, which is inconceivable for traditional methods due to their heuristics
nature. A more detailed, per-molecule result is shown in Appendix E.


T = 200 100 T = 200

10[0] T = 500 T = 500

T = 1000 80 T = 1000

10 2 T = 2000T = 5000 T = 2000T = 5000

60

T = 10000 T = 10000

10 4

40

10 6

Energy error (Hartree)

20

10 8

#(molecules that are not converged)

0

1 5 10 25 50 100 1 5 10 25 50 100

The update interval of F (IF) The update interval of F (IF)


(a) IF -precision relation


(b) IF -convergence relation


Figure 3: IF -performance relation of SCGLED

For SCGLED, we also study the influence of IF, the parameter to control F ’s update interval, to
the performance, shown in Figure 3. The IF -performance relation is shown as a U-shaped curve,
which indicates that IF should be set to a moderate value, neither too large nor too small. The best
value of IF (the bottom of the U-shaped curve) is relevant to the total number of iterations T : IF
should be smaller for small T so F can be updated for reasonable times in very limited iterations,
and should be moderately larger for large T so that the convergence acceleration technique (such as
the damping technique in equation 6) works better.


-----

5.2 SCGLED AS A FULL SOLVER

While gradient-based methods are widely used in machine learning in ways that do not require
extreme precision, it is usually computationally intractable in other fields that require high precision,
especially scientific computing, due to its low convergence rate. However, a counterintuitive result
is that our proposed gradient-like method can produce highly precise solutions of equation 1 in a
reasonable number of iterations.

160

140

120

100

80

60

40

#(molecules that are not solved) 20 SCGLED

SCGLED (with DIIS)

0

10[3] 10[4] 10[5] 10[6] 10[7]

#(iterations)


Figure 4: Performance of SCGLED for solving molecules independently.


(a) Benzene
(C6H6)


(b) n-Pentane
(C5H12)


(c) Chloroethane
(C2ClH5)


(d) Formonitrile
oxide (HCNO)


Figure 5: Visualization of electron density for some molecules, solved via SCGLED.

The result that our proposed method acts as a full solver (without any traditional SCF iteration)
is shown in Figure 4. The convergence criterion is selected as the default setting in PySCF. That
is, the difference of a molecule’s total energy before and after a single-step SCF iteration should
be less than 10[−][10], which is extremely strict. IF is set as 100 since T is relatively large in this
scenario. It shows that our proposed method successfully solved more than half of the molecules
within 20,000 iterations, and 155 out of 160 molecules within 1 million iterations. If we replace the
simple damping technique in Algorithm 2 with a more powerful DIIS method (Pulay, 1980) for the
last 10% iterations[2], we can successfully solve all 160 molecules with very high precision. Some of
the highly precise solutions are visualized in Figure 5 by VMD (Humphrey et al., 1996), in which
the coordinate and type of atoms are given for each molecule, and the electron densities are solved
by SCGLED, represented as isosurfaces shown as grey surface around the molecule.

6 CONCLUSION

In this work, we solve the approximated Schr¨odinger equations from scratch with gradient-like
eigen-decomposition. This work contributes to both the field of computational physics and machine
learning as follows:

2E.g., for T = 10000, damping is applied for the first 9000 iterations, and we apply DIIS for the last 1000
iterations (10 actual DIIS computations are done since IF = 100).


-----

-  For computational physics, beside the performance advantage, this work shows the possibility to solve the approximated Schr¨odinger equations from a purely optimization-based
aspect, without any heuristics method based on prior quantum mechanism knowledge to
bootstrap the solving stage. In this way the solving of approximated Schr¨odinger equations
can be stripped out from its physical background, studied independently as a mathematical
optimization problem.

-  For machine learning, this work explores a brand new area of “self-consistent” eigenvalue problems, especially approximated Schr¨odinger equations, for stochastic eigendecomposition methods such as Oja’s algorithm and EigenGame, which are previously
regarded as specialized methods for k-PCA. While such methods can properly handle
stochasticity, this work shows that they are also capable of handling self-consistency, which
leads to a potential of application in a broader field of scientific computing.

REFERENCES

Zeyuan Allen-Zhu and Yuanzhi Li. First efficient convergence for streaming k-pca: A global, gapfree, and near-optimal rate. In 2017 IEEE 58th Annual Symposium on Foundations of Computer
Science (FOCS), pp. 487–492, 2017. doi: 10.1109/FOCS.2017.51.

R. Ditchfield, W. J. Hehre, and J. A. Pople. Self-consistent molecular-orbital methods. ix. an extended gaussian-type basis for molecular-orbital studies of organic molecules. The Journal of
[Chemical Physics, 54(2):724–728, 1971. doi: 10.1063/1.1674902. URL https://doi.org/](https://doi.org/10.1063/1.1674902)
[10.1063/1.1674902.](https://doi.org/10.1063/1.1674902)

Charlotte Froese Fischer. General hartree-fock program. Computer Physics Communications,
43(3):355–365, 1987. ISSN 0010-4655. doi: https://doi.org/10.1016/0010-4655(87)
90053-1. URL [https://www.sciencedirect.com/science/article/pii/](https://www.sciencedirect.com/science/article/pii/0010465587900531)
[0010465587900531.](https://www.sciencedirect.com/science/article/pii/0010465587900531)

Ian M. Gemp, Brian McWilliams, Claire Vernade, and Thore Graepel. Eigengame: PCA as a nash
equilibrium. In 9th International Conference on Learning Representations, ICLR 2021, Virtual
[Event, Austria, May 3-7, 2021. OpenReview.net, 2021. URL https://openreview.net/](https://openreview.net/forum?id=NzTU59SYbNq)
[forum?id=NzTU59SYbNq.](https://openreview.net/forum?id=NzTU59SYbNq)

D. R. Hartree. The wave mechanics of an atom with a non-coulomb central field. part i. theory
and methods. Mathematical Proceedings of the Cambridge Philosophical Society, 24(1):89–110,
1928. doi: 10.1017/S0305004100011919.

Douglas Rayner Hartree and W. Hartree. Self-consistent field, with exchange, for beryllium. Proceedings of the Royal Society of London. Series A - Mathematical and

Physical Sciences, 150(869):9–33, 1935. doi: 10.1098/rspa.1935.0085. [URL https://](https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1935.0085)
[royalsocietypublishing.org/doi/abs/10.1098/rspa.1935.0085.](https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1935.0085)

Jan Hermann, Zeno Sch¨atzle, and Frank No´e. Deep-neural-network solution of the electronic
Schr¨odinger equation. Nature Chemistry, 12(10):891–897, October 2020. ISSN 1755[4349. doi: 10.1038/s41557-020-0544-y. URL https://www.nature.com/articles/](https://www.nature.com/articles/s41557-020-0544-y)
[s41557-020-0544-y. Number: 10 Publisher: Nature Publishing Group.](https://www.nature.com/articles/s41557-020-0544-y)

John Hertz, Anders Krogh, and Richard G Palmer. Introduction to the theory of neural computation.
CRC Press, 2018.

P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136:B864–B871, Nov
[1964. doi: 10.1103/PhysRev.136.B864. URL https://link.aps.org/doi/10.1103/](https://link.aps.org/doi/10.1103/PhysRev.136.B864)
[PhysRev.136.B864.](https://link.aps.org/doi/10.1103/PhysRev.136.B864)

William Humphrey, Andrew Dalke, and Klaus Schulten. VMD – Visual Molecular Dynamics.

Journal of Molecular Graphics, 14:33–38, 1996.

Amir Karton, Nitai Sylvetsky, and Jan M. L. Martin. W4-17: A diverse and high-confidence
dataset of atomization energies for benchmarking high-level electronic structure methods. Journal
of Computational Chemistry, 38(24):2063–2075, 2017. doi: https://doi.org/10.1002/jcc.24854.
[URL https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.24854.](https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.24854)


-----

M. F. Kasim and S. M. Vinko. Learning the exchange-correlation functional from nature
with fully differentiable density functional theory. Phys. Rev. Lett., 127:126403, Sep 2021.
[doi: 10.1103/PhysRevLett.127.126403. URL https://link.aps.org/doi/10.1103/](https://link.aps.org/doi/10.1103/PhysRevLett.127.126403)
[PhysRevLett.127.126403.](https://link.aps.org/doi/10.1103/PhysRevLett.127.126403)

W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects.

[Phys. Rev., 140:A1133–A1138, Nov 1965. doi: 10.1103/PhysRev.140.A1133. URL https:](https://link.aps.org/doi/10.1103/PhysRev.140.A1133)
[//link.aps.org/doi/10.1103/PhysRev.140.A1133.](https://link.aps.org/doi/10.1103/PhysRev.140.A1133)

Siu Kwan Lam, Antoine Pitrou, and Stanley Seibert. Numba: A llvm-based python jit compiler. In

Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15,
New York, NY, USA, 2015. Association for Computing Machinery. ISBN 9781450340052. doi:
[10.1145/2833157.2833162. URL https://doi.org/10.1145/2833157.2833162.](https://doi.org/10.1145/2833157.2833162)

Susi Lehtola. Assessment of Initial Guesses for Self-Consistent Field Calculations. Superposition
of Atomic Potentials: Simple yet Efficient. Journal of Chemical Theory and Computation, 15
[(3):1593–1604, March 2019. ISSN 1549-9618. doi: 10.1021/acs.jctc.8b01089. URL https:](https://doi.org/10.1021/acs.jctc.8b01089)
[//doi.org/10.1021/acs.jctc.8b01089. Publisher: American Chemical Society.](https://doi.org/10.1021/acs.jctc.8b01089)

Li Li, Stephan Hoyer, Ryan Pederson, Ruoxi Sun, Ekin D. Cubuk, Patrick Riley, and Kieron Burke.
Kohn-sham equations as regularizer: Building prior knowledge into machine-learned physics.
[Phys. Rev. Lett., 126:036401, Jan 2021. doi: 10.1103/PhysRevLett.126.036401. URL https:](https://link.aps.org/doi/10.1103/PhysRevLett.126.036401)
[//link.aps.org/doi/10.1103/PhysRevLett.126.036401.](https://link.aps.org/doi/10.1103/PhysRevLett.126.036401)

Richard M. Martin. Electronic Structure: Basic Theory and Practical Methods. Cambridge University Press, 2004. doi: 10.1017/CBO9780511805769.

Erkki Oja. Simplified neuron model as a principal component analyzer. Journal of mathematical

biology, 15(3):267–273, 1982.

Erkki Oja and Juha Karhunen. On stochastic approximation of the eigenvectors and eigenvalues of
the expectation of a random matrix. Journal of Mathematical Analysis and Applications, 106(1):
[69–84, February 1985. ISSN 0022-247X. doi: 10.1016/0022-247X(85)90131-3. URL https:](https://www.sciencedirect.com/science/article/pii/0022247X85901313)
[//www.sciencedirect.com/science/article/pii/0022247X85901313.](https://www.sciencedirect.com/science/article/pii/0022247X85901313)

David Pfau, James S. Spencer, Alexander G. D. G. Matthews, and W. M. C. Foulkes. Ab initio solution of the many-electron schr¨odinger equation with deep neural networks. Phys. Rev.

Research, 2:033429, Sep 2020. doi: 10.1103/PhysRevResearch.2.033429. [URL https:](https://link.aps.org/doi/10.1103/PhysRevResearch.2.033429)
[//link.aps.org/doi/10.1103/PhysRevResearch.2.033429.](https://link.aps.org/doi/10.1103/PhysRevResearch.2.033429)

P´eter Pulay. Convergence acceleration of iterative sequences. the case of scf iteration.

Chemical Physics Letters, 73(2):393–398, 1980. ISSN 0009-2614. doi: https://doi.org/
[10.1016/0009-2614(80)80396-4. URL https://www.sciencedirect.com/science/](https://www.sciencedirect.com/science/article/pii/0009261480803964)
[article/pii/0009261480803964.](https://www.sciencedirect.com/science/article/pii/0009261480803964)

Terence D. Sanger. Optimal unsupervised learning in a single-layer linear feedforward neural network. Neural Networks, 2(6):459–473, 1989. ISSN 0893-6080. doi: https://doi.org/
[10.1016/0893-6080(89)90044-0. URL https://www.sciencedirect.com/science/](https://www.sciencedirect.com/science/article/pii/0893608089900440)
[article/pii/0893608089900440.](https://www.sciencedirect.com/science/article/pii/0893608089900440)

Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li,
Junzi Liu, James D. McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and
Garnet Kin-Lic Chan. Pyscf: the python-based simulations of chemistry framework. WIREs
Computational Molecular Science, 8(1):e1340, 2018. doi: https://doi.org/10.1002/wcms.1340.
[URL https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1340.](https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1340)

Attila Szabo and Neil S Ostlund. Modern quantum chemistry: introduction to advanced electronic

structure theory. Dover Publications, Inc., 1996.

Cheng Tang. Exponentially convergent stochastic k-pca without variance reduction. In
H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alch´e-Buc, E. Fox, and R. Garnett
(eds.), Advances in Neural Information Processing Systems, volume 32. Curran Asso[ciates, Inc., 2019. URL https://proceedings.neurips.cc/paper/2019/file/](https://proceedings.neurips.cc/paper/2019/file/38faae069a1371784081ea9ad9b279d0-Paper.pdf)
[38faae069a1371784081ea9ad9b279d0-Paper.pdf.](https://proceedings.neurips.cc/paper/2019/file/38faae069a1371784081ea9ad9b279d0-Paper.pdf)


-----

JH Van Lenthe, R Zwaans, Huub JJ Van Dam, and MF Guest. Starting scf calculations by superposition of atomic densities. Journal of computational chemistry, 27(8):926–932, 2006.

Chao Yang, Juan C. Meza, and Lin-Wang Wang. A constrained optimization algorithm for total energy minimization in electronic structure calculations. Journal of Computational Physics, 217(2):
[709–721, 2006. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2006.01.030. URL https:](https://www.sciencedirect.com/science/article/pii/S0021999106000325)
[//www.sciencedirect.com/science/article/pii/S0021999106000325.](https://www.sciencedirect.com/science/article/pii/S0021999106000325)

A PHYSICAL BACKGROUND

In this section, we briefly introduce the physical background of the problem in section 3. We refer
to (Szabo & Ostlund, 1996; Martin, 2004) for more details.

In computational physics, a foundation problem is to solve the Schr¨odinger equation of a many-body
quantum system
_H |Ψ⟩_ = E |Ψ⟩

where H is the Hamilton operator for a system of M nuclei and N electrons described by their
coordinates RA and ri. With the Born-Oppenheimer approximation (nuclei are much heavier than
electrons), we consider the electrons to be moving in the field of fixed nuclei, thus the kinetic energy
of the nuclei (approximated as zero) and the repulsion between the nuclei (approximated to be
constant) can be neglected. In this case we write H as


_N_ _M_

_ZA_

_ri_ _RA_

Xi=1 _AX=1_ _|_ _−_ _|_

coulomb attraction be
tween electrons and nu
{z }

clei


1

_−_ 2 _[∇]i[2]_

_i=1_

X

kinetic energy of the

electrons

| {z }


_H =_


_ri_ _rj_

Xi=1 _j=Xi+1_ _|_ _−_ _|_

repulsion between elec
trons

{z }


and Ψ(r1, _, rN_ ) the electronic wave function, which should be normalized (i.e., Ψ Ψ =
_· · ·_ _⟨_ _|_ _⟩_
ΨM[∗](r1, · · ·ZA, rN )Ψ(r1, · · ·, rN )dr1 · · · drN = 1). For convenience, let operatorN 1 _h(i) = −_ [1]2 _[∇]i[2]_ _[−]_
R _A=1_ _|ri−RA|_ [so][ H][ can be rewritten as][ H][ =][ P]i[N]=1 _[h][(][i][) +][ P]i[N]=1_ _j=i+1_ _|ri−rj_ _|_ [. According to]

the variational principle, the problem of finding solution Ψ for the ground state energy E0 can be

P P

transformed to the following constrained optimization problem

min ⟨Ψ| H |Ψ⟩ _s.t. ⟨Ψ|Ψ⟩_ = 1

Since the multi-electron wave function Ψ(r1, _, rN_ ) is computationally intractable when
_· · ·_
_N is large, a basic way is to approximate it as the product of N orbital wave function_
_ψ1(r1)ψ2(r2) · · · ψN_ (rN ) satisfying ⟨ψi|ψi⟩ = _ψi[∗][(][r][)][ψ][i][(][r][)][dr][ = 1][,][ ∀][i][ ∈]_ [1][ · · ·][ N] [, which is called]
“Hartree approximation”. We also expand ψi(rR) = _u=1_ _[C][ui][φ][(][r][)][ as a linear combination of][ K]_
“basis functions” which is fixed and given, in this case our task becomes to determine the value of
all Cui so as to determine the approximated wave function. To construct a Lagrange multiplier

[P][K]


_L = ⟨Ψ| H |Ψ⟩−_


_ϵi(_ _ψi_ _ψi_ 1)
_⟨_ _|_ _⟩−_
_i=1_

X


we have

_⟨Ψ| H |Ψ⟩_ =


Ψ[∗](r1, _, rN_ )h(i)Ψ(r1, _, rN_ )dr1 _drN_
_· · ·_ _· · ·_ _· · ·_


_i=1_

+
Z

_N_

_i=1_

X


Ψ[∗](r1, _, rN_ )
_· · ·_


_ri_ _rj_
_|_ _−_ _|_ [Ψ(][r][1][,][ · · ·][, r][N] [)][dr][1][ · · ·][ dr][N]


_i=1_ _j=i+1_


_ψi[∗][(][r][)][h][(][i][)][ψ][i][(][r][)][dr][ + 1]_

2


_ψi[∗][(][r][i][)][ψ]j[∗][(][r][j][)]_


_ri_ _rj_
_|_ _−_ _|_ _[ψ][i][(][r][i][)][ψ][j][(][r][j][)][dr][i][dr][j]_


_i=1_


_j≠_ _i_


-----

_CuiCvi_ _φ[∗]u[(][r][)][h][(][i][)][φ][v][(][r][)][dr]_

_i=1_ _u,v_ Z

X X

_N_ _N_

+ [1] _CuiCλjCσjCvi_

2

Xi=1 Xj≠ _i_ _u,v,λ,σX_ Z


_φ[∗]u[(][r][i][)][φ][∗]λ[(][r][j][)]_


_ri_ _rj_
_|_ _−_ _|_ _[φ][σ][(][r][j][)][φ][v][(][r][i][)][dr][i][dr][j]_


and


_ϵi(_ _ψi[∗][(][r][)][ψ][i][(][r][)][dr][ −]_ [1)]
_i=1_ Z

X


_ϵi(_ _ψi_ _ψi_ 1) =
_⟨_ _|_ _⟩−_
_i=1_

X


_φ[∗]u[(][r][)][φ][v][(][r][)][dr][ −]_ [1)]


_ϵi(_ _CuiCvi_
_i=1_ _u,v_

X X


Let _Huv_ = _φ[∗]u[(][r][)][h][(][i][)][φ][v][(][r][)][dr][,]_ _Suv_ = _φ[∗]u[(][r][)][φ][v][(][r][)][dr]_ and _Euvλσ_ =
_φ[∗]u[(][r][1][)][φ][∗]λ[(][r][2][)]_ _r1_ 1 _r2_

_|_ _−_ R| _[φ][σ][(][r][2][)][φ][v][(][r][1][)][dr][1][dr][2][ which are all constants since]R_ _[ {][φ][i][}][ are given functions,]_
we have

R


_CuiCviHuv + [1]_

2

_u,v_

X


_L =_


_CuiCλjCσjCviEuvλσ_

_u,v,λ,σ_

X


_i=1_


_i=1_


_j≠_ _i_


_ϵi(_
_i=1_

X


_CuiCviSuvdr_ 1)
_−_
_u,v_

X


_∂L_

=2
_∂Cui_


_CλjCσjEuvλσ_ _ϵiSuv)_
_λ,σ_ _−_

X


_Cvi(Huv +_


_j(≠_ _i)_


Let _∂C∂Lui_ [= 0][,][ ∀][i][ = 1][ · · ·][ N, u][ = 1][ · · ·][ K][ and we have]


_Cvi(Huv +_ _CλjCσjEuvλσ) = ϵi_ _CviSuv, ∀i = 1 · · · N, u = 1 · · · K_

Xv _jX(≠_ _i)_ Xλ,σ Xv

Let Pλσ = _j(=i)_ _[C][λj][C][σj][,][ [][U][eff][(][P]_ [)]][uv][ =][ P]λ,σ _[P][λσ][E][uvλσ][,][ F][uv][ =][ H][uv][ + [][U][eff][(][P]_ [)]][uv][ and][ Λ =]

_̸_

diag(ϵ1, · · ·, ϵN ), P, Ueff(P ), F ∈ R[K][×][K], we have the matrix form of the above equation

[P] _FC = SCΛ_

which is a generalized eigenvalue problem where the matrix F to be decomposed is defined by the
eigenvectors C.

Actually, there are several improved approximation theories based on the Hartree approximation
mentioned above, in which the definitions of P and Ueff(P ) are different. A most influential one
is the Hartree-Fock theory in which the multi-electron wave function is approximated by a slater
determinant

_χ1(x1)_ _χ2(x1)_ _χN_ (x1)

_· · ·_

Ψ(x1, _, xN_ ) = 1 _χ1(.x2)_ _χ2(.x2)_ _· · ·_ _χN_ (.x2)
_· · ·_ _√N_ ! .. .. ... ..

_χ1(xN_ ) _χ2(xN_ ) _χN_ (xN )

_· · ·_

to conform to the antisymmetry principle Ψ( _, xi,_ _, xj,_ ) = Ψ( _, xj,_ _, xi,_ )

_· · ·_ _· · ·_ _· · ·_ _−_ _· · ·_ _· · ·_ _· · ·_
whose famous representation is Pauli exclusion principle. In this way Puv = 2 _i_ _CλjCσj (or in_
matrix form, P = 2CC _[⊤]) and [Ueff(P_ )]uv = _λ,σ_ _[P][λσ][E][uvλσ][ −]_ 2[1] _λ,σ_ _[P][λσ][E][uλσv][, in which an]_

“exchange term” − [1]2 _λ,σ_ _[P][λσ][E][uλσv][ is added.]_ P [P][N/][2]

[P]

P

B SELF-CONSISTENT FIELD (SCF) METHOD

Self-Consistent Field (SCF) is a standard method to solve equation 1. An initial density matrix
_P0 is generated via heuristics based on prior quantum mechanism knowledge, then the generalized_


-----

eigen-decomposition problem F (Pt 1)Vt = SVtΛt is repeatedly solved to obtain eigenvectors Vt
_−_
and density matrixconvergence threshold. The detailed routine is shown in Algorithm 3. Pt = 2VtVt[⊤] at the t-th iterations, t = 1, 2, · · ·, until |Pt − _Pt−1| is less than the_

To solve the generalized eigen-decomposition problem, it should be transformed to standard form
first. To achieve this, the orthogonalization technique introduced in (Szabo & Ostlund, 1996) is
applied to eliminate the overlap matrix S. First, find a linear transformation X so that X _[⊤]SX = I._
There are several ways to achieve this, a popular one is named “canonical orthogonalization”
which lets X = U diag(s[−][1][/][2]), where U and s are all eigenvectors and eigenvalues of S. Note
that all eigenvalues of S are positive so there is no difficulty of taking square roots. Then let
_V_ _[′]_ = X _[−][1]V (V = XV_ _[′]), we have F_ (V )XV _[′]_ = SXV _[′]Λ. Multiply X_ _[⊤]_ on the left and we
have X _[⊤]F_ (V )XV _[′]_ = X _[⊤]SXV_ _[′]Λ. Let F_ _[′](V ) = X_ _[⊤]F_ (V )X, we have

_F_ _[′](V )V_ _[′]_ = V _[′]Λ_ (9)

which is a standard eigen-decomposition problem.

**Algorithm 3 Self-Consistent Field (SCF) method**

**Input: H, S, Ueff(·) in equation 1 and equation 2. Converge threshold ϵ.**
**Output: V** _[∗], the solution of equation 1_

Obtain initial density matrix P0 via initial guess methods (e.g., SAD or core Hamiltonian).
Find X satisfying X _[⊤]SX = I._
_t ←_ 0
**while** _Pt_ _Pt_ 1 _< ϵ do_

_tFF ←tt |[′] ←[←]t − + 1H[X] +[⊤][F] U−[t][X]eff|_ (Pt)
Obtain precise eigenvectors Vt[′] [of][ F][ ′]t−1 [corresponding to the top-][k][ smallest eigenvalues via]
classical eigen-decomposition methods such as QR iteration.

_VPtt ← ←_ _XV2VttV[′]_ _t[⊤]_

**end while**
_V_ _Vt_

_[∗]_ _←_

C PROOF

Proof of Proposition 1:

_Proof. V_ _[∗]_ being the solution of equation 1 means that V _[∗]_ contains the generalized eigenvectors
corresponding to the top-k smallest generalized eigenvalues of [F (V _[∗]), S]. In Algorithm 1, we_
denote the converged V _[′]_ in the final step as V _[∗′]. With the orthogonalization technique introduced in_
Appendix B, V _[∗]_ being the solution of equation 1 is equivalent to

1. V _[∗′]_ = X _[−][1]V_ _[∗]_ contains the eigenvectors corresponding to the top-k smallest eigenvalues
of F _[∗′]._

2. F _[∗′]_ = X _[⊤]F_ (V _[∗])X._

When V _[′]_ in the iteration of Algorithm 1 is equal to the converged solution V _[∗′], F_ _[∗′]_ = X _[⊤]F_ (V _[∗])X_
holds. Since V _[∗′]_ is already converged, the one-step decomposition of −F _[∗′]_ will produce no update for V _[∗′]. For Oja’s algorithm equation 3, letting M = −F_ _[∗′], this means (v1[∗][,][ · · ·][, v]k[∗][) =]_
_QR(v1[∗]_ [+][ ηMv]1[∗][,][ · · ·][, v]k[∗] [+][ ηMv]k[∗][)][, η >][ 0][. Note that][ v]1[∗][,][ · · ·][, v]k[∗] [are orthonormal guaranteed by the]
QR decomposition. Considering that the QR decomposition is implemented via the Gram-Schmidt
process, starting from i = 1, we have


_β1v1[∗]_ [=][ v]1[∗] [+][ ηMv]1[∗] 1 [=][ β][1][ −] [1] _v1[∗]_

_[⇒]_ _[Mv][∗]_ _η_

_β2v2[∗]_ [=][ v]2[∗] [+][ ηMv]2[∗] _[−]_ [[(][v]2[∗] [+][ ηMv]2[∗][)][⊤][v]1[∗][]][v]1[∗]


-----

= v2[∗] [+][ ηMv]2[∗] _[−]_ [[][v]2[∗⊤][v]1[∗] [+][ ηv]2[∗⊤][M][ ⊤][v]1[∗][]][v]1[∗]

_β1_ 1
= v2[∗] [+][ ηMv]2[∗] 2 _−_ _v1[∗][]][v]1[∗]_

_[−]_ [[][ηv][∗⊤] _η_

= v2[∗] [+][ ηMv]2[∗] 2 [=][ β][2][ −] [1] _v2[∗]_

_[⇒]_ _[Mv][∗]_ _η_

_· · ·_

factors. Therefore,in which β1 = ∥v1[∗] v[+]1[∗][ ηMv][,][ · · ·][, v]1[∗][∥]k[∗][, β][are all eigenvectors of][2] [=][ ∥][v]2[∗] [+][ ηMv]2[∗] _[−][ −][[(][v][F]2[∗][ ∗′][+][.]_ _[ ηMv]2[∗][)][⊤][v]1[∗][]][v]1[∗][∥][,][ · · ·][ are normalization]_

To show that v1[∗][,][ · · ·][, v]k[∗] [are eigenvectors corresponding to the top-][k][ eigenvalues, we mainly follow]
(Oja & Karhunen, 1985) and (Hertz et al., 2018). We start from k = 1 to show that v1[∗] [should]
correspond to the largest eigenvalue λ1 of −F _[∗′]_ to be a stable convergence point.

First, for v1, the QR decomposition degenerates to normalization, that is

_F_ _[′]_ _←_ _X_ _[⊤]F_ (v1)X, _v1[′]_ _[←]_ _[v][1]_ [+][ η][(][−][F][ ′][)][v][1][,] _v1 ←_ _∥vv11[′][′]_ _[∥]_

Assuming η is small enough, it can be expanded as a power series of η. By ignoring O(η[2]) terms,
we have
_F_ _[′]_ _←_ _X_ _[⊤]F_ (v1)X, _v1 ←_ _v1 + η∆v1,_ ∆v1 = (−F _[′])v1 −_ [v1[⊤][(][−][F][ ′][)][v][1][]][v][1]
While we have a converged solution v1[∗][, we already know that it is an eigenvector but do not know]
which eigenvalue it corresponds to. In this case, we assume that after a series of iterations, v1 is
very close to the eigenvector corresponding to the αth largest eigenvalue (denoted as v[α]) of the final
converged matrix −F _[∗′]_ = −X _[⊤]F_ (v1[∗][)][X][, that is]
_v1 = v[α]_ + ϵ
where ϵ is a very small perturbation vector. Since we already know that v1[∗] [is an eigenvector, we]
have v1[∗] [=][ v][α][.]

Let M = −F _[∗′]_ = −X _[⊤]F_ (v1[∗][)][X][ and][ E][ =][ X] _[⊤][F]_ [(][v]1[∗][)][X][ −] _[X]_ _[⊤][F]_ [(][v][1][)][X][ =][ −][F][ ′][ −] _[M][ ⇒−][F][ ′][ =]_
_M + E. From the assumption in the proposition, we have Ev1 = O(ϵ[2]). Then we do one-step_
iteration as follows:
∆v1 = (−F _[′])v1 −_ [v1[⊤][(][−][F][ ′][)][v][1][]][v][1]
= (M + E)v1 [v1[⊤][(][M][ +][ E][)][v][1][]][v][1]
_−_
= M (v[α] + ϵ) − [(v[α] + ϵ)[⊤]M (v[α] + ϵ)](v[α] + ϵ) + O(ϵ[2])

= (M + E)(v[α] + ϵ) − [(v[α] + ϵ)[⊤](M + E)(v[α] + ϵ)](v[α] + ϵ)

= λ[α]v[α] + Mϵ − (v[α][⊤]Mv[α])v[α] _−_ (ϵ[⊤]Mv[α])v[α] _−_ (v[α]Mϵ)v[α] _−_ (v[α][⊤]Mv[α])ϵ + O(ϵ[2])

= Mϵ − 2λ[α](ϵ[⊤]v[α])v[α] _−_ _λ[α]ϵ + O(ϵ[2])_

Now we want to analyse the direction of ∆v1 to see whether it can conquer the perturbation ϵ and
converge to v1[∗] [=][ v][α][ stably. We do it by multiplying another eigenvector][ v][β][⊤] [on the left (ignoring]
the O(ϵ[2]) term)

_v[β][⊤]∆v1 = (M_ _[⊤]v[β])[⊤]ϵ −_ 2λ[α](ϵ[⊤]v[α])δαβ − _λ[α]v[β][⊤]ϵ_

= (λ[β] _λ[α])v[β][⊤]ϵ_ 2λ[α](ϵ[⊤]v[α])δαβ
_−_ _−_
We notice that, if β > α (i.e., λ[β] _> λ[α]), there will always exist a direction ϵ so that both v[β][⊤]ϵ and_
_v[β][⊤]∆v1 are larger than zero. In this case v1 will flip to the other eigenvector thus cannot converge_
to v[α] stably. Therefore, α should be 1 so that v1[∗] [=][ v][α][ corresponds to the largest eigenvalue of]
_−F_ _[∗′]._

For top-k case, take k = 2 as an example, notice that the Gram-Schmidt process not only normalize
the eigenvector v2 but also orthogonalize it towards v1, so (v2 + ∆v2)[⊤]v1 = (∆v2)[⊤]v1 = 0, and
we already know that v1 will converge to the eigenvector corresponding to the largest eigenvalue.
Therefore we can only select other eigenvectors out of v[1] to analyse the direction of ∆v2, which
shows that v2 will converge to the eigenvector corresponding to the second largest eigenvalue.


-----

D FULL-SOLVER PERFORMANCE BENCHMARK

While SCGLED is not targeted to be an extremely precise (convergence threshold less than 10[−][10]),
end-to-end solving technique due to its first-order nature, we still compare its performance toward
traditional SCF methods in a full-solver setting for reference.

We test three full-solving schemes

-  Independent SCGLED: applying our proposed method, SCGLED, as a full solver without
any tranditional SCF methods.

-  SCF + minao: applying the traditional SCF method with “minao” as the initial guess
method, which performs best within traditional methods as shown in subsection 5.1, and
acts as the default setting in PySCF.

-  SCF + SCGLED: applying the traditional SCF method with SCGLED as the initial guess
method.

with three settings for heuristics convergence acceleration techniques

-  Vanilla: no convergence acceleration method is applied.

-  Damping: the damping technique shown in equation 6 is applied in SCGLED or SCF
iteration.

-  DIIS: the DIIS technique (Pulay, 1980) is applied in independent SCGLED or SCF iteration. Note that DIIS is not applied when SCGLED is served as an initial guess method.

So excluding vanilla SCGLED in Algorithm 1 whose empirical performance is too poor for benchmarking, we have 8 test methods in total to be benchmarked. In Table 1, we evaluate these methods
from two aspects

-  Convergence: evaluated via the number of successfully solved molecules within the total
160 singlet molecules in W4-17 dataset.

-  Efficiency: evaluated via the averaged[3] end-to-end time cost of the method along the
molecules that are successfully solved.

For the parameter setting, when SCGLED is served as an initial method, we set T = 5000 and
_IF = 50 for Vanilla and Damping cases, which is relatively large since their performance are more_
sensitive to the quality of initial guess. For the DIIS case which are relatively not so sensitive to
the initial guess, we set T = 200 and IF = 25. When SCGLED is running independently without
traditional SCF iterations, the parameter setting stays the same as in subsection 5.2.

|Col1|Vanilla Damping DIIS|
|---|---|
|Independent SCGLED SCF + minao SCF + SCGLED|— 155 / 574.3 ms 160 / 657.9 ms 141 / 332.6 ms 159 / 244.3 ms 160 / 160.9 ms 149 / 190.3 ms 159 / 207.1 ms 160 / 148.7 ms|



Table 1: Benchmark results for three full-solving schemes with three heuristics convergence acceleration settings. In each cell, the former is the number of successfully solved molecules within 160
test molecules, and the latter is the averaged time cost.

The benchmark result is shown in Table 1 and Figure 6. While independent SCGLED works in the
full-solver setting, there is still a performance gap with other methods that target for full-solving
performance. Instead, SCF + SCGLED works best, especially in the Vanilla case that both the
convergence performance (unsolved molecules from 19 to 11) and efficiency performance (averaged time cost from 332.6 ms to 190.3 ms) are significantly improved. For the Damping and DIIS
case, the performance advantages are narrowed since the heuristics convergence acceleration techniques themselves become a more dominant factor for end-to-end performance. Also note that, as an

3Note that the time cost grows exponentially towards the size of the molecule, so the distribution of time
cost is quite skewed along different size of molecules. For this reason we take the logarithm of the time cost
before the average.


-----

|160|SCGLED (Damping) SCGLED (DIIS) SCF + minao (Vanilla) SCF + minao (Damping) SCF + minao (DIIS) SCF + SCGLED (Vanilla) SCF + SCGLED (Damping) SCF + SCGLED (DIIS)|
|---|---|
|||


10 1 10[0] 10[1] 10[2]

SCGLED (Damping)
SCGLED (DIIS)
SCF + minao (Vanilla)
SCF + minao (Damping)
SCF + minao (DIIS)
SCF + SCGLED (Vanilla)
SCF + SCGLED (Damping)
SCF + SCGLED (DIIS)

1 [0] [1] [2]

time cost (s)

Figure 6: Full-solver Performance Benchmark.

optimization-based method, SCGLED brings more randomness compared with traditional heuristics
initial guess methods that are built on domain knowledge. Therefore, while SCF + SCGLED performs better for most of the molecules, its distribution of time costs is wider, shown as curves with
smaller slope in Figure 6, and the end of the curve in Damping and DIIS settings can correspond to
larger time cost compared with SCF + minao.

E RESULT PER MOLECULE

In this section, we provide the description of the full numerical result listed in the supplementary
material of this paper, for all the 160 singlet molecules in the W4-17 dataset. For each molecule and
method, we list four values as follows:

-  Energy: the total Hartree-Fock energy of the molecule computed with the method (Unit:
Hartree, 1 Hartree = 4.3598 × 10[−][10] Joule).

-  Energy error: the difference between the computed energy above and the exact HartreeFock energy listed in the second column of the table (expected to be as small as possible).
Note that we only keep 6 decimals due to the space limit, while the converge criterion is
stricter, so there will be some cases that the error is 0.000000 while the number of iterations
in the fourth line is not equal to zero.

-  Wall time: the wall time to run the method for the corresponding molecule.

-  Convergence: whether the computed result can reach to a converged solution with the
vanilla SCF iteration shown in Algorithm 3 (cannot reach a converged solution via SCF:
; can reach a converged solution via SCF: ; already a converged solution: ). If it can
reach to a converged solution, we also show the number of SCF iterations it costs. If the
number of SCF iterations is equal to 0, it means that the method already reaches a precise
(converged) solution before the SCF iteration.

In the result table, we bold the result if it is a dominated result over the baselines (the energy error
and wall time are both smaller than the best of the baseline methods while the number of iterations is
not larger than the best baseline). We underline the result if it takes more time than the best baseline
while the other target values are dominated. The Hcore method are not taken into account here since
both the energy error and the convergence performance are very unsatisfactory, even though it is
extremely fast.

The table shows that

-  75 out of 160 molecules contain at least one dominated result (in bold).

-  155 out of 160 molecules contain at least one more precise result (marked with underline).


-----

-  72 out of 160 molecules contain at least one result that already reached the converged solution (marked with , the number of SCF iterations is equal to 0) within 10,000 iterations
of SCGLED. That is, we already obtained the converged, precise solution. No further SCF
iterations are needed.

-  14 out of 160 molecules are successfully led to convergence by SCGLED while all the four
baseline methods failed to do so.

The full table of numerical result can be found in “full numerical result.pdf” as a supplementary material of this paper.


-----

