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.

# 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]

# 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]

# 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:

• CartesianIndices

• iterator

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?

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

matgen (generic function with 1 method)

# 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)))]

2-element Vector{Float64}:
 5.099433819233062
 5.1031036307982856

# 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)))]

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:

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]
        )]
    ]

2-element Vector{Vector{Float64}}:
 [-0.010475500501098162, -0.010915195936659329]
 [-0.0019012438540132104, -0.0016228593278076148]