Why is the 'active-set' algorithm being removed from LSQLIN when the alternative 'interior-point' algorithm lacks the same performance?

I'm providing an example of a problem that is typical in our use of LSQLIN. We solve a physics based problem where all the constraints are required to ensure a physically real solution. Solving problems like this are at the core of a much larger problem, where LSQLIN is called ~10,000 times per problem. As the LSQLIN ‘interior_point’ algorithm doubles computational time, the removal of ‘active-set’ makes our data processing speed slow down significantly. I understand that 'active-set' may not be available in the future. Does Mathworks plan on addressing what seems to be a downgrade of the LSQLIN to an inferior state? Will ‘interior-point’ ever offer the same performance as ‘active-set’?

Variables C, d, A and b can be loaded from the attached MAT-file. C: is full and 45x6 d: 45x1 A: 9x6 each row has 3 non-zero values. b: 9x1 zeros array

load example
tic; x1 = lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'active-set','Display','off')); toc
tic; x2 = lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'interior-point','Display','off')); toc
Elapsed time is 0.002902 seconds.
Elapsed time is 0.006580 seconds.

[Titus Edelhofer]

Hi David,

which version are you using? I tried with R2017a and the difference was not that large:

x1 = @() lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'active-set','Display','off'));
x2 = @() lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'interior-point','Display','off'));
t1 = timeit(x1)
t2 = timeit(x2)

t1: 0.0039 t2: 0.0049

So there is a difference on my machine, but not as significant ...

Titus

[David Sinex]

Hi Titus,

I was using R2016b. I have rerun a more comprehensive test using R2017a that better relates to our overall problem where LSQLIN is run many 1000s of times. I suspect in your case the small test you ran was a slight outlier. Try the code snippet below and see if you get a plot similar to mine.

for kk = 1:10000
    tic; x1 = lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'active-set')); t1(kk) = toc;
    tic; x2 = lsqlin(C, d, A, b, [], [], [], [], [], optimoptions(@lsqlin,'Algorithm', 'interior-point')); t2(kk) = toc;
end

rat = t1./t2;

figure(301)
clf;
ax = subplot(3,1,1);
histogram(t1)
ax.XLim = [0.001 0.005];
title('active-set')
xlabel('Time (s)')

ax = subplot(3,1,2);
histogram(t2)
ax.XLim = [0.001 0.005];
title('interior-point')
xlabel('Time (s)')

ax = subplot(3,1,3);
histogram(rat)
ax.XLim = [0.3 0.65];
title('t1 / t2')
xlabel('Time (s)')

Answer 1

[Sailesh Sidhwani]

Interior-point Algorithm:

'interior-point' handles large, sparse problems, as well as small dense problems. The algorithm satisfies bounds at all iterations, and can recover from NaN or Inf results. It is a large-scale algorithm. The algorithm can use special techniques for large-scale problems.

Active-set Algorithm:

'active-set' can take large steps, which adds speed. The algorithm is effective on some problems with nonsmooth constraints. It is not a large-scale algorithm.

Large Scale vs Medium Scale Algorithms:

An optimization algorithm is large scale when it uses linear algebra that does not need to store, nor operate on, full matrices. This may be done internally by storing sparse matrices, and by using sparse linear algebra for computations whenever possible. Furthermore, the internal algorithms either preserve sparsity, such as a sparse Cholesky decomposition, or do not generate matrices, such as a conjugate gradient method.

In contrast, medium-scale methods internally create full matrices and use dense linear algebra. If a problem is sufficiently large, full matrices take up a significant amount of memory, and the dense linear algebra may require a long time to execute.

Please refer to the "Reasoning Behind the Recommendations" section on the link below:

https://www.mathworks.com/help/optim/ug/choosing-the-algorithm.html

[Walter Roberson]

David Sinex comments,

This does not address the question. This is simply a copy of a segment of the documentation with no further explanation offered by the author. I agree the content relates to the material in my question, but thats all. I would appreciate a thoughtful answer.

[Walter Roberson]

David: the volunteers do not know Mathworks' plans about this -- or if we do know, then we are not permitted to talk about it as we would have only found out under Non-Disclosure Agreement.

You should be creating a support case.

Answer 2

[C.T. Kelley]

This is a problem for me. I have an ill-conditioned problem and the active set method does not use the normal equations as part of the KKT system. That is a big deal for me.

Is there a public-domain implementation of the Golub-Saunders 1969 algorithm in Matlab out there?

[John D'Errico]

Please don't resurrect an old question with a question of your own, and don't pose it as an answer.