@alex_01

You wrote, "later comment-2... how do you solve Min(Norm(R.W)) with LP?", and incidentally you were using the 2-norm.

Well, that is not a linear objective so it isn't a Linear Programming problem, so `LPSolve` won't solve that formulation. I showed that your posted problem can be taken as a Quadratic Programming problem with a quadratic objective and linear constraints, and you agreed. What now are your grounds for thinking that it can also be reformulated as an equivalent LP problem!?

You also wrote, "Unfortunately, I already know that the QP and the NLP optimization produces the same allocations..." which makes it sound as if you are unhappy with the computed solution, and are hoping for a better one. But this (according to evidence from you) is a convex QP problem and the objective is bounded from below in the feasibility region, so the computed solution should be globally optimal. How then can you hope for a better solution!?

The only way that I can currently envision getting a more desirable solution is if you have additional qualifying restrictions or constraints in mind (but as yet not stated in this thread). For example, perhaps you want an integer-valued solution (well, binary valued, for these ranges), which is something that I just make up as an example. Do you have additional constraints?

So, you want to construct the Matrix form for the NLP formulation of this QP class problem. Why is that? Is it so that you can form procedures for objective, constraints, objective gradient, and constraint jacobian which all accept purely numeric vector arguments, compute results and fill some arguments with the results inplace, and that is evalhf'able or ideally even compilable? And you want them highly optimized and tuned to the quadratic & linear natures of the relevant formulas? It can be done (though keeping down overhead from callbacks through main Maple to get from external Optimization wrappers to Compiled procs is awkward).

But you know who's already done that, without any bypass overhead? The Numerical Algorithms Group. They wrote all the specialized QP objective, linear constraint evaluation, and gradient/jacobian/hessian stuff and compiled it, together. And Maple's Optimization:-QPSolve passes the float[8] data of the QP Matrix form to NAG's e04nfa via Optimization:-External:-E04NFA.

Are your problems so large and sparse that e04nfa is not good enough? How much specialized optimization to the procs and code do you think you can muster, to make up for what appears to be a speed difference factor of about 7 between QPSolve and NLPSolve used plainly.

Look, I don't want to be rude, but it's a little frustrating that you still seem so very disinclined to believe my advice. You wanted to do QP, so I set up the Matrix form and explained it. You wanted to do LP, so I set up the Matrix form and explained it. Why not believe what I suggest about this route into NLP?

I'm also still confused about how much it seems to matter to you whether a 7 second (nstock=200) code section is 100% doubleplus-interior-point-semidefinite-cone-fantastic when the preamble before it is running 80 seconds instead of 30 or so. But then, you have to date left unanswered queries as to how large your typical problems will be.

@alex_01

You wrote, "later comment-2... how do you solve Min(Norm(R.W)) with LP?", and incidentally you were using the 2-norm.

Well, that is not a linear objective so it isn't a Linear Programming problem, so `LPSolve` won't solve that formulation. I showed that your posted problem can be taken as a Quadratic Programming problem with a quadratic objective and linear constraints, and you agreed. What now are your grounds for thinking that it can also be reformulated as an equivalent LP problem!?

You also wrote, "Unfortunately, I already know that the QP and the NLP optimization produces the same allocations..." which makes it sound as if you are unhappy with the computed solution, and are hoping for a better one. But this (according to evidence from you) is a convex QP problem and the objective is bounded from below in the feasibility region, so the computed solution should be globally optimal. How then can you hope for a better solution!?

The only way that I can currently envision getting a more desirable solution is if you have additional qualifying restrictions or constraints in mind (but as yet not stated in this thread). For example, perhaps you want an integer-valued solution (well, binary valued, for these ranges), which is something that I just make up as an example. Do you have additional constraints?

So, you want to construct the Matrix form for the NLP formulation of this QP class problem. Why is that? Is it so that you can form procedures for objective, constraints, objective gradient, and constraint jacobian which all accept purely numeric vector arguments, compute results and fill some arguments with the results inplace, and that is evalhf'able or ideally even compilable? And you want them highly optimized and tuned to the quadratic & linear natures of the relevant formulas? It can be done (though keeping down overhead from callbacks through main Maple to get from external Optimization wrappers to Compiled procs is awkward).

But you know who's already done that, without any bypass overhead? The Numerical Algorithms Group. They wrote all the specialized QP objective, linear constraint evaluation, and gradient/jacobian/hessian stuff and compiled it, together. And Maple's Optimization:-QPSolve passes the float[8] data of the QP Matrix form to NAG's e04nfa via Optimization:-External:-E04NFA.

Are your problems so large and sparse that e04nfa is not good enough? How much specialized optimization to the procs and code do you think you can muster, to make up for what appears to be a speed difference factor of about 7 between QPSolve and NLPSolve used plainly.

Look, I don't want to be rude, but it's a little frustrating that you still seem so very disinclined to believe my advice. You wanted to do QP, so I set up the Matrix form and explained it. You wanted to do LP, so I set up the Matrix form and explained it. Why not believe what I suggest about this route into NLP?

I'm also still confused about how much it seems to matter to you whether a 7 second (nstock=200) code section is 100% doubleplus-interior-point-semidefinite-cone-fantastic when the preamble before it is running 80 seconds instead of 30 or so. But then, you have to date left unanswered queries as to how large your typical problems will be.

@pet1 If you reread that last post above by Doug, you should be able to see that it was a wish: a description of extended functionality for dsolve/numeric with parameters which Doug would like to see implemented in a future version of Maple.

Part of the point was that it does not currently work like that (ie. so simply) at present (Maple 15).

If you look at a post I did at the beginning of that month of May, 2010, you can see a slightly different approach which can give something pretty close to what Doug describes, but of course it too needs a little setup. Doug expressed a wish for dsolve/numeric to handle to it all more nicely so that setup (whether his way, or mine, or other) can be sidestepped.

There's not much different between our approaches. I used `subs` on a proc because I like to avoid global variables, and I didn't bother to have it return unevaluated on nonnumeric input (hence my uneval quotes when calling it below).

Eg,

restart:
with( plots ):

deqs := { diff(x(t),t) = h*sin(t) + y(t),
          diff(y(t),t) = x(t),
          x(0)=0, y(0)=0}:

integ := dsolve(deqs, numeric, parameters=[h]):

master:=subs(dummy = eval(integ),
             proc(h) dummy('parameters' = [h]);
             dummy end proc):

# A bit like Doug described...
animate(odeplot, ['master'(H), [[t, x(t)], [t, y(t)]], t = 0 .. 3], H = -1 .. 5);

# ...the way I happened to use `master`
display([seq(odeplot(master(H), [[t, x(t)], [t, y(t)]], 0..3,
             'numpoints' = 500, 'thickness' = 2),H=-1..5)],
        'axes' = 'box', 'insequence' = true);

The above is about Doug's `animate` call. Something similar to Doug's `plot3d` call above could also be made to work, also with a little (different) setup.

@pet1 If you reread that last post above by Doug, you should be able to see that it was a wish: a description of extended functionality for dsolve/numeric with parameters which Doug would like to see implemented in a future version of Maple.

Part of the point was that it does not currently work like that (ie. so simply) at present (Maple 15).

If you look at a post I did at the beginning of that month of May, 2010, you can see a slightly different approach which can give something pretty close to what Doug describes, but of course it too needs a little setup. Doug expressed a wish for dsolve/numeric to handle to it all more nicely so that setup (whether his way, or mine, or other) can be sidestepped.

There's not much different between our approaches. I used `subs` on a proc because I like to avoid global variables, and I didn't bother to have it return unevaluated on nonnumeric input (hence my uneval quotes when calling it below).

Eg,

restart:
with( plots ):

deqs := { diff(x(t),t) = h*sin(t) + y(t),
          diff(y(t),t) = x(t),
          x(0)=0, y(0)=0}:

integ := dsolve(deqs, numeric, parameters=[h]):

master:=subs(dummy = eval(integ),
             proc(h) dummy('parameters' = [h]);
             dummy end proc):

# A bit like Doug described...
animate(odeplot, ['master'(H), [[t, x(t)], [t, y(t)]], t = 0 .. 3], H = -1 .. 5);

# ...the way I happened to use `master`
display([seq(odeplot(master(H), [[t, x(t)], [t, y(t)]], 0..3,
             'numpoints' = 500, 'thickness' = 2),H=-1..5)],
        'axes' = 'box', 'insequence' = true);

The above is about Doug's `animate` call. Something similar to Doug's `plot3d` call above could also be made to work, also with a little (different) setup.

What I wrote about Matrix form in another thread (from which this NLPSolve Question was branched) pertained specifically to LPSolve and the cost of problem setup.

For NLPSolve there may be quite a different set of purposes to obtaining the Matrix form. For NLPSolve the benefits are not always about efficiency of problem setup (nonlinear problems are often smaller, on practical grounds) but are often more about the solving stage (both efficiency and even solvability).

None of this relates to whether NLPSolve is the best thing to use here, of course.

acer

@alex_01 I understand. You yourself would like to see a pattern visually (perhaps sometimes from a smaller example). People who work regularly with linear algebra and math can often see such patterns from the notation itself, regardless of the size or number of individual equations forming the patterns.

I converted your example from the top of this post, to the Matrix form that I showed, in my head and directly typed into the worksheet. I did not use that internal routine to do it, or to detect the pattern. It was very easy, and took only about ten minutes.

The Matrix form target ([A,b,Aeq,beq],[bl,bu]) is extremely straightforward. It was quite clear how to obtain it from the Matrix-Vector multiplications which you did using W.

The proof is in the pudding. The fact that your Matrix-Vector multiplications were near trivial to convert is reflected in the simple fact that the data Matrices R and ER appear directly as the major portion of the Matrix form Matrix A.

Note that I had also shifted the ur and dr to the left-hand-side since that makes it clea that they are variables just like the w[i]. The rhs's are just numbers (mostly zero). The bottom left parts of A, Aeq are for the W, and handling only two more variable ur and dr required slightly larger matrices and a only few more lines of code.

If you don't think that the two-fold speedup at size n=500 is good enough then you have a simple choice, moving forward: if it's too hard, and the benefit too slight, then simply don't bother with it.

Careful optimization of code often requires several kinds of modifications, picking up a factor of 2 here, and some more there, sometimes leading to more impressive overall improvement. You are the one who originally wrote that you want the code here "as fast and efficient as possible", and that usually means that any modest improvements (taken piecemeal or not) are wanted. If your examples already run fast enough and you don't need them to scale up in size then you aren't forced down that road. Nobody else will care.

[I recall now that we have discussed this kind of thing before...]

@alex_01 I understand. You yourself would like to see a pattern visually (perhaps sometimes from a smaller example). People who work regularly with linear algebra and math can often see such patterns from the notation itself, regardless of the size or number of individual equations forming the patterns.

I converted your example from the top of this post, to the Matrix form that I showed, in my head and directly typed into the worksheet. I did not use that internal routine to do it, or to detect the pattern. It was very easy, and took only about ten minutes.

The Matrix form target ([A,b,Aeq,beq],[bl,bu]) is extremely straightforward. It was quite clear how to obtain it from the Matrix-Vector multiplications which you did using W.

The proof is in the pudding. The fact that your Matrix-Vector multiplications were near trivial to convert is reflected in the simple fact that the data Matrices R and ER appear directly as the major portion of the Matrix form Matrix A.

Note that I had also shifted the ur and dr to the left-hand-side since that makes it clea that they are variables just like the w[i]. The rhs's are just numbers (mostly zero). The bottom left parts of A, Aeq are for the W, and handling only two more variable ur and dr required slightly larger matrices and a only few more lines of code.

If you don't think that the two-fold speedup at size n=500 is good enough then you have a simple choice, moving forward: if it's too hard, and the benefit too slight, then simply don't bother with it.

Careful optimization of code often requires several kinds of modifications, picking up a factor of 2 here, and some more there, sometimes leading to more impressive overall improvement. You are the one who originally wrote that you want the code here "as fast and efficient as possible", and that usually means that any modest improvements (taken piecemeal or not) are wanted. If your examples already run fast enough and you don't need them to scale up in size then you aren't forced down that road. Nobody else will care.

[I recall now that we have discussed this kind of thing before...]

gem 18 is a nice keystroke saver.

acer

@alex_01 You appear to have a fundamental misunderstanding about the purpose of passing an LP problem to LPSolve in Matrix form.

The evidence for this is in your Comment above, where you first form the algebraic form, then convert it, and then pass it in. You even write that you are amazed that anyone would ever try and form the Matrix-form directly (without just using a utility to convert it).

That works. But there is no benefit to doing that conversion. The whole point of using Matrix form is to never form the algebraic equalities/inequalities/bounds-as-ranges.

The way you have done it (algebraic form, then conversion by command, then pass in the Matrix form) is exactly what LPSolve would do for you if you just passed it the algebraic form.

If you have formed the algebraic form, ever, at all, with symbolic variable names like w[i] etc, then you have already wrecked all the advantaages of using the Matrix form for calling LPSolve.

If you cannot understand this then I apologize but I'm afraid that I cannot help you further.

@alex_01 You appear to have a fundamental misunderstanding about the purpose of passing an LP problem to LPSolve in Matrix form.

The evidence for this is in your Comment above, where you first form the algebraic form, then convert it, and then pass it in. You even write that you are amazed that anyone would ever try and form the Matrix-form directly (without just using a utility to convert it).

That works. But there is no benefit to doing that conversion. The whole point of using Matrix form is to never form the algebraic equalities/inequalities/bounds-as-ranges.

The way you have done it (algebraic form, then conversion by command, then pass in the Matrix form) is exactly what LPSolve would do for you if you just passed it the algebraic form.

If you have formed the algebraic form, ever, at all, with symbolic variable names like w[i] etc, then you have already wrecked all the advantaages of using the Matrix form for calling LPSolve.

If you cannot understand this then I apologize but I'm afraid that I cannot help you further.

Are you sure that the current (serial, non-parallel) version of your code is well optimized?

Given (only) that overview of the problem I might guess that Grid would be easier, provided that you have enough memory to hold the multiple full instances of the whole problem instantiated at various parameter value sets. That would also alleviate your having to ensure that your code is Thread-safe.

But you might consider uploading your code here (green arrow in the editor). We have a little experience with optimizing code, though of course there are no promises. Feel free to contact me, if you'd prefer to not post your code publicly.

acer

@alex_01 I think that if you are having trouble understanding how to formulate an NLPSolve problem in Matrix-form then it would be better if you asked in a separate question (giving a full example that explains exactly what you are trying to do).

That way, the two threads remain clear to all late-comers.

@alex_01 I think that if you are having trouble understanding how to formulate an NLPSolve problem in Matrix-form then it would be better if you asked in a separate question (giving a full example that explains exactly what you are trying to do).

That way, the two threads remain clear to all late-comers.

@Joe Riel For this particular problematic file, the two characters are in stranges places. Doesn't one appear between otherwise conjoined numerals 2 and a 4, in a text section?

I think that the GUI should be able 1) to remove such when pasting, 2) not save with such, and 3) remove such when re-opening. Right now it seems to just halt all further reading/parsing and give up, when re-opening, which is weak.

@Joe Riel For this particular problematic file, the two characters are in stranges places. Doesn't one appear between otherwise conjoined numerals 2 and a 4, in a text section?

I think that the GUI should be able 1) to remove such when pasting, 2) not save with such, and 3) remove such when re-opening. Right now it seems to just halt all further reading/parsing and give up, when re-opening, which is weak.