Least absolute deviations solving using the Barrodale-Roberts-algorithm: Premature termination?

Question

Please excuse the longish question, it just needs some explanation to get down to the actual problem. Those familiar with the mentioned algorithms probably could jump directly to the first simplex tablau.

---

To solve least absolute deviation problems (a.k.a $L_1$-optimization), the Barrodale-Roberts-algorithm is a special purpose simplex method that needs far less storage and computational efforts to find a suitable minimum.

My implementation of the algorithm terminates at a simple example before a proper minimum is reached. However, probably let me state the problem in a more elaborated manner first:

Given data $(x_i,y_i)$, $L_1$-optimization tries to find $c\in m$ that minimizes $$ \sum_{i=1}^n |y_i-f(x_i)| \quad\text{with}\quad f(x):=A_x\cdot \phi $$ where $A_x$ is a $n\times m$ matrix that depends in some way on $x$. This problem can be stated as a linear program and therefore among others be solved using simplex-like methods.

Barrodale and Roberts suggested an (apparently widely used) modification of the simplex method that radically simplifies the simplex method using the special structure of $L_1$-problems. Most notably, this is that an optimal solution interpolates at least $\mathop{rank}(A)$ of the given datapoints. Those with Jstor access may find the corresponding article here.

Lei and Anderson in 2002 proposed a small modification that is supposed to increase numerical stability and therefore to overcome known problems with the simplex algorithm.

Basically, this algorithm assumes that you start with a given set of points that have to be interpolated, use the given procedures to build a simplex tableau and then use the rules of Barrodale and Roberts to decide on which basis variables to change and therefore modify the set of datapoints that is approximated.

Barrodale and Roberts give a small example that I tried to reproduce. It tries to approximate the points $\{(1,1), (2,1), (3,2), (4,3), (5,2)\}$ by a function $a_1+a_2x$. The finish their algorithm with the following condensed simplex tableau:

$$ \begin{array}{l|l|ll} \hline \text{Basis}&R& u_1 & u_3\\ \hline b_1 & 1/2 & 3/2 & -1/2 \\ v_2 & 1/2 & 1/2 & 1/2 \\ b_2 & 1/2 & -1/2 & 1/2 \\ u_4 & 1/2 & 1/2 & -3/2 \\ v_5 & 1 & -1 & 2 \\ \hline \text{Marginal cost} &2 & -1 & 0\\ \hline \end{array} $$

Most importantly, the first and third point are interpolated and the overall error is equal to $2$. They conclude that

Since all of the nonbasic vectors have nonpositive marginal cost [...]

the iteration is finished and and optimum is reached.

If I use the algorithm of Lei and Anderson, I can reproduce that simplex tableau for the interpolation set {1,3}, as it is expected. However, if I start the algorithm with the set $\{2, 5\}$ (which clearly is not optimal), I get the following simplex tableau:

\begin{array}{l|l|ll} \hline \text{Basis}&R& u_2 & u_5\\ \hline u_1 & 1/3 & -4/3 & 1/3 \\ b_1 & 1/3 & 5/3 & -2/3 \\ u_3 & 2/3 & -2/3 & -1/3 \\ u_4 & 4/3 & -1/3 & -2/3 \\ b_2 & 1/3 & -1/3 & 1/3 \\ \hline \text{Marginal cost} &7/3 & -10/3 & -5/3\\ \hline \end{array}

This result is puzzling me, though. If I understand the quote above correctly, having no positive marginal cost indicates that the optimum is reached. The function value of about 2.33 certainly is not optimal, though. Exchanging $u_2$ with $u_1$ would yield a result that is on par with the solution of Barrodale and Roberts and therefore optimal.

Additional info: If I start with the initial tableau given by Barrodale and Roberts, I am also able to reproduce the tableau above by ordinary simplex steps, thus I am fairly confident that the actual numerical values are correct and my interpretation of the pivot selection rule is faulty.

Any thoughts on this?

I realize that the question itself is quite complicated and probably requires knowledge of at least the Barrodale and Roberts algorithm to be answered sufficiently. The algorithm as a whole is to long to repeat it here in full detail. However, if you have additional questions on steps I took or on missing bits of information, feel free to ask and I will gladly augment the question.

Answer 1

Solved it. Actually, Barrodale and Roberts solved it and I just did not read carefully.

In my question I left it to the reader to understand that the variables Barrodale and Roberts labeled $u_i$ stand for the positive residuals the $i$-th datapoint in relation to the current fit. If the residual is negative, $u_i=0$ and $v_i$ takes the corresponding value. As only one of them may be within the basis and the coefficients in the simplex tableau are just the negatives of each other, it is not necessary to explicitly state them in the simplex tableau. Barrodale and Roberts mention in their article:

[...] and that the sum of the marginal (or reduced) costs of $b_j$ and $c_j$ is zero and of $u_i$ and $v_i$ is -2.

Thus, my simplex tableau above has to be thought to look as follows:

\begin{array}{l|l|ll| ll} \hline \text{Basis}&R& u_2 & u_5 & v_2 & v_5\\ \hline u_1 & 1/3 & -4/3 & 1/3 & 4/3 & -1/3 \\ b_1 & 1/3 & 5/3 & -2/3 & -5/3 & 2/3 \\ u_3 & 2/3 & -2/3 & -1/3 & 2/3 & 1/3 \\ u_4 & 4/3 & -1/3 & -2/3 & 1/3 & 2/3 \\ b_2 & 1/3 & -1/3 & 1/3 & -1/3 & -1/3 \\ \hline \text{Marginal cost} &7/3 & -10/3 & -5/3 & 4/3 & -1/3\\ \hline \end{array}

where we clearly see that $v_2$ can be brought into the basis to archive a better result. Doing this the algorithm terminates while interpolating the first and fifth datapoint with an overall error of 2 - which is the best solution.

Thanks for reading and giving me some place to write my problem down, which usually helps to narrow the solution down significantly. Hopefully, this answer might sometimes be helpful for anybody else trying to implement Barrodale & Roberts.