Dealing with errors in non-linear least square problem

Question

I am currently working with a optimization problem involving a non-linear least square problem. I have chosen to use lsqnonlin in Matlab. What follows is a simplified short description of the problem.

Let $P_i=(x_i,y_i)$ $i=1,...,n$,

$A_1=(a_1,b_1)$, $A_2=(a_2,b_2)$ be a set of points which are known. Let $\tilde{R}_1$ ans $\tilde{R}_2$ be an unknown set of variables.

$$d_i = \Vert P_i-A_1\Vert-\Vert P_i-A_2\Vert-(\Vert P_i-\tilde{R}_1\Vert -\Vert P_i-\tilde{R}_2\Vert)$$

And non-linear-least square problem is then defined as:

$$\min\sum_{i=1}^{n} d_{i}^2$$

This initial setup works very well hovever I am now in need to introduce random errors on the $P_i$ variables. This makes the system very difficult to solve and my question is if any of you have any advise on how to best adapt for these errors? Is there maybe an alternative better way to pose the problem?

Let me know if anything is unclear.

Answer 1

Your problem can be described as a regression problem where not only the function values $y_i$ but also the points $x_i$ where the values are given are noisy. The classical approach to solve such a problem is called total least squares, which basically amounts to fitting the pairs $(x_i,y_i)$ using regular least squares (in a higher-dimensional space). The classical reference is Golub, van Loan: An analysis of the total least squares problem. SIAM J Numerical Analysis 17 (1980), 883-893 (Charles van Loan's website has a -- rather poorly, but freely available -- scanned version).

For a simple regression problem, total least squares can be formulated as minimizing $$\sum_{i=1}^n\frac{d_i^2}{1+\|x\|_2^2}.$$

In general, total least squares applies to finding a vector $x\in\mathbb{R^n}$ satisfying $Ax=b$ for a given matrix $A\in\mathbb{R}^{m\times n}$ and vector $b\in \mathbb{R}^m$ such that $Ax=b$. (Here, $x$ could contain the coefficients of a regression polynomial, $A$ the polynomial basis evaluated in the given regression points and $b$ the given values.) If both $A$ and $b$ are contain errors, you cannot have $Ax=b$, but there are some (unknown) $E$ and $e$ such that $(A+E)x=(b+e)$. Total least squares then amounts to solving $$\min_{E,e} \|E\|_F^2 +\|b\|_2^2 \quad\text{subject to}\quad (A+E)x=(b+e),$$ where $\|E\|_F$ is the Frobenius norm. This is a standard least squares problem for the extended matrix $\hat E:=[E|e]\in \mathbb{R}^{(n+1)\times m}$, which can be solved (if it can be solved at all; this is not always the case) using the singular value decomposition of $\hat E$: If $v_{n+1}\in\mathbb{R}^{n+1}$ is a singular vector corresponding to the smallest singular value $\sigma_{n+1}$ with $v_{n+1} = [y^T|\alpha]^T$ for some $y\in\mathbb{R}^n$ and $\alpha\in\mathbb{R}$, then the solution is given by $$ x = -\frac{1}{\alpha} y.$$

Answer 2

You could consider trying to solve this problem using some more global optimization approach, such as a particle swarm optimization.

Given that you set this up sufficiently well, this approach may be able to traverse the solution field better than a gradient based approach since you may be able to avoid local minima better.

Answer 3

A possible way of solving this is with Bayesian inference. Since you have already a model, you need to define the priors to $R_1$ and $R_2$ as well as the distribution for the error in $P_i$. For this approach, fhe following steps have to be performed

1. Set up a full probability model defining the joint probability distribution for all observable ($P_i$) as well as unobservable quantities ($R_1, R_2$).

2. Define the conditioning on the observed data, calculate the posterior distribution given the observed data. $p(R_1,R_2|P_i)\propto p(P_i|R_1,R_2)p(R_1,R_2)$.

3. Estimate the unobserved quantities $(R_1,R_2)$ of interest given the observed data $(P_i)$ from the posterior distribution. Estimate their confidence intervals.