# Dealing with errors in non-linear least square problem

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.

• Probably define the partial derivatives of the objective function first, followed by gradient descent based on the Jacobian. You provided essentially the sum of squares $d^2$, and stated that the SS needs to be minimized, however, what is the function you are fitting? – user15259 May 3 '15 at 21:08
• I think what you are looking for is total least squares. – Christian Clason May 6 '15 at 12:04
• Yes thank you @ChristianClason. Total least squares is performing significantly better, thank you. – Seb May 16 '15 at 12:10
• Excellent, glad to hear that! (If it solved your problem and you think it would be useful, I could write a short answer summarizing the approach.) – Christian Clason May 16 '15 at 12:24
• Yes please do if you want to @ChristianClason. I think that it would be useful. – Seb May 16 '15 at 14:12

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.$$

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.

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.