What are the "smart algorithms" applied to solve the "curse of dimensionality"?

Question

The "curse of dimensionality" is an ubiquitous issue arising in both electronic structure and quantum molecular dynamics, which refers to the exponential scaling of computational cost with the number of degrees of freedom of the systems of interest.

This problem is manifested in many applications of computational studies (e.g. determining transition states of chemical reactions, geometry optimization of large molecular systems, simulating spectra of molecular systems, etc). The "curse of dimensionality" dramatically restricts the practical size of systems that can be studied computationally.

Assuming that only classical computers are available, what are the "smart algorithms" that have been devised to partially solve this "curse of dimensionality" problem and make the subset of the full problem be close to linear scaling?

Answer 1

This is a very broad question, so I am going to give a very brief overview of typical exponentially-scaling problems. I am not an expert in most of these areas, so any suggestions or improvements will be welcome.

Solving the Schrödinger equation

In order to solve the Schrödinger equation numerically, you need to diagonalise a rank $3N$ tensor -- as you can see, a pretty impossible operation not only in terms of CPU power, but also in terms of memory. The main problem in fact is that the wavefunction has to be antisymmetric with respect to all electrons, which is the main reason for combinatorial explosion. An alternative way is to expand the wavefunction as a multivariable Taylor series of antisymmetric functions (determinants), and if you were to do it exactly (full configuration interaction), it also scales exponentially. So at this point you can either solve the equation by ignoring most of the correlation between different degrees of freedom (Hartree-Fock, Moller-Plesset perturbation theory, truncated configuration interaction), project the $3N$-dimensional problem onto a 3-dimensional one, where the exact solution is unknown, but able to be approximated (density functional theory), or solve the correlated problem exactly for an idealised approximate infinite sum (coupled cluster theory). Another way to solve the equation is to convert it to a sampling problem (diffusion quantum Monte Carlo), which is exact for bosons, but needs an approximation for fermions (fixed node approximation), so that it doesn't scale exponentially. There is a lot of literature on making a lot of the above methods linear-scaling using clever approximations or making the formally exact full configuration interaction method more efficient (full configuration interaction quantum Monte Carlo), but in general, the more computational time you throw in, the larger the class of problems your method can tackle and some of the above approximations are better (and slower) than others.

Exploring potential energy surfaces

This is related to the sampling problem which I will address later. Here you convert a $3N$-dimensional sampling problem into a 1, 2 or 3-dimensional one, where you only care about particular nonlinear degrees of freedom (reaction coordinates, collective variables). This gets rid of the exponential scaling, but also needs a certain knowledge of the best/relevant collective variables, which are typically unknown. So this approach is similar in spirit to density functional theory - you convert your problem into a simple one, for which you don't know the exact method and you have to make an educated guess. In terms of sampling nuclear quantum effects, the problem is particularly badly scaling and common methods to estimate typical correlation functions/constants of interest is to either approximate them as simpler classical problems (semi-classical transition state theory), or to convert them into a sampling problem (ring polymer molecular dynamics). The latter is very similar in spirit to diffusion Monte Carlo for electronic structure.

Geometry optimisation

As with all optimisation algorithms, finding a global minimum is an exponentially scaling problem, so to my knowledge, most minimisation algorithms in computational chemistry provide local minima, which scale much better but are also more approximate. In classical computational chemistry you could afford to go one step further and explore much wider conformational space by heating your system up and slowly cooling it down to find some other better minima (simulated annealing). However, as you can see, the result you obtain from this will be highly dependent on chance and convergence will still be exponentially scaling -- there is no way around this.

Sampling

This is one of the biggest unsolved problems in classical computational chemistry. As usual, local sampling is straightforward and typically scales as $3N\log 3N$ (Markov chain Monte Carlo, leapfrog/any other integrator), whereas enhanced sampling either resorts to using collective variables (metadynamics, umbrella sampling) or providing "locally global" sampling, by smoothening kinetic barriers (replica exchange, sequential Monte Carlo). Now, kinetic barriers slow down local sampling exponentially, but the above methods smoothen these linearly, resulting in cheaper locally enhanced sampling. However, there is no free lunch and global convergence will still be exponential, no matter what you do (e.g. protein folding problem).

Partition function calculation

The partition function is a $3N$-dimensional integral (I am going to focus on the classical case, as the quantum one is even more difficult). One way is to try to estimate the partition function (nested sampling, sequential monte carlo), where your convergence will typically scale exponentially but still much, much more efficient than regular quadrature (see exact diagonalisation of the Schrödinger equation, similar problem). This is very difficult, so we typically only try to calculate ratios of partition functions, which are much more nicely behaved. In these cases you can convert the integration problem into a sampling problem (free energy perturbation, thermodynamic integration, nonequilibrium free energy perturbation) and all above sampling issues still apply, so you never really escape the curse of dimensionality, but you get some sort of local convergence, which is still better than nothing :)

So in conclusion, there is no free lunch in computational chemistry and there are various classes of approximations suitable for different problems and in general, the better scaling your problem is, the more approximate and less applicable in general it is. In terms of "best value" nearly exact methods, my vote is on path integral methods (diffusion Monte Carlo, ring polymer molecular dynamics, sequential Monte Carlo), which convert the exponentially scaling problems into polynomially scaling ones (but still with convergence problems) -- although not perfect, at least you won't need all the atoms of the universe to run these and you won't need to know the answer to get the answer, which is sadly an overwhelming problem in many subfields of computational chemistry.

Answer 2

The curse of dimensionality is indeed a huge problem in quantum chemistry, since the possible ways N electrons can occupy K orbitals is a binning problem whose computational cost grows factorially (almost as fast as x^x!) with the size of the system. Moreover, for accurate results you need K>>N in order to account for the so-called dynamical correlation, highlighting the computational challenge of the problem.

A huge breakthrough to the curse of dimensionality was proposed by Walter Kohn: instead of the exponentially difficult problem of describing the anti-symmetric wave function, density functional theory (DFT) shows that it is enough to describe just the electron density n(r), which is just a scalar function. The only problem is that we don't know the exact exchange-correlation functional, which describes how the movement of the electrons is correlated. Still, DFT has been hugely successful in both chemistry and materials science, since in many cases it yields sufficiently accurate results. You can also make DFT linear scaling, if you are smart about the algorithm; however, as far as I am aware, many people are still using the polynomially scaling O(N^3) algorithms since for many systems the lower-order terms are still dominating the cost...

The main problem with DFT is that you don't know the accuracy a priori, and DFT doesn't allow a systematic approach to the exact solution. Wave function based methods to the rescue! It turns out that by being smart, in many cases you can avoid the exponential scaling of exact wave function theory. The exact solution is given by diagonalizing the Hamiltonian in the basis of the possible electronic configurations (given by distributing the N electrons into K orbitals, or the K choose N problem); the size of this Hamiltonian is then (K choose N) x (K choose N) although it is extremely sparse. This is known in chemistry as the configuration interaction problem, and in physics as exact diagonalization.

The problem is extremely hard even for K=N. For example, the 16 electrons in 16 orbitals problem, or (16e,16o), if you are looking at the singlet state you have 8 spin-up and 8 spin-down electrons, yielding (16 choose 8)^2 = 165 million possible configurations. If you go to (18e,18o), you get 2.4 billion configurations. (20e,20o) has 34 billion configurations. (22e,22o) has 500 billion configurations. (24e,24o) has 73 trillion configurations. The (18e,18o) is still practical on a desktop computer, but the (24e,24o) is extremely hard even with a huge supercomputer.

The coupled-cluster method re-expresses the problem with an exponential ansatz, which yields a much more rapidly converging expansion for the wave function; you go down from exponential scaling to polynomial cost - assuming that you don't need to include all possible "excitations". The "gold standard" of quantum chemistry, the CCSD(T) method, scales as O(N^7). It's not cheap, but it yields amazingly accurate results for well-behaved molecules. The density matrix renormalization group could also be mentioned here; it is polynomially scaling for "easy" systems, but reduces to exponential scaling for hard ones....