Why does fmincon using the sqp algorithm need a full matrix to specify linear constraints?

I need to specify that some of the variables in my optimization need to be ordered, so v_i-v_(i+1)<0 for i = N_1 to N_2. This is easy to set up in fmincon, in the form Av<b, but, as you can see, the algorithm can't take a sparse input for A. In a large problem like the one I'm solving, this gives me about a 1000 by 4000 matrix that has only 2000 nonzero entries, which seems crazy to me.

I know I could reformulate the problem so that I use the differences in variables instead (and I will probably have to end up doing this), which gets around the problem but (i) this makes the calculation of the gradient more awkward and (ii) I can see no good reasonnot to allow a sparse A.

So, my question is, is there a good reason for this 'feature'?

Answer 1

[Matt J]

Probably because the SQP algorithm must extract subsets of rows of the constraint matrice quite often. This can be much faster for full matrices than sparse (see below). Also, even if your constraints are sparse, there is no reason to expect that the quadratic programming sub-problems of the algorithm will be. So it would be of questionable benefit to have constraints alone in sparse form.

A=sprand(4000,4000,4000/4e5)+speye(4000);
Af=full(A);


tic
for i=1:size(A,1)
 a=A(i,:);
end
toc

Elapsed time is 0.676821 seconds.

 
tic
for i=1:size(A,1)
 a=Af(i,:);
end
toc

Elapsed time is 0.155082 seconds.

Finally, even if everything in your problem is sparse, any matrix inversions that must be done will often be faster in full form than sparse for small scale problems, which sqp problems are expected to be. Example,

b=rand(4000,1);
tic
 A\b;
toc

Elapsed time is 1.480658 seconds.

 
tic
 Af\b;
toc

Elapsed time is 0.507184 seconds.

[John D'Errico]

I think I agree with Matt. But if so, that would make the problem difficult to resolve by any means other than working in a difference formulation for the unknowns.

[John Billingham]

Hmmmm. OK. Thanks. I suppose it depends on whether I'm running out of time or running out of memory (neither as yet). It would however be nice to have to option to make it sparse.

[Matt J]

If you are running out of memory, you shouldn't be using SQP. You should be using one of the large scale algorithms.

[John Billingham]

Probably, but none of them works as well as sqp on the problem I have to deal with.

Answer 2

[John D'Errico]

This sounds like a good excuse for a feature request to me. However, my gut tells me there is a technical reason in the algorithm itself (possibly in terms of how it was implemented) that precludes a sparse matrix for the linear constraints. And, unfortunately, we cannot look into the code itself, since sqpInterface is p-coded. If I had to make a wild guess, the issue comes down to the fact that QR (used to) behave subtly differently for sparse versus full matrices. I recall that in the past, you could not use QR to get a pivoted economy QR factorization when the matrix was sparse. This is not the case today. So this may possibly be only an issue that was true for past releases, and need not be true today, but someone would need to dive into the fmincon code and make the necessry changes.

Anyway, the best way to get a good answer for this would need to come from TMW directly. And the best way to get that is to contact them directly.