# Proxy Integration
We are interested in building a "proxy" for purposes of integrating nasty functions in 3D dimensional space that come up in physics contexts. One might think that we are talking about a numerical integration scheme. This is different and can be viewed as something like a low rank approximation of functions. To be precise about the problem we want to solve, let's suppose you are given a function $f:\mathbb{R}^2\to\mathbb{\mathbb{C}}$, and suppose you were interested in knowing the function $g:\mathbb{R}\to\mathbb{C}$ such that $$g(k)=\int_{\Omega}f(x,k)dx,$$

where $\Omega\subset\mathbb{R}$. This, in general, can be computed numerically for as many values of $k$ one has the computational resources for. However, it is clear that having access to $g(k)$ in closed form has obvious advantages.

Let's take the following example of calculating
$$
f(k)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-x^2}e^{ikx}dx.
$$
This is nothing more than computing the Fourier transform of a Gaussian, which is another Gaussian, that is, $f(k)=\frac{1}{\sqrt{2}}e^{-k^2/4}$. This, of course, can be numerically computed by a discrete Fourier transform and quickly by using the Tukey and Cooley fast Fourier transform (FFT) algorithm.

But suppose we didn't know how to calculate this integral in closed form and we did not have access to the FFT algorithm (because we want to build something that is integrand agnostic and not tailored to complex exponentials). We assume we can only compute the integral at a finite number of points in $\Omega$. Moreover, suppose these points are chosen at random, since we have no idea what the landscape actually looks like. I will exploit just one fact; $f(k)\in L^2(\mathbb{R})$ and, without loss of generality, it has unit $L^2(\mathbb{R})$ norm. In this case, I know that I can approximate $f(k)$ using a Hermite basis in the sense that
$$
\lim_{n\to\infty}\int_{\mathbb{R}}|f(k)-H_n(k)|^2dk=0,
$$
where the Hermite functions of order $n$
$$
H_n(k)=\frac{1}{\pi^{1/4}\sqrt{2^kk!}}e^{-k^2/2}h_n(k)
$$
are written in terms of the (physicist's) Hermite polynomials given by the recursion
$$
h_{n+1}(x)=2 k h_n(k)-2 n h_{n-1}(k)
$$
with base cases $h_0=1,$ $h_1=2k$.

Now we aim to solve the following regression problem:
$$
\min_{\xi\in\mathbb{R}^N}J[\xi]=\frac{1}{N_{\rm samples}}\sum_{j=1}^{N_{\rm samples}}\left|f(k_j)-\sum_{n=1}^N\xi_nH_n(k_j)\right|^2
$$

The idea is that by solving this optimization problem, the Hermite function of order $n$ with computed weights $\xi$ will serve as a reduced order model of the integral. Then, in whatever physics applications, we have an interpretable model of the integral that can aid in furthering some analytical calculation.

To solve the optimization problem, I used a standard line search method implemented by MATLAB's built in fminunc function. Here's are numerical examples of a twelve term approximation of $f(k)=\sqrt[4]{\frac{2a}{\pi}}e^{-a(k-\mu)^2}$ with 20 sampling points and $a=.8,\ \mu=2$


This is not reliable. However, If I double the sampling rate to 40 points, the approximation is fairly consistent and looks like the second image (the proxy is good). There is a subtle balance between the number of Hermite terms one should take (under and over fitting) and the sampling rate. I believe there is interesting mathematics that can be done even at this elementary level (something akin to the Shannon-Nyquist theorem could be leveraged). I believe we may be able to reframe this as a compressed sensing problem as well. Of course, we should keep in mind how this will scale to larger classes of function spaces as well as increasing the spatial dimension.

## Matlab Code
Just 61 lines

% wavenumber grid
l = 8;
N = 2^10;
k = linspace(-l, l, N+1);

% get the Hermite basis (each basis vector is stored as a row)
numHerm = 10;
H = getHerm(k, numHerm);

% get the test f(k)
a = .8; mu = 2;
fk = exp(-a * abs(k - mu).^2);
fk = fk / sqrt(sqrt(pi / 2 / a));

% get the randomly sampled points sorted in increasing order
ind = randi(N, 40, 1);
ind = sort(ind);

% define the loss function for the regression
J = @(xi) getLoss(H, fk, xi, numHerm, ind);

% solve it
bestxi = fminunc(J, rand(1, numHerm + 2));

% build the proxy and plot it
Hk = getSuper(H, bestxi, numHerm);
plot(k, fk), hold on, plot(k, Hk, '--')
plot(k(ind), fk(ind), 'k*'), axis tight, hold off
legend('Exact $f(k)$', 'Proxy Hermite Basis $H_n(k)$', 'Sample Points', 'Location', 'best')
xlabel('$k$')

function H = getHerm(x, numHerm)
    % gets you the Hermite basis and stores it in a matrix as the rows
    N = length(x);
    % base cases
    H(1, :) = ones(1, N);
    H(2, :) = 2 * x;
    % get the Hermite polynomials
    for k = 1:numHerm
        H(k + 2, :) = 2 * x .* H(k + 1, :) - 2 * k * H(k, :);
    end
    % normalize and hit them with a Gaussian
    for k = 1:numHerm + 2
        H(k, :) = H(k, :) .* exp(-x.^2 / 2) / sqrt(2^(k - 1) * factorial(k - 1)) * pi^(-1/4);
    end
end

function J = getLoss(H, fk, eta, numHerm, ind)
    % get the superposition
    Hk = getSuper(H, eta, numHerm);
    % compute the loss as the empirical mean of the l2 distance
    J = sum(abs(fk(ind) - Hk(ind)).^2) / length(ind);
end

function Hk = getSuper(H, eta, numHerm)
    % gets me that superposition
    Hk = 0;
    for k = 1:numHerm + 2
        Hk = Hk + eta(k) * H(k, :);
    end
end