My goal is to robustly estimate a general 2D line from $n$ data points, where the line is parameterized by $\rho > 0$, the distance from the origin to the line and $\varphi$, the angle between the $x$-axis and the perpendicular line (that goes through the origin). This parametrization has been chosen as it requires two parameters only and can also represent vertical lines.

To robustly estimate the two parameters $\rho > 0$ and $\varphi$, I can make use of a general robust linear least-squares solver (RLLS), which solves $\min_\mathbf{c} \vert\vert \mathbf{X}\,\mathbf{c}-\mathbf{y}\vert\vert_2$, where $\mathbf{X}$ is the $n\times p$ regressor matrix and $\mathbf{c}$ is the $p$-dimensional parameter vector.

Is there an elegant way (a special transformation that allows to express the estimation of $\varphi$ and $\rho$ in terms of a linear regression setting) to make use of the RLLS that can handle arbitrary (also vertical) lines?


I figured out, that the distance between an arbitrary data point $\left(x_i, y_i\right)$ from a line with $\rho$ and $\varphi$ is given as

$$ x_i\,\cos\left(\varphi\right) + y_i\,\sin\left(\varphi\right) - \rho $$

To get this into the format the RLS solver takes as input, I rearranged the equation to

$$ \begin{bmatrix}\dots&\dots\\x_i&1\\\dots&\dots\end{bmatrix}\begin{bmatrix}-\frac{\cos\left(\varphi\right)}{\sin\left(\varphi\right)}\\\frac{\rho}{\sin\left(\varphi\right)}\end{bmatrix} = \begin{bmatrix}\dots\\y_i\\\dots\end{bmatrix}, $$ The parameters $\mathbf{c} = \left[c_1, c_2\right]^\top$, were then used to compute $\varphi = \mathrm{atan}\left(-\frac{1}{c_1}\right) $ and $\rho = c_2\,\sin\left(\varphi\right)$.

This works but requires some hacks: The solver may return values of $\rho < 0$, which can be handled by taking its negative and shifting $\varphi$ by $\pi$.

However, the greatest issue are straight lines, where $\varphi \approx 0$. To handle horizontal lines I do the following: I solve the problem twice, normally and with swapped $x$ and $y$ axis. The solution with smaller residual means squared error is accepted. This is not elegant and requires solving the least squares problem twice.

enter image description here

enter image description here


Suppose we have a collection of points in the plane and we want to draw a straight line through these points, in such a way that the line is a "best fit".
An equation for a straight line can always be set up as follows: $$ \cos(\theta) (x - p) + \sin(\theta) (y - q) = 0 $$ Here $\theta$ is the angle of the line's normal with the x-axis and $(p,q)$ is a point supporting the line. The distance of an arbitrary point $\vec{r}_k = (x_k,y_k)$ to the line is given by the length of the projection of the point's vector $\vec{r}_k$ onto the line. The normal $\vec{n}$ of the line is given by $\vec{n} = (\cos(\theta),\sin(\theta))$. Hence the length of the projection is: $$ \left| \frac{(\vec{r}_k \cdot \vec{n})}{(\vec{n} \cdot \vec{n})} \vec{n} \right| = \left| \cos(\theta) (x_k - p) + \sin(\theta) (y_k - q) \right| $$ For the straight line to be a "best fit", it will be required that the sum of the weighted squares of all distances to the line shall be a minimum: $$ \sum_k w_k \left[ \cos(\theta) (x_k - p) + \sin(\theta) (y_k - q) \right]^2 = \mbox{minimum}(p,q,\theta) $$ Here the weights $w_k$ are chosen in such a way that $\sum_k w_k = 1$ . Working out a bit: $$ \cos^2(\theta) \sum_k w_k (x_k - p)^2 + \sin^2(\theta) \sum_k w_k (y_k - q)^2 + $$ $$ 2 \sin(\theta) \cos(\theta) \sum_k w_k (x_k - p) (y_k - q) = \mbox{minimum}(p,q,\theta) $$ Let us solve just one part of the puzzle, namely: how the points $(p,q)$ must be selected in such a way that a minimum may be reached with respect to this choice. For certain parts of the above expression this would mean that: $$ \sum_{k} w_k (x_k - p)^2 = \mbox{minimum} \quad \mbox{and} \quad \sum_{k} w_k (y_k - q)^2 = \mbox{minimum} $$ It is sufficient to consider only the expression in $x$. Differentiate to $p$: $$ \frac{d}{dp} \left[\sum_k w_k x_k^2 - 2 p \sum_k w_k x_k + p^2 \sum_k w_k \right] = - 2 \sum_k w_k x_k + 2 p \sum_k w_k = 0 $$ $$ \quad \Longrightarrow \quad p = \mu_x = \sum_k w_k x_k $$ And in very much the same way: $$ q = \mu_y = \sum_k w_k y_k $$ Define second order momenta with respect to the midpoint as usual: $$ \sigma_{xx} = \sum_k w_k (x_k - \mu_x)^2 \quad \mbox{and} \quad \sigma_{yy} = \sum_k w_k (y_k - \mu_y)^2 $$ $$ \sigma_{xy} = \sum_k w_k (x_k - \mu_x) (y_k - \mu_y) $$ We can concentrate now on minimalization with respect to the angle $\theta$: $$ \cos^2(\theta) \sigma_{xx} + \sin^2(\theta) \sigma_{yy} + 2 \sin(\theta) \cos(\theta) \sigma_{xy} = \mbox{minimum}(\theta) $$ Extreme values may be found by differentiation to the independent variable: $$ - 2 \sin(\theta) \cos(\theta) \sigma_{xx} + 2 \cos(\theta) \sin(\theta) \sigma_{yy} + 2 \cos^2(\theta) \sigma_{xy} - 2 \sin^2(\theta) \sigma_{xy} = 0 $$ Which leads to the equation: $$ - \sin(2 \theta) (\sigma_{xx} - \sigma_{yy}) + \cos(2 \theta) 2 \sigma_{xy} = 0 $$ A solution is: $$ \tan(2\theta) = \frac{2 \sigma_{xy}}{\sigma_{xx} - \sigma_{yy}} \quad \Longrightarrow \quad \theta = \frac{1}{2} \left[ \arctan\left(\frac{2 \sigma_{xy}}{\sigma_{xx} - \sigma_{yy}}\right) + k\cdot\pi \right] $$ with $k = 0,1,2, \cdots$ . But the extreme value can still be a minimum or a maximum. In order to decide about this, we have to evaluate the expression: $$ M(\theta) = \cos^2(\theta) \sigma_{xx} + \sin^2(\theta) \sigma_{yy} + 2 \sin(\theta)\cos(\theta) \sigma_{xy} $$ for two different values of $\,\theta$ , namely with $\,k=0,1$ .

  • $\begingroup$ Thank you very much for your answer! I will give it a try. $\endgroup$
    – Flo Ryan
    Sep 15 '21 at 14:52
  • $\begingroup$ For the sake of completeness, a special case should be distinguished, namely $\sigma_{xx}=\sigma_{yy}$, making the denominator of the $\arctan$ argument zero. In that case all lines through $(p,q)$ are a best fit. $\endgroup$ Sep 16 '21 at 9:08
  • 1
    $\begingroup$ Also it should be noted that the answer doesn't address the entire problem. A robust estimate requires a reasonable choice for the weights. To do so, I computed an initial solution with uniform weights. After this, I computed the distance of every point to the current best fit. Based on the magnitude of that distance, I recomputed the weights (using Huber / Bisquare weights) and recomputed the the best fit with the new weights. This weigh -> compute -> weigh ... procedure was repeated until convergence. Just like it's done in IRLS. $\endgroup$
    – Flo Ryan
    Sep 16 '21 at 9:37

The equation of the best fit line is

$ (\cos \phi, \sin \phi) \cdot (r_i - r_0) = E_i \hspace{24pt} (1)$

Where $r_0 = (x_0, y_0) $ is a fixed point on the best fit line, and $r_i = (x_i, y_i)$, and $E_i$ is the equation error, introduced because points $r_i$ will not, in general, fit the equation perfectly.

Now we want to minimize of the sum of squared equation errors. That is we want to minimize $f(\phi, r_0)$ given by

$ f(\phi, r_0) = \displaystyle \sum_{i=1}^n E_i^2 = E^T E$

where the $i$-th entry of vector $E$ is $E_i$.

At the minimum we have the partial derivatives with respect to $\phi$ and to $r_0$ equal to zero.

Taking the partial derivative of $f(\phi, r_0) $ with respect to the vector $r_0$, we get,

$ \dfrac{\partial f} {\partial r_0} =\displaystyle -2 [\cos \phi, \sin \phi ]^T \sum_{i=1}^{n} E_i$

Therefore, by equating $\displaystyle \sum_{i=1}^n E_i $ to $0$, we have,

$\begin{equation} \begin{split} 0 &= \displaystyle \sum_{i=1}^n (\cos \phi, \sin \phi) \cdot (r_i - r_0) =(\cos \phi, \sin \phi) \cdot \sum_{i=1}^n (r_i - r_0) \\ & = (\cos \phi, \sin \phi) \cdot ((\sum_{i=1}^n r_i) - n r_0) \end{split} \end{equation}$

Hence, it will be sufficient to take

$ r_0 = \displaystyle \frac{1}{n} \sum_{i=1}^n r_i $

So, now, we have a point on the best fit line. What is left is to determine $\phi$. Equation $(1)$ can be written in matrix-vector form as

$ A x = E $

where $ x = [\cos \phi, \sin \phi]^T$, and the $i$-th row of matrix $A$ is $(r_i - r_0)^T$.

Here, again, we want vector $E$ to have the smallest norm possible, which is the same as requiring the minimum sum of squared errors.

$E^T E = x^T A^T A x $, and the unknown vector $x$ is a unit vector. It is well-known that the minimum of $x^T A^T A x $ is $\lambda_{\text{min}}(A^T A) $. This means that all we have to do is compute the $2 \times 2 $ matrix $A^T A$ , determine its two eigenvalues and corresponding eigenvectors. Then the required vector $x$ is just the (unit) eigenvector corresponding to the smaller eigenvalue, and our job is done.

  • $\begingroup$ Thank you for your answer! I tried something very similar initailly. But due to the presence and effect of outliers, I was unable to follow this track. Since there are very established methods for Robust Least Squares (e.g. IRLS) I tried to get it into a Ax=b shape. I tried your method, computing a robust mean $r_0$ and using IRLS on $\mathbf{Ax}=\mathbf{0}$ (Internally it does a SVD and I use the last column of $V$ as a non-trivial solution to $\mathbf{Ax}=\mathbf{0}$). This works but I am not all too pleased. $\endgroup$
    – Flo Ryan
    Sep 10 '21 at 8:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.