# Numerical methods for inverting integral transforms?

I'm trying to numerically invert the following integral transform:

$$F(y) = \int_{0}^{\infty} y\exp{\left[-\frac{1}{2}(y^2 + x^2)\right]} I_0\left(xy\right)f(x)\;\mathrm{d}x$$

So for a given $F(y)$ I need to approximate $f(x)$ where:

• $f(x)$ and $F(y)$ are real and positive (they are continuous probability distributions)
• $x,y$ are real and positive (they are magnitudes)

I have a very messy and brute force method for doing this at the minute:

I define $f(x)$ and the spline over a series of points, the values of the splined points are 'guessed' by random sampling, which yields a predicted $F(y)$. A basic genetic algorithm I wrote up minimises the difference between the predicted and measured $F(y)$ array. I then take the $f(x)$ which the algorithm converges to as my answer for the inversion.

This approach works fairly well for some simple cases, but it feels messy to me and not particularly robust.

Can anyone give me guidance on better ways of solving this problem?

Thanks for your time & help!

[x-posted at computerscience]

A fairly simple method would be to choose a basis in function space and convert the integral transformation to a matrix. Then you can just invert the matrix.

Mathematically, here's how that works: you need some set of orthonormal basis functions $T_i(x)$. (You can get away without them being normalized too, but it's easier to explain this way.) Orthonormal means that the inner product $\langle T_i,T_j\rangle = \delta_{ij}$, where

$$\langle T_i,T_j\rangle \equiv \int_a^b W(x) T_i(x)T_j(x)\,\mathrm{d}x = \delta_{ij}\tag{1}$$

Here $W(x)$ is some weight function. That and the limits $a$ and $b$ are tied to your choice of $T_i$. Once you pick which set of basis functions to use, you can hard code the limits and weight function into your program.

Using the orthonormality, you can express any function, such as $f(x)$ and $F(y)$, as linear combinations of these basis functions:

\begin{align} f(x) &= \sum_{i} c_i T_i(x) & F(y) &= \sum_{j} C_j T_j(y) \tag{2} \end{align}

where the coefficients are calculated as

\begin{align} c_i &= \langle f,T_i\rangle = \int_a^b W(x) f(x)T_i(x)\,\mathrm{d}x \tag{3}\\ C_j &= \langle F,T_j\rangle = \int_a^b W(y) F(y)T_j(y)\,\mathrm{d}y \tag{4} \end{align}

You can verify that these expressions are consistent with the definitions of the coefficients, eq. (2), and the orthonormality, eq. (1).

Now, compute the transform of each of the basis functions; let's call it $\tilde{T}_i(y)$.

$$\tilde{T}_i(y) \equiv \int_0^\infty y \exp\biggl[-\frac{1}{2}(y^2 + x^2)\biggr]I_0(xy)T_i(x)\,\mathrm{d}x$$

$\tilde{T}_i(y)$ is a function, and so you can express it as a linear combination of the basis functions just as we did with $f(x)$ and $F(y)$:

$$\tilde{T}_i(y) = \sum_k A_{ik}T_k(y)$$

where the matrix elements $A_{ik}$ are determined the same way we found $c_i$ and $C_j$ above:

$$A_{ik} = \langle \tilde{T}_i, T_k\rangle = \int_a^b W(y) \tilde{T}_i(y)T_k(y)\,\mathrm{d}y\tag{5}$$

In practice, this is a rather icky double integral, but you only have to do it once (ever) for each combination of $i$ and $k$. You can do the integrals numerically and then hard-code the resulting values in your program. (Side note: With a clever choice of $T_i(x)$ and $W(x)$, you might be able to make it so the integral can be done symbolically. Whether this is possible depends on your transform. You can do it with the Fourier transform, but I'm inclined to think it's not possible for the transform you're asking about here.)

In terms of the matrix elements $A_{ik}$ and coefficients $c_i$ and $C_j$, the relation between $f(x)$ and $F(y)$ reduces to a linear system

\begin{align} \overbrace{\sum_{j} C_j T_j(y)}^{F(y)} &= \int_0^\infty y \exp\biggl[-\frac{1}{2}(y^2 + x^2)\biggr]I_0(xy)\overbrace{\sum_{i} c_i T_i(x)}^{f(x)}\,\mathrm{d}x\\ &= \sum_{i} c_i \int_0^\infty y \exp\biggl[-\frac{1}{2}(y^2 + x^2)\biggr]I_0(xy)T_i(x)\,\mathrm{d}x\\ &= \sum_{i} c_i \sum_{k} A_{ik}T_k(y) \end{align}

Given the orthogonality of the basis functions, you can isolate any particular coefficient $C_{\ell}$ by taking the inner product of both sides with $T_{\ell}$:

\begin{align} \left\langle \biggl(\sum_{j} C_j T_j\biggr), T_\ell\right\rangle &= \left\langle \biggl(\sum_{i} c_i \sum_{k} A_{ik}T_k\biggr), T_\ell\right\rangle \\ \int_a^b W(y) \sum_{j} C_j T_j(y)T_\ell(y)\,\mathrm{d}y &= \int_a^b W(y) \sum_{i} c_i \sum_{j} A_{ik}T_k(y)T_\ell(y)\,\mathrm{d}y \\ \sum_{j} C_j \int_a^b W(y) T_j(y)T_\ell(y)\,\mathrm{d}y &= \sum_{i} c_i \sum_{k} A_{ik}\int_a^b W(y) T_k(y)T_\ell(y)\,\mathrm{d}y \\ \sum_{j} C_j \delta_{j\ell} &= \sum_{i} c_i \sum_{k} A_{ik}\delta_{k\ell} \\ C_\ell &= \sum_{i} c_i A_{i\ell} \end{align}

Of course $\ell$ is just a dummy index so I'll go back to calling it $C_j$.

This is just a linear algebra problem. $C_j$ is a component of a vector, as is $c_i$, and $A_{ij}$ are components of a matrix. You can calculate $c_i$ from the function you're trying to transform using eq. (3), and you know $A_{ij}$ from the one-time calculation you did for this particular integral transform, eq. (5), so you can get the $C_j$s by doing a matrix multiplication (which computers are very good at), and then reconstruct the function $F(y)$ from those using eq. (2).

Conversely, to do the inverse transform, you start with the function $F(y)$, calculate the $C_j$s from it using eq. (4), and then you would need to solve the linear system

$$C_j = \sum_{i} c_i A_{ij}$$

This can be done by multiplying both sides by the inverse of the matrix $A$, but in practice there are more efficient ways to do it. Use the linear system solver of whatever linear algebra library you have at hand.

Note that everything I've written so far has left the limits of the sums over $i$, $j$, etc. unspecified. In practice, you will need to choose some limits, say $1$ to $N$ (you choose $N$), such that any $f(x)$ you may encounter can be sufficiently well approximated by a linear combination of $T_1(x),\ldots,T_N(x)$, and also a range like $1$ to $M$ such that any $F(y)$ can be sufficiently well approximated by a linear combination of $T_1(y),\ldots,T_M(y)$. It's probably easiest to choose $M = N$. Here, what "sufficiently well" means is up to you. Generally, the larger you make $M$ and $N$, the better your approximation will be, but the more memory (and time) you'll need to do the calculation. Remember that you need to calculate $N$ coefficients $c_i$, and for the matrix $A$ you need to have $M\times N$ coefficients calculated, from $A_{11}$ up to $A_{NM}$.

For functions defined on a finite interval which can be linearly rescaled to $[-1,1]$, the Chebyshev polynomials are a common choice for $T_i$. There is a sense in which Chebyshev polynomials produce the most accurate approximation of a given function using a fixed number of terms, which makes them especially well suited for this kind of application. If you can truncate the domain of your probability distributions, you might be able to use those. They are associated with the weight $W(x) = \frac{1}{\sqrt{1 - x^2}}$ and the limits $a = -1$ and $b = 1$. (Note that in the form they're often given, the normalization is such that $\langle T_i, T_j\rangle = \delta_{ij}\pi/2$ for $i = j\neq 0$ and $\langle T_0, T_0\rangle = \pi$, so you may have to account for this in your code.)