# Objective function scaling in an Inverse Problem

I am trying to solve a large scale inverse problem using the Bayesian formulation. To estimate the Maximum a Posteriori Estimation (MAP) solution I will have to minimize the following objective function:

$F(m) = F_d + F_\text{prior} = \frac{(f(m) - d)^T(f(m) - d)}{\sigma^2} + \frac{(m - m_\text{prior})^T(m - m_\text{prior})}{\alpha^2}$

where d is the observed data, and $\sigma$ is the uncertainty in $d$. $m$ is the optimization parameter, $\alpha$ represents the confidence in prior.

In the current problem setup, the number of data points, $d$, are $O(n)$ whereas the number of parameters, $m$, are $O(n^2)$. Everything else is dimensionless, therefore order of magnitude of f(m) is same as m. This results to an objective function which is inherently biased towards the prior term($F_\text{prior}$), unless $\sigma << \alpha$, further resulting in poor convergence of $F_d$ in the optimization process.

In such a case can I scale $F_d$ and $F_\text{prior}$ by the number of terms they contain? From what I understand such scaling will change the interpretation of $\sigma$ and $\alpha$.

Note that I get better "match" between $f(m)$ and $d$ without the prior term. But I need to include the prior in order to get the bounds on the posterior solution.

• What do you observe if you divide the regularizing term by $n$? The choice of $\alpha$ is almost arbitrary, depending only on your idea of what a good regularized solution would look like, and the relative scaling of the two terms also depends on the magnitude of $f,d$. – Kirill Jul 16 '15 at 23:18
• I never tried dividing the regularization term with n, but it is somewhat equivalent to a large $\alpha$. The choice of $\alpha$ cannot be arbitrary as the posterior distribution in my problem depends entirely on $\alpha$ except for a small region where $d$ is specified, further the solver, $f(m)$ does not converge for a high $\alpha$ when calculating posterior bounds. Order of magnitude of f ~ d ~ m. – maverick Jul 17 '15 at 2:05
• Regularization is usually done to prevent over-fitting; you seem to be saying that you can't fit your model at all. Re $\alpha$: with regularization you always get a family of solutions, one for each $\alpha$, and have some other method for picking which value of $\alpha$ works for you. – Kirill Jul 17 '15 at 2:08
• @ChristianClason Sorry, unless I'm missing something, what you said doesn't contradict it being arbitrary. It can be set to any value in this formulation, and there will a corresponding solution. It's up to OP to pick the most useful value (the most useful prior). OP talks about scaling the $F$'s relative to each other, which means the prior is not determined in advance, so the value of $\alpha$ isn't fixed. – Kirill Jul 17 '15 at 8:12
• @ChristianClason I understand that, but OP has written down the prior $N(m_{\mathrm{prior}},\alpha^2)$, and AFAICT has a free choice of $\alpha$ (corresponding to free choice of the prior). I think "rescaling afterwards" wouldn't make much sense, since it just means picking a different prior to start with (which is fine if the prior isn't pre-determined somewhere else). – Kirill Jul 17 '15 at 8:35

First, a disclaimer: I'll answer specifically within the context of Bayesian inverse problems, not the wider statistical theory of Bayesian inference (which tends to devolve into philosophy at some point...)

Second, a general point: If you are only computing a MAP estimate and are not trying to extract higher order moments from the posterior distribution, the only (meaningful) difference between Bayesian and classical inverse problems is in the modeling, i.e., the choice of discrepancy and regularization terms (including their scaling). To put it bluntly: If you're computing a MAP estimate and you're not doing Bayesian modeling (i.e., based on objective statistical considerations), all you're doing is exchanging the labels "discrepancy term" and "regularization term" with "likelihood" and "prior", respectively.

Since you didn't give any details on where your objective comes from, I see three possibilities:

1. Your modeling is based on proper statistical considerations, i.e., you know (from somewhere) that the data $d_i$ are independently and identically normally distributed with mean $f(m)_i$ and variance $\sigma^2$, and in the absence of any further information, you'd expect the $m_j$ to also be independently and identically normally distributed with mean $m_{\text{prior},j}$ and variance $\alpha^2$. In this case the scaling you are currently observing and are uncertain about is in fact the correct scaling: Bayesian theory tells you that in this case, you need to give a very strong weight to the prior. (Similar to what Brian Borchers was telling you.)

2. Your modeling is based on improper statistical considerations. The most common instance is when your problem comes from a discretization of a parameter identification problem in a partial differential equation. In this case, your model tells you, e.g., that the unknown parameter function is normally distributed with variance $\alpha^2$, but that doesn't necessarily mean that the entries of the discrete parameter vector are. If you consider the vector $m$ as the unknown expansion coefficients of a given basis $\{\phi_i\}$ of piecewise linear functions (as in a standard finite element discretization) and the corresponding discrete parameter function $m_h := \sum_i m_i\phi_i$ as normally distributed with variance $\alpha^2$, your likelihood term would in fact be $\alpha^{-2}(m_h-m_p)^T M_h(m_h-m_p)$, where $M_h$ is the usual mass matrix. (And similarly for the likelihood.) Now, the mass matrices already account for the different scalings of both terms due to the discretization, and you can use the same parameters as given by your model. You can find discussions on the discretization of Bayesian inverse problems in http://www.siltanen-research.net/LassasSaksmanSiltanen.pdf or the works of Georg Stadler.

3. Your choice of likelihood and prior (including $\sigma$ and $\alpha$) is in fact not based on firm a priori knowledge, but chosen ad hoc. (Not that there's anything wrong with that -- people in inverse problems have been doing this for years without being ashamed of it.) Then all that's left of the Bayesian framework is the interpretation of the terms, but since there is no external objective information to relate this to, it is (from the mathematical point of view) meaningless. But this means that you are free to scale the parameters (or rather, the ratio $\sigma/\alpha$, since that is all that counts here) to your liking, without anybody being able to tell you you're wrong. Still, it is always good practice make sure your discrepancy and regularization terms are discretization-independent, so you can change the discretization without having to fiddle with the parameters again. You can then also treat the parameter as an additional unknown variable and include it in the Bayesian formulation with the help of a so-called hyper-prior.

Assuming it's either 2. or 3., the answer is therefore Yes, you can, and you should. (This will also help with the interpretation, because otherwise you mix confidence/noise covariance with a purely numerical scaling.)

Changing the relative weight of the two terms is equivalent to changing your prior.

The fundamental issue here is that you've got a problem with far more parameters to estimate than data points. You will need to use a fairly strong prior to get tight bounds on the fitted parameters, and the choice of that prior will have a strong influence on your solution. .

• So there is now way to get a tighter bound with an uninformative prior? Even in the region where $d$ is specified? – maverick Jul 17 '15 at 2:07
• Generally no. For example, in the linear inverse problem (f(m)=A*m+b) you can have an unbounded set of models that all exactly fit the data. If there's something very special about the problem that you haven't told us, then maybe something can be done. You really need to step back and consider the nature of your inverse problem first. – Brian Borchers Jul 17 '15 at 2:34

As a sidenote, the physical dimensions of your scaling parameters $\sigma,\alpha$ are different. So a statement of the form $\sigma\ll\alpha$ as given in your question doesn't make any sense: are 3kg of apples much less or much more than 5 micrometers of yarn?

In addition to the references others have given before, let me also recommend the book by Somersalo and Kaipio on "Statistical Inverse Problems". It also contains many discussions of the practicalities of inverse problems.

• Everything is dimensionless in the current problem setup. – maverick Jul 18 '15 at 2:45
• But that still doesn't mean that you can compare two numbers for size. A Reynolds number of 100 is relatively small. A Knudsen number of 10 is enormous. – Wolfgang Bangerth Jul 18 '15 at 10:28