cleanup conclusion

This commit is contained in:
Pierre-Francois Loos 2022-09-20 15:09:42 +02:00
parent 5e7edd8b37
commit db7e8089c5

132
g.tex
View File

@ -557,7 +557,7 @@ Note that it is possible only if the dimension of the domains is not too large (
%\label{sec:time}
%--------------------------------------------%
In this section we generalize the path-integral expression of the Green's matrix, Eqs.~\eqref{eq:G} and \eqref{eq:cn_stoch}, to the case where domains are used.
For that we introduce the Green's matrix associated with each domain
To do so, we introduce the Green's matrix associated with each domain as follows:
\be
G^{(N),\cD}_{IJ}= \mel{ J }{ T_I^N }{ I }.
\ee
@ -565,30 +565,30 @@ For that we introduce the Green's matrix associated with each domain
%\be
%G^{(N)}_{i_0 i_N}= \sum_{i_1,...,i_{N-1}} \prod_{k=0}^{N-1} \langle i_k| T |i_{k+1} \rangle.
%\ee
Starting from Eq.~\eqref{eq:cn} and using the representation of the paths in terms of exit states and trapping times, we get
Starting from Eq.~\eqref{eq:cn} and using the representation of the paths in terms of exit states and trapping times [see Eq.~\eqref{eq:eff_series}], we get
\be
G^{(N)}_{I_0 I_N} = \sum_{p=0}^N
\sum_{\cC_p} \sum_{(i_1,...,i_{N-1}) \in \cC_p}
\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} }
\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} } ,
\ee
where $\cC_p$ is the set of paths with $p$ exit states, $\ket{I_k}$, and trapping times $n_k$ with the constraints that $\ket{I_k} \notin \cD_{k-1}$ (with $1 \le n_k \le N+1$ and $\sum_{k=0}^p n_k= N+1$).
We then have
It follows that
\begin{multline}
\label{eq:Gt}
G^{(N)}_{I_0 I_N}= G^{(N),{\cal D}}_{I_0 I_N} +
\sum_{p=1}^{N}
\sum_{|I_1\rangle \notin {\cal D}_{I_0}, \hdots , |I_p\rangle \notin {\cal D}_{I_{p-1}} }
\sum_{|I_1\rangle \notin {\cal D}_{I_0}, \ldots , |I_p\rangle \notin {\cal D}_{I_{p-1}} }
\sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 1}
\delta_{\sum_{k=0}^p n_k,N+1}
\\
\times
\qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} } ]
G^{(n_p-1),{\cal D}}_{I_p I_N}
G^{(n_p-1),{\cal D}}_{I_p I_N},
\end{multline}
where $\delta_{i,j}$ is a Kronecker delta.
This expression is the path-integral representation of the Green's matrix using only the variables $(\ket{I_k},n_k)$ of the effective dynamics defined over the set of domains.
The standard formula derived above [see Eq.~\eqref{eq:G}] may be considered as the particular case where the domain associated with each state is empty,
The standard formula derived above [see Eq.~\eqref{eq:G}] may be considered as the particular case where the domain associated with each state is empty.
In this case, $p=N$ and $n_k=1$ (for $0 \le k \le N$) and we are left only with the $p$th component of the sum, that is, $G^{(N)}_{I_0 I_N}
= \prod_{k=0}^{N-1} \mel{ I_k }{ F_{I_k} }{ I_{k+1} }$, with $F=T$.
@ -616,7 +616,7 @@ and using the effective transition probability defined in Eq.~\eqref{eq:eq3C}, w
\qty( \prod_{k=0}^{p-1} W_{I_k I_{k+1}} )
\bar{G}^{(n_p-1), {\cal D}}_{I_p I_N} },
\ee
where the average of a given function $F$ is defined as
where, in this context, the average of a given function $F$ is defined as
\begin{multline}
\expval{F}
= \sum_{\ket{I_1} \notin \cD_{I_0}, \ldots , \ket{I_p} \notin \cD_{I_{p-1}}}
@ -624,18 +624,18 @@ where the average of a given function $F$ is defined as
\delta_{\sum_k n_k,N+1}
\\
\times
\prod_{k=0}^{N-1}\cP_{I_k \to I_{k+1}}(n_k-1) F(I_0,n_0;\ldots.;I_N,n_N).
\prod_{k=0}^{N-1}\cP_{I_k \to I_{k+1}}(n_k-1) F(I_0,n_0;\ldots;I_N,n_N).
\end{multline}
In practice, a schematic DMC algorithm to compute the average is as follows.\\
i) Choose some initial vector $\ket{I_0}$\\
\titou{i) Choose some initial vector $\ket{I_0}$\\
ii) Generate a stochastic path by running over $k$ (starting at $k=0$) as follows.\\
$\;\;\;\bullet$ Draw $n_k$ using the probability $P_{I_k}(n)$ [see Eq.~\eqref{eq:PiN}]\\
$\;\;\;\bullet$ Draw the exit state, $\ket{I_{k+1}}$, using the conditional probability $$\frac{\cP_{I_k \to I_{k+1}}(n_k)}{P_{I_k}(n_k)}$$\\
iii) Terminate the path when $\sum_k n_k=N$ is greater than some target value $N_{\rm max}$ and compute $F(I_0,n_0;...;I_N,n_N)$\\
iii) Terminate the path when $\sum_k n_k=N$ is greater than some target value $N_{\rm max}$ and compute $F(I_0,n_0;\ldots;I_N,n_N)$\\
iv) Go to step ii) until some maximum number of paths is reached.\\
\\
At the end of the simulation, an estimate of the average for a few values of $N$ greater but close to $N_\text{max}$ is obtained.
At large $N_\text{max}$ where the convergence of the average as a function of $p$ is reached, such values can be averaged.
At large $N_\text{max}$ where the convergence of the average as a function of $p$ is reached, such values can be averaged.}
%==============================================%
\subsection{Domain Green's Function Monte Carlo}
@ -643,30 +643,26 @@ At large $N_\text{max}$ where the convergence of the average as a function of $p
%==============================================%
The aim of this section is to show that it is possible to go further by integrating out the trapping times, $n_k$, of the preceding expressions, thus defining a new effective stochastic dynamics involving now only the exit states.
Physically, it means that we are going to compute exactly within the time-evolution of all stochastic paths trapped within each domain.
Physically, it means that we are going to compute exactly the time evolution of all stochastic paths trapped in each domain.
We shall present two different ways to derive the new dynamics and renormalized probabilistic averages.
The first one, called the pedestrian way, consists in starting from the preceding time-expression for $G$ and make the explicit integration over the $n_k$'s.
The first one, called the pedestrian way, consists in starting from the preceding time expression for $G$ and make the explicit integration over the $n_k$'s.
The second, more direct and elegant, is based on the Dyson equation.
%--------------------------------------------%
\subsubsection{The pedestrian way}
%--------------------------------------------%
Let us define the quantity\\
$$
G^E_{ij}= \tau \sum_{N=0}^{\infty} \mel{ i }{ T^N }{ j}.
$$
By summing over $N$ we obtain
Let us define the energy-dependent Green's matrix
\be
G^E_{ij}= \mel{i}{ \qty( H-E \Id )^{-1} }{j}.
G^E_{ij}= \tau \sum_{N=0}^{\infty} \mel{ i }{ T^N }{ j} = \mel{i}{ \qty( H-E \Id )^{-1} }{j},
\ee
This quantity, which no longer depends on the time step, is referred to as the energy-dependent Green's matrix.
which does not depends on the time step.
Note that, \titou{in a continuous space}, this quantity is essentially the \titou{Laplace transform of the time-dependent Green's function}.
Here, we then use the same denomination.
We shall use the same denomination in the following.
The remarkable property is that, thanks to the summation over $N$ up to the infinity the constrained multiple sums appearing in Eq.~\eqref{eq:Gt} can be factorized in terms of a product of unconstrained single sums as follows
\be
\sum_{N=1}^\infty \sum_{p=1}^N \sum_{n_0 \ge 1} ...\sum_{n_p \ge 1} \delta_{n_0+...+n_p,N+1} \titou{??}
= \sum_{p=1}^{\infty} \sum_{n_0=1}^{\infty} ...\sum_{n_p=1}^{\infty} \titou{??}.
\sum_{N=1}^\infty \sum_{p=1}^N \sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 1} \delta_{n_0+...+n_p,N+1} \titou{??}
= \sum_{p=1}^{\infty} \sum_{n_0=1}^{\infty} \cdots \sum_{n_p=1}^{\infty} \titou{??}.
\ee
It is then a trivial matter to integrate out exactly the $n_k$ variables, leading to
\begin{multline}
@ -691,9 +687,9 @@ In fact, there is a more direct way to derive the same equation by resorting to
where $H_0$ is some arbitrary reference Hamiltonian, we have the Dyson equation
\be
\label{eq:GE}
G^E_{ij}= \titou{G^E_{0,ij}} + \sum_{kl} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj}
G^E_{ij}= \titou{G^E_{0,ij}} + \sum_{kl} G^{E}_{0,ik} (H_0-H)_{kl} G^E_{lj}.
\ee
Let us choose as $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$.
Let us choose $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$.
Then, The Dyson equation \eqref{eq:GE} becomes
\begin{multline}
\mel{ i }{ \qty(H-E \Id)^{-1} }{ j }
@ -701,16 +697,16 @@ Then, The Dyson equation \eqref{eq:GE} becomes
\\
+ \sum_k \mel{ i }{ P_i \qty(H-E \Id)^{-1} P_i(H_0-H) }{ k } \mel{ k }{ \qty(H-E \Id)^{-1} }{ j }.
\end{multline}
Now, we have
Using the following identity
\be
\begin{split}
P_i \qty(H-E \Id)^{-1} P_i(H_0-H)
& = P_i \qty(H-E \Id)^{-1} P_i (P_i H P_i - H)
\\
& = P_i \qty(H-E \Id)^{-1} P_i (-H) (1-P_i)
& = P_i \qty(H-E \Id)^{-1} P_i (-H) (1-P_i),
\end{split}
\ee
and the Dyson equation may be written under the form
the Dyson equation may be written under the form
\begin{multline}
\mel{ i }{ \qty(H-E \Id)^{-1} }{ j }
= \mel{ i }{ P_i \qty(H-E \Id)^{-1} P_i }{ j}
@ -718,7 +714,6 @@ and the Dyson equation may be written under the form
+ \sum_{k \notin \cD_i} \mel{ i }{ P_i \qty(H-E \Id)^{-1} P_i (-H)(1-P_i) }{ k } \mel{k}{\qty(H-E \Id)^{-1}}{j}
\end{multline}
which is identical to Eq.~\eqref{eq:eqfond} when $G$ is expanded iteratively.
\\
Let us use as effective transition probability density
\be
P_{I \to J} = \frac{1} {\PsiG(I)} \mel{ I }{ P_I \qty(H^+ - \EL^+ \Id)^{-1} P_I (-H^+) (1-P_I) }{ J } \PsiG(J)
@ -730,7 +725,7 @@ and the weight
{\mel{ I }{ \qty(H^+ - \EL^+ \Id)^{-1} P_I (-H^+)(1-P_I) }{ J} }
\ee
Using Eqs.~\eqref{eq:eq1C}, \eqref{eq:eq3C} and \eqref{eq:relation}, one can easily verify that $P_{I \to J} \ge 0$ and $\sum_J P_{I \to J}=1$.
Finally, the probabilistic expression writes
Finally, the probabilistic expression reads
\begin{multline}
\label{eq:final_E}
\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ I_N }
@ -742,22 +737,27 @@ Finally, the probabilistic expression writes
%----------------------------%
\subsubsection{Energy estimator}
%----------------------------%
To calculate the energy we introduce the following quantity
To calculate the energy we introduce the following energy estimator
\be
\cE(E) = \frac{ \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ H\PsiT } } {\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \PsiT } }.
\cE(E) = \frac{ \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ H\PsiT } } {\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \PsiT } },
\label{calE}
\ee
and search for the solution $E=E_0$ of $\cE(E)= E$.
Using the spectral decomposition of $H$ we have
Using the spectral decomposition of $H$, we have
\be
\label{eq:calE}
\cE(E) = \frac{ \sum_i \frac{E_i c_i}{E_i-E}}{\sum_i \frac{c_i}{E_i-E}}
\ee
with \titou{$c_i = \braket{ I_0 }{ \Phi_0 } \braket{ \Phi_0}{ \PsiT }$}.
with $c_i = \braket{ I_0 }{ \Phi_i } \braket{ \Phi_i}{ \PsiT }$, \titou{where $\Phi_i$ are eigenstates of $H$}.
It is easy to check that, in the vicinity of $E = E_0$, $\cE(E)$ is a linear function of $E - E_0$.
Therefore, in practice, we compute several values of \titou{$\cE(E^{k})$} and fit the data using some function of $E$ close to the linearity to extrapolate the exact value of $E_0$.
%Let us describe the 3 functions used here for the fit.
Therefore, in practice, we compute the value of $\cE(E)$ for several values of $E$, and fit these data using a linear or quadratic function of $E$ in order to obtain, via extrapolation, the exact value $E_0$.
In order to have a precise extrapolation of the energy, it is best to compute $\cE(E)$ for values of $E$ as close as possible to the exact energy.
However, as $E \to E_0$, both the numerators and denominators of Eq.~\eqref{eq:calE} diverge.
This is reflected by the fact that one needs to compute more and more $p$-components with an important increase of statistical fluctuations.
Thus, from a practical point of view, a trade-off has to be found between the \titou{quality} of the extrapolation and the amount of statistical fluctuations.
%Let us describe the 3 functions used here for the fit.
%i) Linear fit\\
%We write
%\be
@ -785,10 +785,6 @@ Therefore, in practice, we compute several values of \titou{$\cE(E^{k})$} and fi
%\ee
%Here, the variational parameters used for the fit are ($c_0,E_0,c_1,E_1)$.\\
In order to have a precise extrapolation of the energy, it is interesting to compute the ratio $\cE(E)$ for values of $E$ as close as possible to the exact energy.
However, in that case the numerators and denominators computed diverge.
This is reflected by the fact that we need to compute more and more $p$-components with an important increase of statistical fluctuations.
So, in practice a trade-off has to be found between the possible bias in the extrapolation and the amount of simulation time required.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical application to the Hubbard model}
@ -1137,28 +1133,27 @@ $N$ & Size Hilbert space & Domain & Domain size & $\alpha,\beta$ &$\bar{t}_{I_0}
\label{sec:conclu}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In this work it has been shown how to integrate out exactly within a DMC framework the contribution of all
stochastic trajectories trapped in some given domains of the Hilbert space and the corresponding general equations have been derived.
In this way a new effective stochastic dynamics connecting only the domains and not the individual states is defined.
A key property of the effective dynamics is that it does not depend on the shape of the domains used for each state, so rather general
domains overlapping or not can be used. To obtain the effective transition probability giving the probability of going from
one domain to another and the renormalized estimators,
the Green's functions limited to the domains are to be computed analytically, that is, in practice, matrices of size the number
of states of the sampled domains need to be inverted. This is the main computationally intensive step of the method.
The efficiency of the method is directly related to
the importance of the average time spent by the stochastic trajectories within each domain. Being able to define domains with large average trapping times
is the key aspect of the method since it may lead to some important reduction of the statistical error as illustrated in our numerical application.
So a tradeoff has to be found in the complexity of the domains maximizing the average trapping time and the cost of computing
the local Green's functions associated with such domains. In practice, there is no general rule to construct efficient domains.
For each system at hand, we need to determine on physical grounds which regions of the space are preferentially sampled by
the stochastic trajectories and to build domains of minimal size enclosing such regions. In the first application presented here on the 1D-Hubbard model
we exploit the physics of the large $U$ regime which is known to approach the Heisenberg limit where double occupations have a small weight.
This simple example has been chosen to illustrate the various aspects of the approach. However, our goal is of course to
tackle much larger systems, like those treated by state-of-the-art methods, such as selected CI, (for example, \cite{Huron_1973,Harrison_1991,Giner_2013,Holmes_2016,
Schriber_2016,Tubman_2020}, QMCFCI,\cite{Booth_2009,Cleland_2010}, AFQMC\cite{Zhang_2003}, or DMRG\cite{White_1999,Chan_2011} approaches).
Here, we have mainly focused on the theoretical aspects of the approach. To reach larger systems will require some
more elaborate implementation of the method in order to keep under control the cost of the simulation.
Doing this was out of the scope of the present work and will be presented in a forthcoming work.
In this work, it has been shown how to integrate out exactly --- within a DMC framework --- the contribution of all stochastic trajectories trapped in some given domains of the Hilbert space.
The corresponding equations have been derived in a general context.
In such a way, a new effective stochastic dynamics connecting only domains (and not the individual states) is defined.
A key property of this effective dynamics is that it does not depend on the ``shape'' of the domains used for each state.
Therefore, rather general domains (with or without overlap) can be considered.
To obtain the effective transition probability (which provides the probability of going from one domain to another) and the corresponding renormalized estimators, the domain Green's functions must be computed analytically, that is, in practice, matrices of the size of the number of states for the sampled domains have to be inverted.
This is the main computationally intensive step of the present approach.
The efficiency of the method is directly related to the importance of the average time spent by the stochastic trajectories in each domain.
Being able to define domains with large average trapping times is the key aspect of the method since it may lead to some important reduction of the statistical error as illustrated in our numerical applications.
Therefore, a trade-off has to be found between maximizing the average trapping time and minimizing the cost of computing the local Green's functions.
In practice, there is no general rule to construct such domains.
For each system at hand, one needs to determine, on physical grounds, which \titou{spatial} regions are preferentially sampled by the stochastic trajectories and to build domains of minimal size enclosing such regions.
In the first application presented here on the 1D Hubbard model we exploit the physics of the large-$U$ regime that is known to approach the Heisenberg limit where double occupations have small weights.
This simple example has been chosen to illustrate the various aspects of the approach.
Our goal is, of course, to tackle much larger systems, like those treated by state-of-the-art methods, such as selected CI, \cite{Huron_1973,Harrison_1991,Giner_2013,Holmes_2016,Schriber_2016,Tubman_2020}, FCIQMC, \cite{Booth_2009,Cleland_2010}, AFQMC, \cite{Zhang_2003} or DMRG. \cite{White_1999,Chan_2011}
Here, we have mainly focused on the theoretical aspects of the approach.
In order to consider larger systems, an elaborate implementation of the present method is necessary in order to keep under control the cost of the simulation.
This is outside the scope of the present study and will be presented in a forthcoming work.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
@ -1180,16 +1175,17 @@ under grant agreement no.~952165.
\label{app:A}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
For the simplest case of a two-state system the fundamental equation (\ref{eqfond}) writes
$$
For the simplest case of a two-state system the fundamental equation \eqref{eq:eqfond} reads
\be
{\cal I}=
\langle I_0|\frac{1}{H-E}|\Psi\rangle = \langle I_0|P_0\frac{1}{H-E} P_0|\Psi\rangle
+ \sum_{p=1}^{\infty} {\cal I}_p$$
+ \sum_{p=1}^{\infty} {\cal I}_p
\ee
with
$$
\be
{\cal I}_p=
\sum_{I_1 \notin {\cal D}_0, \hdots , I_p \notin {\cal D}_{p-1}}
$$
\ee
\be
\Big[ \prod_{k=0}^{p-1} \langle I_k| P_k \frac{1}{H-E} P_k (-H)(1-P_k)|I_{k+1} \rangle \Big]
\langle I_p| P_p \frac{1} {H-E} P_p|\Psi\rangle