Change in linprog.m output?

I use the following code to fit multiple data sets using linprog.m.

for k=1:nsets
      b = [data(:,k)+2*noise; -data(:,k)+2*noise];
      [fitpars(:,k), ~, flag(k)] = linprog(w,A,b,[],[],lb,ub,[],opts);
end

A solution is not found for each data set, so I save the exitflag. The code executes successfully in R2016b, but in R2018a I get the error message "Unable to perform the assignment because the size on the left side is 320-by-1 and the size on the right side is 0-by-0." whenever the exitflag does not equal 1. (fitpars(:,k) is size 320-by-1) In both R2016b and R2018a I use the default algorithms. Is there a way to get R2018a to provide output compatible with the output from R2016b?

Answer 1

[Matt J]

You should save to a cell array.

 [fitpars{k}, ~, flag(k)] = linprog(w,A,b,[],[],lb,ub,[],opts);

It was just random luck that your original code ever worked. You cannot stick an empty result into the column of a non-empty matrix.

[Walter Roberson]

I have seen in the past that sometimes when the optimization routines terminated, they returned whatever the current objective point happened to be, even if the routine had identified that the point was not a valid solution.

In some cases that sort of makes sense: if, for example, you ended up in a local minima, you would probably want to know where the minima was.

[Matt J]

I can see that. But in the case of linprog, I don't think you can trust that a non-empty output will always be returned, even in older versions of Matlab. It has always been possible that a "current objective point" may never be generated. The linprog code does certain pre-analysis of the problem's feasibility and, if it pre-determines that the problem is infeasible, will abort right away without ever attempting the search.

Possibly, R2018a is better than older versions at detecting these cases and is doing early aborts more often than R2016. This might explain the OP's experiences.

[Walter Roberson]

I would have expected that the early aborts would error(), giving something to try/catch on ?

[Matt J]

Not currently, but you could go back to older Matlab versions and test how situations like the following were handled,

 [x,fval,exitflag,output] = linprog(1,[1;-1],[1;-2]); 
 x,iters=output.iterations

      x =

         []

    iters =

         0

Personally, I don't know if I would like for linprog to treat an infeasible problem as an error.

[Joel Miller]

Thanks to Matt for the workaround; however, it is not random luck that my code works in previous versions of Matlab. In R2016b, at least when exitflag = -2 (no feasible point was found), linprog.m outputs a numerical array of the proper dimension. In fact, there is a bug report https://www.mathworks.com/support/bugreports/698379?s_tid=srchtitle from 2011 saying that the bug "linprog reports a matrix dimension mismatch error on some infeasible problems." was fixed in R2011a. It appears the bug has returned.

[Matt J]

In R2016b, at least when exitflag = -2 (no feasible point was found), linprog.m outputs a numerical array of the proper dimension.

Interesting. What is the output of my example code above in R2016b, just out of curiosity?

As I said to Walter, though, I don't know why one would prefer a non-empty output in cases when we know the problem is infeasible. It's unnecessary extra memory consumption and also an extra hazard for users who are not careful enough to check the exitflag.

[Joel Miller]

This is what I get when I run your code in R2016B.

[x,fval,exitflag,output] = linprog(1,[1;-1],[1;-2]); 
 x,iters=output.iterations
Warning: Your current settings will run a different algorithm ('dual-simplex') in a future release.

> In linprog (line 204) 
Exiting: One or more of the residuals, duality gap, or total relative error
 has stalled:
         the primal appears to be infeasible (and the dual unbounded).
         (The dual residual < OptimalityTolerance=1.00e-08.)

x =

      1.0000

iters =

      15

I agree that it doesn't seem sensible to provide numerical output that is not a solution to the problem, but since that was the behavior, the two lines of code in my original question handled the output simply. I am surprised I have not been able to find any documentation on the change in output.

[Matt J]

In fact, there is a bug report ... "linprog reports a matrix dimension mismatch error on some infeasible problems." was fixed in R2011a

I don't think it is the same issue. In your case, it is not linprog that is throwing the dimension mismatch error. That error happens in the workspace that calls linprog.

This is what I get when I run your code in R2016B.

What about something even sillier, like

 [x,fval,exitflag,output] = linprog(1,[],[],[],[],2,1); 

x, iters=output.iterations

Do you still get non-empty output, x?

[Walter Roberson]

When I test that in R2016b, I get fval = [] and x 0×1 empty double column vector

[Joel Miller]

Same here, and iters = 0.

[Matt J]

OK, well I think that proves my earlier conjecture. There was never any guarantee, even in R2016b, that the output of linprog will be non-empty, so storing results to a cell array was always the only safe thing to do.