@alex_01 

You write that you worry about the performance of your LS approach, by which you likely mean its robustness and correctness since your concern covers the two cases where its error norm is much higher. But you have not acknowledged that using the explictly formed Matrix inverse of (A^%T.A) is not the best numerical way to compute (A^%T.A)^(-1).A^%T which was illustrated by numeric example above. That exponent-(-1) is just a convenient way to denote the Matrix inverse for publication. That does not mean that it explicit inverse calculation is always numerically appropriate. But you can use LinearSolve and still validly claim that what you're doing is computing that same formula. And so you ought to, as it might clear up one of the two cases of your concern over your LS approach.

You write that you are concerned about the performance of the SVD method as used by the LeastSquares command for the over-determined case (m>n), and you give an example where it (and the other) approaches all produce a noticable error norm. But what value did you expect? It won't typically be anywhere close to zero, for the kinds of random entries you're using in the Matrix, because the m>n (df>0) system will very likely be overdetermined. That term means that there cannot be a solution which solves the system exactly -- an overdetermined system is inconsistent, and its least squares solution will have a non-zero error norm.

You may have to show the existence of significantly a better solution, before the method=SVD solution is cast in serious doubt.

You asked about eigenvalues. If you mean you want the eigenvalues of (A^%T.A) then those are the singular values of Matrix A. Note that one does not necessarily use the time and memory to explicit form (A^%T.A) as a Maple Matrix in order to compute the singular values of floating-point Matrix A. The computational work is done in compiled C functions, and clever workspace use/re-use may also mean that some forms are not produced in any stand-alone/returnable fashion.

In your first comment you claim that your Maple 13 does not have an SVD method implemented within LeastSquare. But that is untrue.

> infolevel[LinearAlgebra]:=1:

> LinearAlgebra[LeastSquares](X, Y, method=SVD):

SingularValues: calling external function
SingularValues: NAG hw_f08kef
SingularValues: NAG hw_f08kff
SingularValues: NAG hw_f08kff
SingularValues: NAG hw_f08mef
LeastSquares: calling external function
LeastSquares: NAG hw_f06paf
LeastSquares: NAG hw_f06pbf
LeastSquares: NAG hw_f06paf

> LinearAlgebra[LeastSquares](X, Y):

LeastSquares: calling external function
LeastSquares: NAG hw_f08aef
LeastSquares: NAG hw_f08agf
LeastSquares: NAG hw_f07tef

> kernelopts(version);
          Maple 13.02, X86 64 WINDOWS, Oct 12 2009, Build ID 436951

What was more recently added was the ability of LeastSquares to make use of the so-called "thin SVD" which makes it much more efficient in the greatly overdetermined case m>>n. (See here for a stand-alone implementation of that, outside of LinearAlgebra.)

If you want to see a more user-level implementation of the SVD approach then you could do worse than look over several questions posed by member herclau in the second half of 2010. eg. 1, 2, 3.

Note that in some cases the QR method may be faster, while the SVD method might be more robust. This is a commom dichotomy: better speed vs. better assurance of correctness.

You might think that everyone else should be able to figure out, without explanation beforehand, your intended meanings that df = num.rows - num.cols, or that "ey" means "Estimated Y" and is computed via A.x, and that data1 is the solution in the code while `x` is the solution elsewhere, and X is the Matrix in the code while A is the Matrix elsewhere, and that sometimes Y means lowercase y while X does not mean lowercase x. But it does make it harder to read. I think that it isn't as polite as can be, for your audience.

acer

@Markiyan Hirnyk You misunderstand my points, I believe. Maple has not had session dependent set ordering based on memory addresses in several major releases, so it is not relevent that you get the same order upon multiple reexecution here.

What I was trying to say is that the order will depend upon the lexicographic order (although that too can vary with option `setsort` as Maple is launched). So, choose some other name for the parameter, ie. Z instead of b, and its position in the set or list will be different. So choosing B[2] is ad hoc and will not reliably produce the same kind of result for other examples.

And also, one can choose a formal parameter x and along with additional (scoped) b, or a formal parameter b along with additional (scoped) x. And both of this two could produce the same list B, in which case the code that you showed, alone, will be insufficient to determine which is the additional parameter. That is why, when manipulating the inert form, I deliberately excluded the formal parameters (of an operator present) before calling FromInert.

This is just another instance of the fact that lexicographic set ordering is, essentially, about as indeterminate as was the old address-based set ordering. The latter is indeterminate due to what may or may not have occurred previosuly in the session, while the former is indeterminate because the user has free reign to choose his or her variable names at will.

Post-conversion to a list does not pertain to or alter the above.

@Markiyan Hirnyk You misunderstand my points, I believe. Maple has not had session dependent set ordering based on memory addresses in several major releases, so it is not relevent that you get the same order upon multiple reexecution here.

What I was trying to say is that the order will depend upon the lexicographic order (although that too can vary with option `setsort` as Maple is launched). So, choose some other name for the parameter, ie. Z instead of b, and its position in the set or list will be different. So choosing B[2] is ad hoc and will not reliably produce the same kind of result for other examples.

And also, one can choose a formal parameter x and along with additional (scoped) b, or a formal parameter b along with additional (scoped) x. And both of this two could produce the same list B, in which case the code that you showed, alone, will be insufficient to determine which is the additional parameter. That is why, when manipulating the inert form, I deliberately excluded the formal parameters (of an operator present) before calling FromInert.

This is just another instance of the fact that lexicographic set ordering is, essentially, about as indeterminate as was the old address-based set ordering. The latter is indeterminate due to what may or may not have occurred previosuly in the session, while the former is indeterminate because the user has free reign to choose his or her variable names at will.

Post-conversion to a list does not pertain to or alter the above.

@Markiyan Hirnyk Let's also recall that the positioning in the indets result may depend upon set ordering.

So determining the placement is problematic.

> restart:
> h := Statistics:-Distribution(PDF = (x-> Dirac(x-b))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "arrow", "b", "x", "operator"}

> h := Statistics:-Distribution(PDF = (x-> Dirac(x-A))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "A", "Dirac", "arrow", "x", "operator"}

> h := Statistics:-Distribution(PDF = (x-> Dirac(x-Z))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "Z", "arrow", "x", "operator"}

Since this may have been intended for user-defined distribution, the formal parameter name(s) or the PDF operator may also be ad hoc,

> h := Statistics:-Distribution(PDF = (b-> Dirac(b-x))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "arrow", "b", "x", "operator"}

In this last example, we would need to determine the parameter name as well, so as to not misidentify.

This all relates to why I originally sieved out the formal parameter names and some other atomics, when processing the inert form.

But I suppose that all this is moot now, as member icegood may choose another interface for his code's users.

@Markiyan Hirnyk Let's also recall that the positioning in the indets result may depend upon set ordering.

So determining the placement is problematic.

> restart:
> h := Statistics:-Distribution(PDF = (x-> Dirac(x-b))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "arrow", "b", "x", "operator"}

> h := Statistics:-Distribution(PDF = (x-> Dirac(x-A))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "A", "Dirac", "arrow", "x", "operator"}

> h := Statistics:-Distribution(PDF = (x-> Dirac(x-Z))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "Z", "arrow", "x", "operator"}

Since this may have been intended for user-defined distribution, the formal parameter name(s) or the PDF operator may also be ad hoc,

> h := Statistics:-Distribution(PDF = (b-> Dirac(b-x))):
> B := indets(ToInert(eval(h:-PDF)), atomic);           
               B := {1, "Dirac", "arrow", "b", "x", "operator"}

In this last example, we would need to determine the parameter name as well, so as to not misidentify.

This all relates to why I originally sieved out the formal parameter names and some other atomics, when processing the inert form.

But I suppose that all this is moot now, as member icegood may choose another interface for his code's users.

Thanks for uploadind your worksheet. Please note that we cannot yet run it, as it relies on some procedures like procBombcurve and (presumably assigned names of expressions like) PreBombcurve and PostBombcurve. Is it possible for you to upload everything required to run the worksheet? (Uploads and links to multiple files could be included, if you at liberty to do so.)

Also, I notice that the lines of the worksheet must be executed out of order, even if the missing bits were available. For example it executes a test of procedures `p` and `procRsum` before it defines/assigns procedures `procPatrientres` and `procC14incell`. That is not critical, but may confuse someone trying to run it.

It's slightly unweildy for procedures such as procPatientres to assign to globals like `residual`, when it might be more straightforward to just have it return the computed value for that. (And, in turn, any other procedure which called that could assign the result to its own local `residual`. Doing it in that way can reduce mix-ups, and keep the global name space clean and avoid need for unassignments at the top-level to avoid inadvertant mistakes.

acer

I really like purely exact approaches to such problems: no hidden float log, no `floor` or `ceil`, no outright or hidden `evalf` calls at all. The above method seems to be like that -- just exact multiplication and `op`. A reduced risk of the task's being liable to mistaken exclusion/inclusion due to roundoff or other numerical difficulty.

And it's probably quite fast.

And trailing-zero removal may be applied afterwards.

I'm supposing that `union` itself is not invoking evalf, and is relying only on evalb here.

Of course, for such tasks the speed can depend heavily on details. Is it for cases where a) there are many duplicates of just a few "unique" values, b) few duplicates of many "unique" values, or c) some mixed scenario.

acer

I really like purely exact approaches to such problems: no hidden float log, no `floor` or `ceil`, no outright or hidden `evalf` calls at all. The above method seems to be like that -- just exact multiplication and `op`. A reduced risk of the task's being liable to mistaken exclusion/inclusion due to roundoff or other numerical difficulty.

And it's probably quite fast.

And trailing-zero removal may be applied afterwards.

I'm supposing that `union` itself is not invoking evalf, and is relying only on evalb here.

Of course, for such tasks the speed can depend heavily on details. Is it for cases where a) there are many duplicates of just a few "unique" values, b) few duplicates of many "unique" values, or c) some mixed scenario.

acer

@alex_01 I could be mistaken in thinking that your problem is necessarily nonlinear. And you are right: solving it could likely be faster if it could be reformulated as LP. (There are constraints to accomodate, albeit simple bounds.) I'll try and give it some proper thought.

Ideally, the help-pages and Example worksheets are supposed to explain how to get Matrix form, with examples that explain it all. I agree, it can be messy. But as you say, seeing patterns is still not enough: one has to fully understand the structures in order to be able to get it for any new problem.  If I can find time, I'll convert this example, forcing use of NLPSolve.

Cheers. Sorry if I've been overbearing in this thread.

@alex_01 I could be mistaken in thinking that your problem is necessarily nonlinear. And you are right: solving it could likely be faster if it could be reformulated as LP. (There are constraints to accomodate, albeit simple bounds.) I'll try and give it some proper thought.

Ideally, the help-pages and Example worksheets are supposed to explain how to get Matrix form, with examples that explain it all. I agree, it can be messy. But as you say, seeing patterns is still not enough: one has to fully understand the structures in order to be able to get it for any new problem.  If I can find time, I'll convert this example, forcing use of NLPSolve.

Cheers. Sorry if I've been overbearing in this thread.

@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