basic Gaussian Integrals¶
The final form (for now)¶
$\int dx_1 dx_2 \ldots dx_N \, e^{-\vec{x}^t A \vec{x} + \vec{h}^t \vec{x}} = \sqrt{\frac{\pi^N}{det A}} \, e^{\frac{1}{4} \, \vec{h}^t \, A^{-1} \, \vec{h}}$.
Dangerously close to field theory.
In [1]:
# 1D Gaussian
N = 40
rv = range(-1.0, 1.0, length=N)
dr = 2.0/(N-1)
# map (-1, 1) to (-Inf, Inf)
xvf = r -> r/(1-r^2)
xwf = r -> (1+r^2)/(1-r^2)^2
f = x -> exp(-x^2)
ret = 0.0
for rr in rv[1:N-1]
ret += dr * xwf(rr+dr/2) * f(xvf(rr+dr/2))
end
println([ret, sqrt(pi)])
[1.7724472204698454, 1.7724538509055159]
In [2]:
# using Gaussian Quadrature
using QuadGK
N = 40
rv, rw = QuadGK.gauss(N) # -1:1
# map (-1, 1) to (-Inf, Inf)
xvf = r -> r/(1-r^2)
xwf = r -> (1+r^2)/(1-r^2)^2
f = x -> exp(-x^2)
ret = 0.0
# no need for dr
ret = sum(@. rw * xwf(rv) * f(xvf(rv)) )
println([ret, sqrt(pi)])
[1.772453426851392, 1.7724538509055159]
In [3]:
# 2D Gaussian
N = 40
rv = range(-1.0, 1.0, length=N)
dr = 2.0/(N-1)
xvf = r -> r/(1-r^2)
xwf = r -> (1+r^2)/(1-r^2)^2
f = (x, y) -> exp(-x^2-y^2)
ret = 0.0
for rr1 in rv[1:N-1], rr2 in rv[1:N-1]
ret += dr^2 * xwf(rr1+dr/2) * xwf(rr2+dr/2) * f(xvf(rr1+dr/2), xvf(rr2+dr/2))
end
println([ret, sqrt(pi)^2])
[3.141569149351278, 3.1415926535897927]
use of iterator¶
either the more modern construction: CartesianIndices
or with Iterators
For details, consult the manual:
The final form (for now)¶
$\int dx_1 dx_2 \ldots dx_N \, e^{-\vec{x}^t A \vec{x} + \vec{h}^t \vec{x}} = \sqrt{\frac{\pi^N}{det A}} \, e^{\frac{1}{4} \, \vec{h}^t \, A^{-1} \, \vec{h}}$.
Dangerously close to field theory.
How to generate a random positive definite matrix?¶
In [4]:
using LinearAlgebra
function matgen(ndim)
# making a positive definite matrix
# get an orthogonal matrix U from QR factorization
R = rand(ndim, ndim)
U = qr(R).Q
# A = diagm(rand(ndim))
A = diagm([(i1 / ndim * pi) for i1 = 1:ndim])
A = U' * A * U
return A
end
Out[4]:
matgen (generic function with 1 method)
In [5]:
# demo: using CartesianIndices
ndim = 5
# get a positive definite matrix
A = matgen(ndim)
N = 20
# rv goes from -1 -> 1
small = 1E-10
rv = range(-1.0+small, 1.0-small, length=N)
dr = 2.0/(N-1)
# xvf = rr -> atanh(rr)
# xwf = rr -> 1 / (1 - rr^2)
xvf = r -> r/(1-r^2)
xwf = r -> (1+r^2)/(1-r^2)^2
xv, xw = (xvf.(rv), xwf.(rv))
# nice device to generate (1:N, 1:N, 1:N, ...) # ndim of them
imat = CartesianIndices(ntuple(i1 -> N, ndim))
ret = 0.0
for II in imat
# .I to extract the indices
xv1 = [xv[i1] for i1 in II.I]
xw1 = [xw[i1] for i1 in II.I]
ret += dr^ndim * prod(xw1) * exp(-xv1' * A * xv1)
end
[ret, sqrt(pi^ndim / (det(A)))]
Out[5]:
2-element Vector{Float64}: 5.099433819233062 5.1031036307982856
In [6]:
# demo: using iterator
ndim = 5
# get a positive definite matrix
A = matgen(ndim)
N = 20
# rv goes from -1 -> 1
small = 1E-10
rv = range(-1.0+small, 1.0-small, length=N)
dr = 2.0/(N-1)
# xvf = rr -> atanh(rr)
# xwf = rr -> 1 / (1 - rr^2)
xvf = r -> r/(1-r^2)
xwf = r -> (1+r^2)/(1-r^2)^2
xv, xw = (xvf.(rv), xwf.(rv))
# iterator to generate a grid of all pairs
pair_collection = Iterators.product([zip(xv, xw) for i1 in 1:ndim]...)
ret = 0.0
for pair in pair_collection
xv = [pair[i1][1] for i1 in 1:ndim]
xw = [pair[i1][2] for i1 in 1:ndim]
ret += dr^ndim * prod(xw) * exp(-xv' * A * xv)
end
[ret, sqrt(pi^ndim / (det(A)))]
Out[6]:
2-element Vector{Float64}: 5.099822580147237 5.103103630798289
Integration with Monte Carlo¶
Monte Carlo method is perfect for high dimensional integral, here is an illustration:
In [7]:
ndim = 8
A = matgen(ndim)
Nmax = 2888888
# 0:1 -> -Inf:Inf
xv = r -> tan(pi * (r-0.5))
xw = r -> pi/cos(pi * (r-0.5))^2
ret0 = 0.0
ret1 = 0.0
ret2 = 0.0
for ii in 1:Nmax
rr = rand(ndim)
xv_arr = xv.(rr)
xw_arr = xw.(rr)
basic = prod(xw_arr) / Nmax * (
exp(-xv_arr' * A * xv_arr)
)
ret0 += basic
ret1 += basic * xv_arr[2] * xv_arr[3]
ret2 += basic * (xv_arr[1] * xv_arr[2] * xv_arr[3] * xv_arr[4])
end
[
[ret1 / ret0, (1 / 2) * inv(A)[2, 3]],
[ret2 / ret0, (1 / 4) * (
inv(A)[1, 2] * inv(A)[3, 4] +
inv(A)[1, 3] * inv(A)[2, 4] +
inv(A)[1, 4] * inv(A)[2, 3]
)]
]
Out[7]:
2-element Vector{Vector{Float64}}: [-0.010475500501098162, -0.010915195936659329] [-0.0019012438540132104, -0.0016228593278076148]