Advice on the regularisation of a linear problem

Question

I'm numerically inverting an integral transform using a method suggested by a scicomp user from an earlier question. The problem is as follows:

I wish to estimate $f(x)$ for a given $F(y)$, both of which are expressed as a linear combination of some set of orthogonal basdis functions $\phi_n$:

$$ f(x) = \sum_{i}^{M} c_i \phi_i(x) \quad \quad \quad F(y) = \sum_{j}^{M} b_j \phi_j(y) $$

The coefficients needed to reconstruct $f(x)$ are given by the vector $\mathbf{c}$, which can be obtained by solving the linear linear system:

$$ \mathbf{b} = \mathbf{A}\mathbf{c} $$

This method works well in some theoretical test cases, but when applying it to real data I find that the solution to this system produces an $f(x)$ which is highly oscillatory, with much more power in the higher order basis functions. Adding small amounts of random noise to test cases which do work also produces incorrect oscillatory solutions.

I know the following about the true solution for $f(x)$:

• $f(x)$ is not oscillatory
• $f(x)$ has a smooth, decaying form - often similar to a half-Cauchy distribution
• $f(x)$ is positive everywhere (it's a continuous probability distribution)

I've done some reading around the problem and I think I need to use regularisation to constrain the solution of the linear system such that it produces the correct $f(x)$, but I've seen there are many types of regularisation.

Can anyone give me advice on which approaches to regularisation would be appropriate and how to apply them?

Thanks for your time and help!

Answer 1

Brian already provides you with the general outline of how regularization works (apart from the fact that choosing the regularization parameter $\lambda$ is difficult). However, you can make the regularization work better if you encode your knowledge that you expect the highly oscillatory terms to be small, whereas the low frequency terms can be large.

The way you do this is to recognize that when you solve $$ \min \|Ac-b\|_2^2 + \lambda\|c\|_2^2 $$ you essentially penalize large coefficients $c_i$. In other words, what you are doing is trying to find those coefficients $c$ so that $Ac-b$ is small but also so that $c$ is small (and the balance between the two terms is given by $\lambda$). But that's not what you want, according to your problem description: you only want the high frequency terms to be small. Consequently, you may want to minimize the following instead: $$ \min \|Ac-b\|_2^2 + \sum_i\lambda_i|c_i|^2 $$ where $\lambda_i$ is an increasing sequence so that higher coefficients are penalized more. What form you choose for the $\lambda_i$ is something you need to decide based on what your expectations are for the sizes of the coefficients. An example may be something along the lines $$ \lambda_i = \bar\lambda \sqrt{i}. $$

There is a general theory on how regularization should be chosen to match what we expect solutions to look like. If you want to spend a significant amount of time on these issues, I suggest to read the books by Tarantola ("Inverse Problem Theory") or by Kaipio and Somersalo ("Statistical and computational inverse problems").

Answer 2

Because your basis functions are orthogonal polynomials (and assuming that they've been properly normalized), $\| f \|_{2}=\| c \|_{2}$. Thus you could regularize the solution of the linear system of equations in the 2-norm and effectively regularize $f$ in the 2-norm as well.

To do this, simply minimize $\| Ac-b\|_{2}^{2}+ \lambda \| c \|_{2}^{2}$. This can be set up as a least squares problem

$\min \left\| \left[ \begin{array}{c} A \\ \sqrt{\lambda} I \end{array} \right] b- \left[ \begin{array}{c} b \\ 0 \\ \end{array} \right] \right\|_{2}^{2} $

The normal equations for this damped least squares problem give the solution

$\hat{c}=(A^{T}A+\lambda I)^{-1}A^{T}b$

You could easily add a constraint to force $f$ to integrate out to 1, if this isn't already implied by the equations $Ac=b$.

However,

An important aspect of your problem is that $f$ is inherently nonnegative (it could be 0 at some points, so its not necessarily positive everywhere.) Unfortunately, the orthogonal basis functions that you're using aren't nonnegative everywhere, and there's no reasonable way to enforce nonnegativity of the solution in terms of the coefficients in the orthogonal polynomial expansion of $f$.

Thus you might want to consider getting rid of the orthogonal polynomials and solving the problem in terms of basis functions that allow you to enforce the nonnegativity of $f$.

You might also consider expressing your problem in terms of the CDF rather than the pdf of your probability distribution.