# I Numerical sign problem software?

1. Dec 18, 2015

### mollwollfumble

I'm an applied mathematician with lots of experience in classical physics but none in quantum physics.
What's the simplest piece of code (eg C++ on Cygwin) that illustrates the https://en.wikipedia.org/wiki/Numerical_sign_problem
in quantum physics?

2. Dec 22, 2015

### mollwollfumble

Have posted question on 3 forums but no replies. Perhaps I'm asking the wrong question. To simplify answer just choose a,b,c,... From this list.
A. There is no such code because you can't get a probability from complex numbers.
B. Ask the author of a journal paper because it's a very advanced topic.
C. Use a full blown lattice qcd package.
D. Can't run on PC because it's too slow.
E. Not in c++. (I'll accept any language).
F. Write from scratch using ideas from "lattice qcd for novices" or similar.
G. The question is ambiguous.

3. Dec 23, 2015

### A. Neumaier

Feynman path integrals (over $e^{iS/\hbar}$) in quantum mechanics (and even more in QFT) are so highly oscillatory that they can hardly be evaluated numerically. The remedy is to work in the Euclidean formulation where (due to the time integration in the action $S$) the constant $i/\hbar$ responsible for the wild oscillations becomes a harmless negative real constant, and analytically continue the results back to real time. In particular, this is used in all of lattice gauge theory, which is feasible numerically only in its Euclidean version.

4. Dec 23, 2015

### atyy

5. Dec 24, 2015

### A. Neumaier

see [above post] - its not simple to code, though. But as a toy example you can easily encode the trapeziodal rule for approximating $\int_{-a}^a e^{i x^2} dx$ in any language of your choice, and observe its misbehavior as $a$ increases!

Last edited by a moderator: Dec 24, 2015
6. Dec 24, 2015

### mollwollfumble

> See p8 of http://www.ncts.ncku.edu.tw/phys/qis/091219/files/07Naoki_Kawashima.pdf
LOL at that. I can see from that how $(3-2)^{10}$ becomes difficult to accurately evaluate using the binomial theorem.
Triangular lattices would be interesting. In my computational fluid dynamics we use finite element or finite volume on lattices where a lattice cell can be any shape; that makes it easy to locally refine it in areas of interest.

Thanks. Will try that. Also will try for $a$ tends to infinity, because I've been looking at equations like this from a pure maths perspective.

7. Dec 25, 2015

### A. Neumaier

You can also try to evaluate $e^{-100}$ using the power series expansion of $e^x$. (But these examples have nothing to do with quantum physics.)

8. Dec 25, 2015

### A. Neumaier

The natural setting for such limits of oscillating integrals in finite dimensions is the Henstock integral.

Last edited: Dec 25, 2015
9. Dec 26, 2015

### mollwollfumble

10. Dec 27, 2015

### A. Neumaier

This still corresponds to a noninteracting situation, hence is physically trivial (and can be solved exactly for $a=\infty$). For a nontrivial interacting version consider $\int e^{ip(x)}dx$ with an even degree polynomial $p(x)$ such as $p(x)=x^2+ \epsilon x^4$ and multivariate generalizations.

11. Dec 29, 2015

### mollwollfumble

I'm still reading about the Henstock integral. I haven't yet convinced myself that it is a) consistent in the sense that the answer is independent of the gauge, b) necessary for QM.

On the pure maths side, I looked at historical systems of infinite numbers. There are a dozen or so, including transseries, Hahn series, Cantor cardinals, Cantor ordinals and Conway's surreals. These are all different, but related. The two systems of relevance to QM are those of Du Bois Reymond's infinitary calculus (1877) and Robinson's hyperreals (1960s). Du Bois Reymond's infinitary calculus ties in with the QM cancellation of infinities and renormalization/regularization. Robinson's hyperreals tie in with the QM concepts of "spooky action at a distance" and the "many-world interpretation".

Starting from Du Bois Reymond, (whose theory generates limits for smoothly diverging sequences) I further developed the theory to include limits of oscillating divergent sequences (and by extension divergent series and Riemann integrals). Hence my interest in the numerical sign problem, I ought to be able to generate unique answers even when the oscillations are divergent. My past experience also includes quadrature and low-discrepancy sequences (for speeding up Monte-Carlo).

I'll follow that advice. I'll try the multivariate generalization.
$\int_{-\infty}^{\infty} e^{i((x_1-1)^2+(x_2-x_1)^2+(x_3-x_2)^2+(x_4-x_3)^2+(x_5-x_4)^2+(x_6-x_5)^2+(x_7-x_6)^2+(1-x_7)^2 +\sum_{i=1}^7p(x_i))}dx_1dx_2dx_3dx_4dx_5dx_6dx_7$
This may take some time, but if I can manage that then I'll be well on the way to handling any integral numerical sign problem. If I succeed then the next hurdle would be in tackling cases which involve a ratio of complex integrals.

12. Dec 30, 2015

### A. Neumaier

If you succeed, you'd tell us something about the techniques you employ; they are interesting even for the toy examples.
But if you want to get really serious about applications in quantum mechanics, you should consider the true Feynman path integral, not the toy versions we have been discussing so far! See, e.g., the following paper:

Albeverio, S., & Mazzucchi, S. (2006). The time-dependent quartic oscillator—a Feynman path integral approach. Journal of Functional Analysis, 238(2), 471-488.

See also the 2008 book Mathematical Theory of Feynman Path Integrals by Albeverio et al..

I don't see how these relate to oscillating integrals. The generalization of the number concept most useful in the present context seems to be nonstandard analysis; see ref [10] in the above paper.

Last edited: Dec 30, 2015
13. Dec 30, 2015

### A. Neumaier

This works in the particular case since the oscillations are equidistant, a property peculiar to the noninteracting situation. A general purpose method cannot rely on such accidents.

14. Jan 8, 2016

### mollwollfumble

Trying that again with equations, and with nonlinearity added to the end.

Here's how the numerical solution of an oscillating divergent integral in QM might work. Equal intervals are used but are not an essential part of the solution method. The first is a toy problem, far too simplistic.
Solve $$\int_{-\infty}^\infty e^{i(x_1^2+...+x_7^2)}dx_1...dx_7$$

The real part equals $$2\int_0^\infty \cos(x_1^2+...+x_7^2)\,dx_1...dx_7$$

Convert to hyperspherical coordinates, that equals $$2S_7\int_0^\infty r^6\cos(r^2)\,dr$$
where $$S_7= 16\pi^3/15$$ is the surface ofr a 7-ball with unit radius.

Substitute $$y=r^2$$ to get the integral component
$$\int_0^\infty y^{2.5}\cos(y)dy$$

To solve this divergent oscillating integral numerically, start by integrating by parts on a finite domain.
$$\int_0^x y^{2.5}\cos(y)dy=x^{2.5}\sin(x)+2.5x^{1.5}\cos(x)-\frac{15}{4}x^{0.5}\sin(x)+\frac{15}{8}\int_0^xy^{-0.5}\sin(y)dy$$

Here's where the study of infinite numbers comes in. It's OK to redefine the infinite limit of a pure oscillating function to be equal to zero. Examples of pure oscillating functions include $$f(x)\cos(g(x))$$ for a range of non-oscillating functions f and g including non-negative powers, finite polynomials, exponential and logarithmic functions. So
$$\lim_{x\to\infty}(x^{2.5}\sin(x)+2.5x^{1.5}\cos(x)-\frac{15}{4}x^{0.5}\sin(x))=0$$
and
$$\int_0^\infty y^{2.5}\cos(y)\,dy=\frac{15}{8}\int_0^\infty y^{-0.5}\sin(y)\,dy$$

The integral on the right hand side is oscillating convergent, so can be directly solved numerically. For faster numerical convergence, do the direct integral only from 0 to π and carry out two more applications of integration by parts to get the term $$\int_{\pi}^{\infty}y^{-2.5}\sin(y)\,dy$$

This can, of course, all be considered moot because the original integral has an analytic solution. But my hope is that a similar numerical method can be used for oscillating integrals when no analytic solution is available.

A word of warning. In general,
$$\lim_{x\to\infty}\frac{\int^xA}{\int^xB}\ne\frac{\int^\infty A}{\int^\infty B}$$
When calculatating expected values, the left hand side must be used not the right hand side.

The nonlinear oscillator is not much more difficult. Start with
$$\int_{-\infty}^\infty e^{i\sum_{n=1}^7p(x_n)}dx_1...dx_7$$
This is separable
$$\int_{-\infty}^\infty e^{i\sum_{n=1}^7p(x_n)}dx_1...dx_7=\prod_{n=1}^7{\int_{-\infty}^\infty e^{ip(x_n)}dx_n}$$

Separate out the real and imaginary component of the one-dimensional integral then, after solving each, multiply the seven complex numbers. For the real component (the imaginary component is treated in an identical way) this is
$$2\int_0^\infty \cos(p(x_n))dx_n$$
Change the variable to $$y=p(x_n)$$
$$dx_n=dy/p'^{-1}$$

When $$p(x_n)=x^2+ax^4$$ that gives
$$p'^{-1}=\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}$$
so we're integrating
$$\int_0^\infty \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$

We might as well split that into two parts, direct numerical integration for
$$\int_0^\pi \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$
and integrate by parts twice (to suppress numerical oscillations) for
$$\int_\pi^\infty \cos(y)/\sqrt{(2/a)*(1+4ay)*(\sqrt(1+4ay)-1)}dy$$

$$\frac{d}{dy}\frac{\sin(y)}{p'^{-1}}=\frac{\cos(y)}{p'^{-1}}+\sin(y)A$$
where
$$A=\frac{4-6\sqrt{1+4ay}}{((2/a)(1+4ay)(\sqrt{1+4ay}-1)^{3/2}}$$

Integrating by parts once more
$$\frac{d}{dy}(-\cos(y)A)=\sin(y)A-\cos(y)\frac{dA}{dy}$$
where
$$\frac{dA}{dy}=\frac{\sqrt{1+4ay}(84ay+33)-120ay-30}{\sqrt{2}\sqrt{1+4ay}(((1+4ay)^{3/2}-1-4ay)/a)^{5/2}}$$

Using our definition of the limit of pure oscillatories at infinity.
$$\lim_{y\to\infty}\sin(y)/p'^{-1}=\lim_{y\to\infty}\cos(y)A=0$$

Then just integrate numerically in 1-D out to whichever multiple of π is sufficient to give the accuracy desired, for two decimal places, 6π may be large enough to use instead of ∞.

In the general nonlinear oscillator where $$p(x_n)\ne x^2+ax^4$$ the procedure is the same, but A and dA/dy are obtained by numerically differentiating 1/p'-1 rather than analytically differentiating it.

15. Jan 8, 2016

### A. Neumaier

Thanks for providing the details.

Yes, the separable case is not very difficult since it factors into univariate integrals. Also the general quadratic case is not difficult since after a linear transformation it becomes separable.

But in the multivariate interacting case, the exponent has the form $i(x^TKx/2m+V(x))$ where $m>0$ is a mass parameter, the potential $V(x)$ is separable and bounded below but $K$, the stiffness matrix of the oscillating system, is symmetric positive definite and non-diagonal. The tridiagonal case is already nontrivial since there is no simple transformation to separable form.

Thus you need to extend your current techniques.

16. Jan 8, 2016

### A. Neumaier

In practice, expectation values are most often calculated via derivatives of the free energy $F=i\log Z$, where $Z$ denotes the integral, after approximating the free energy by an appropriate analytic expression.

17. Jan 9, 2016

### mollwollfumble

Got it. I said before that I wanted to be able to numerically integrate the path integral
$$I=\int_{-\infty}^\infty e^{i((x_1-1)^2+(x_2-x_1)^2+(x_3-x_2)^2+(x_4-x_3)^2+(x_5-x_4)^2+(x_6-x_5)^2+(x_7-x_6)^2+(1-x_7)^2+\sum_{n=1}^7p(x_n))}\,dx_1dx_2...dx_7$$
Write this as
$$I=\int_{-\infty}^\infty e^{iq}\,dx_1dx_2...dx_7$$
The antisymmertic parts of $q$ don't contribute to $I$ so can be omitted, leaving
$$I=\int_{-\infty}^\infty e^{i(2+\sum_{n=1}^7(2x_n^2+p(x_n)))}\,dx_1dx_2...dx_7$$
This is separable, giving
$$I=e^{2i}(\int_{-\infty}^\infty e^{i(2x^2+p(x))}\,dx)^7$$
Only an integral in 1-D is needed. Change variable using $y=2x^2+p(x)$ to give
$$\int_{-\infty}^\infty e^{i(2x^2+p(x))}\,dx=2\int_0^\infty e^{iy}(dx/dy)dy$$
Here for simple $p(x)$ the formula for $dx/dy$ can be found analytically and for messy $p(x)$ the formula for $dx/dy$ can be found numerically. This then resolves to the case of the nonlinear oscillator described above and the same methods are used. In a nutshell, keep integrating by parts until the oscillations are suppressed enough to get an accurate answer using standard quadrature methods. The differentials for the integral by parts $$d^2x/dy^2,\,d^3x/dy^3$$ etc. can be evaluated either analytically or numerically, depending on the original $p(x)$.

This method also goes a long way towards solving the more general integral
$$\int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x})}\,d\mathbf{x}$$
when $f(\mathbf{x})$ is separable and $g(\mathbf{x})$ does not contain symmetric cross correlations such as $x^2y^2$.
Here, it's $f(x)(dx/dy)$ that is differentiated again for each iteration of integration by parts.

Excellent. I'll see what I can do on that basis, as said before this may take some time. (sorry about loss of equations in quote, I don't know why).

18. Jan 9, 2016

### A. Neumaier

There is no antisymmetric part in a non-diagonal quadratic form. $\int dx^n e^{ix^TKx}$ is proportional to $(\det K)^{-1/2}$, but your formula with $p(x)=0$ gives something different!

By the way, at PF, use two hashes # in place of each \$. Morever, when quoting formulas, you need to edit out the artifacts created by the system! (You can edit it even days afterwards).

19. Jan 10, 2016

### mollwollfumble

Oops. Thank you very very much for correcting my error. (It's not by accident that "fumble" is part of my username).

All along, from the very beginning, I've had three ultimate goals in mind. Three general-purpose numerical methods for evaluating the non-separable asymmetric n-dimensional integral
$$\int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x})}\,d\mathbf{x}$$
(And by extension, expected values using the ratio of two such integrals).
All three start the same way. Find the minimum of $g(\mathbf{x})$. Call the xector at that minimum $\mathbf{x}_0$. This minimum can be found either analytically, or by some variation of the conjugate gradient method, Polak-Ribiere if analytical derivatives of $g(\mathbf{x})$ are available or Powell's method if they are not.
From here on evaluate
$$I=\int_{-\infty}^\infty f(\mathbf{x})e^{ig(\mathbf{x}-\mathbf{x}_0)}\,d\mathbf{x}$$

Method 1, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ is dominated by the quadratic term $(\mathbf{x}-\mathbf{x}_0)^TK(\mathbf{x}-\mathbf{x}_0)$. This has high speed and low accuracy.

The conjugate gradient method generates, as a byproduct, the principal axes of the ellipsoid in n dimensions around the minimum point. Using that, or the principal axes found analytically, numerically integrate $I$ along each of the n principal axes in turn, using integration by parts to reduce oscillations. Multiply the resulting n integrals to get the final answer $I$.

Method 2, suitable when $g(\mathbf{x}-\mathbf{x}_0)$ has a significant quadratic term component. This has intermediate speed and intermediatre accuracy.
Numerically integrate $I$ along $n(n+1)/2$ one dimensional vectors given in cartestian coordinates by $x_1-x_0, x_2-x_0, x_2-x_1, x_3-x_0, x_3-x_1, x_3-x_2, x_4-x_0,...,x_n-x_{n-1}$, using integration by parts in each case to reduce oscillations.

From this derive an approximate matrix $K$ using $x_ix_j=(1/2)((x_i-x_0)^2+(x_j-x_0)^2-(x_i-x_j)^2)$ for off-diagonal terms $i\ne j$. Evaluate $I$ from $\sqrt{\det(K)}$.

Method 3, suitable for nasty $f$ or $g$. This is slow and has an accuracy that ranges from poor to good depending on the number of vectors chosen.

Define on the surface of a unit sphere in n-D, the endpoints of $m>>n$ unit vectors using an algorithm that spaces these endpoints $\mathbf{x}_1$ to $\mathbf{x}_m$ as far apart from each other as possible. One way to do this would be to define a grid similar to that used to define the nodes of a geodesic dome. A second way would be use use additive quasirandom numbers (see wikipedia for algorithm on a hypercube and then map to the hypersphere surface using the Sag-Szekeres transformation).

Then use a higher-dimensional extension of the method of Veronoi polyhedra (or a faster approximation) to assign a weight $w_1,...,w_m$ to each vector, endpoints further away from their neighbours are allocated a greater weight.

Numerically integrate $I$ along the $m$ one-dimensional vectors, using integration by parts in each case to reduce oscillations.

Evaluate $I$ as the nth root of the product of the $m$ integrations, each weighted using the assigned weights $w_1,...,w_m$.

A correction for kurtosis is necessary, there's no correction when the 1-D distribution $\int_{-\infty}^\infty f(x)e^{-g(x)}\,dx$ is Gaussian (zero excess kurtosis). When the excess kurtosis is -1.2 (uniform distribution) the correction is $S_n/(n2^n)$ and when the excess Kurtosis is 3 (Laplace distribution) the correction is $S_n(n-1)!/2^n$ where $S_n$ is the surface of an n-ball.

20. Jan 10, 2016

### A. Neumaier

This is a variant of the saddle point method used to find approximations to the integral. it is much used in quantum mechanics. The assumption is not satisfied when there are quartic terms, which always dominate the quadratic term for large $x$. It is often used nevertheless as an asymptotic series.

I don't see why this should give an approximation of the integral. Please explain for $n=2$.

This is best done via spherical $t$-designs and related constructions; see my paper
Measures of strength 2e, and e-optimal designs, Sankhya 54 (1992), 299-309.
and its citations.

Again I don't see how you can replace the integration over $R^n$ by $m$ univariate integrations.