# A Regularization Method for Inverting Fredholm Integral Equations of the First Kind

## Context
In this note, we are interested in Fredholm equations of the type
$$
G(\tau)=\int_{\mathbb{R}}A(\omega)B(\omega,\tau)d\omega,
$$
where $B(\omega,\cdot)\in C^{\infty}(\mathbb{R})$ is some known function with exponential decay for $\tau\in[\varepsilon,\beta-\varepsilon],$ where $\beta>\varepsilon>0.$ In the context of condensed matter physics, the parameter $\beta$ is understood to have units proportional to inverse temperature.

The basic desire here is, given $G(\tau)$ on the time domain $[0,\tau],$ can we recover the function $A(\omega)?$ The answer is, of course as with all inverse problems, somewhat. It is often the case that, in many body physics, $G(\tau)$ is calculated from monte carlo methods. This naturally introduces some uncertainty in its resolution, and for this reason, it seems that the most widely adopted methods in the field of condensed matter to solve the inverse problem are based Bayesian inference.

In this note, we propose, perhaps, the most basic least squares regression method and find that these bread and butter techniques have metrit in their application to this classical inverse problem

## Numerical Method

We first assume an expansion of $A(\omega)$
$$
A(\omega)=\sum_{n=0}^{\infty}a_nH_n(\omega)
$$

where the Hermite functions of order $n$
$$
H_n(\omega)=\frac{1}{\pi^{1/4}\sqrt{2^nn!}}e^{-n^2/2}h_n(\omega)
$$
are written in terms of the (physicist's) Hermite polynomials given by the recursion
$$
h_{n+1}(\omega)=2 \omega h_n(\omega)-2 n h_{n-1}(\omega)
$$
with base cases $h_0=1,$ $h_1=2\omega$. By truncating the expansion at order $N-1$ and given numerical data $G_{\rm data}(\tau)$, we thus aim to resolve the vector $\mathbf{a}\in\mathbb{R}^N$ through the numerical solution of the inverse problem

$$
\min_{\mathbf{a}\in\mathbb{R}^N}\frac{1}{2}\int_{\varepsilon}^{\tau-\beta}\left[\left(G_{\rm data}(\tau)-\tilde{G}(\tau)\right)^2+\gamma\left(\frac{\partial \tilde{A}}{\partial\omega}\right)^2\right]d\tau
$$

subject to

$$
\tilde{G}(\tau)=\int_{-\Omega(\varepsilon)}^{\Omega(\varepsilon)}\tilde{A}(\omega)B(\omega,\tau)d\omega,\ \ \tilde{A}(\omega)=\sum_{n=0}^{N-1}a_nH_n(\omega)
$$
where the cutoff $\Omega$ is introduced to avoid the issue that $B(\omega,\tau)$ does not have finite support on $\omega\in\mathbb{R}$ for all $\tau\in[0,\beta]$. The cutoff satisfies $\Omega(\varepsilon)\to\infty, \mathrm{as}\ \varepsilon\to0$.

## Numerical Example 1 (Test)
For this problem, we test our methodology on the data set
$$
G_{\rm data}(\tau)=-3(\tau-\beta/2)^2-2(\tau-\beta/2)+.3
$$
with $\beta=5$, time domain cutoff $\varepsilon=1/2$, frequency cutoff $\Omega=20$, and regularization parameter $\gamma=0$. The function $B(\omega,\tau)$ is given by
$$
B(\omega,\tau)=\frac{e^{-\omega\tau}}{2\pi(1+e^{-\beta\omega})}
$$

We perform a study over well we can resolve the function $A(\omega)$ by varying the number of Hermite modes $N$
$N=2$:
$N=4$:
$N=6$:
$N=8$:

The final loss reported in the $N=8$ case was 0.0000853884215417. Here is the numerical implementation in Matlab

% make the spectral domain
N=2^10;
W=20;
w=linspace(-W,W,N);

% inverse temperature parameter
b=5;

% this is where I feed in the data
eps=.5;
t=linspace(eps,b-eps,2^10);
Gdata=-3*(t-b/2).^2-2*(t-b/2)+.3;

% this is a Fermionic distribution
rho=@(w,t)exp(-w*t)./(2*pi*(1+exp(-b*w)));

% this is our initial parameter guess
M=15;
a0=2*rand(1,M)-1;

% fminunc options
options=optimoptions('fminunc','MaxFunctionEvaluations',1e16,...
    'OptimalityTolerance',1e-15,"StepTolerance",1e-15,"MaxIterations",1e6,'Display','off');

% regularization parameter
gamma=1e-16;

% objective function
J=@(a)getLoss(a,w,rho,t,Gdata,gamma);
a=fminunc(J,a0,options);

% this is the learned A
A=buildA(w,a);

% Visualize solution to inverse problem
G=getG(rho,w,t,A);
subplot(2,1,1)
plot(w,A), axis tight, 
xlabel('$\omega$'),ylabel('$A(\omega)$')
subplot(2,1,2)
plot(t,G), hold on, plot(t,Gdata,'--'), axis tight
xlabel('$\tau$'), ylabel('$G(\tau)$')
legend('Learned','Data','Location','best')
hold off
fprintf('This loss is %1.16f.\n',J(a))

function J=getLoss(a,w,rho,t,Gdata,gamma)

% build A(w) from the parameter vector a
A=buildA(w,a);

% this is the Green's function
G=getG(rho,w,t,A);

% try to match to the data
Jloss=trapz(t,abs(G-Gdata).^2)/2;

% introduce regularization
Jreg=gamma/2*trapz(w,gradient(A,w).^2);

% this is the total objective
J=Jloss+Jreg;
end

function G=getG(rho,w,t,A)

% initialize
N=length(t);G=zeros(1,N);

% this is the Green's function at each t
for k=1:N
G(k)=trapz(w,A.*rho(w,t(k)));
end

end

function A=buildA(w,a)

% number of modes
M=length(a);

% get the Hermite basis
H = getHerm(w, M-2);

% build the approximation
A=0;
for j=1:M
    A=A+a(j)/j*H(j,:);
end

end

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

## Numerical Example 2 (Real MCMC Data)
Here $\beta=40$ is used along with the same $B(\omega,\tau)$ as before. I used 10 modes, and chose a fairly large $\varepsilon=5$, with the same frequency cutoff $\Omega=20$, and with quite a small regularization $\gamma=10^{-12}$. Here's an example of what I see when I solve the inverse problem.


A hyperparameter study involving the triple $(\varepsilon,\Omega,\gamma)$ should be performed until satisfactory results are found. Here again is the computer implementation

load('tau_GData.mat'); 

start=2000;
t=tau(start:end-start);
Gdata=GData(start:end-start);

% make the spectral domain
N=2^10;
W=20;
w=linspace(-W,W,N);

% this is a Fermionic distribution
b=40;
rho=@(w,t)exp(-w*t)./(2*pi*(1+exp(-b*w)));

% this is our initial parameter guess
M=10;
a0=2*rand(1,M)-1;

% fminunc options
options=optimoptions('fminunc','MaxFunctionEvaluations',1e16,...
    'OptimalityTolerance',1e-15,"StepTolerance",1e-15,"MaxIterations",1e6,'Display','off');

% regularization parameter
gamma=1e-12;

% objective function
J=@(a)getLoss(a,w,rho,t,Gdata,gamma);
a=fminunc(J,a0,options);

% this is the learned A
A=buildA(w,a);

% Visualize solution to inverse problem
G=getG(rho,w,t,A);
subplot(2,1,1)
plot(w,A), axis tight, 
xlabel('$\omega$'),ylabel('$A(\omega)$'), xlim([-8 8])
subplot(2,1,2)
plot(t,G), hold on, plot(t,Gdata,'--'), axis tight
xlabel('$\tau$'), ylabel('$G(\tau)$')
legend('Learned','Data','Location','best')
hold off
fprintf('This loss is %1.16f.\n',J(a))

function J=getLoss(a,w,rho,t,Gdata,gamma)

% build A(w) from the parameter vector a
A=buildA(w,a);

% this is the Green's function
G=getG(rho,w,t,A);

% try to match to the data
Jloss=trapz(t,abs(G-Gdata).^2)/2;

% introduce regularization
Jreg=gamma*trapz(w,gradient(A,w).^2)/2;

% this is the total objective
J=Jloss+Jreg;
end

function G=getG(rho,w,t,A)

% initialize
N=length(t);G=zeros(1,N);

% this is the Green's function at each t
for k=1:N
G(k)=trapz(w,A.*rho(w,t(k)));
end

end

function A=buildA(w,a)

% number of modes
M=length(a);

% get the Hermite basis
H = getHerm(w, M-2);

% build the approximation
A=0;
for j=1:M
    A=A+a(j)/j*H(j,:);
end

end

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