Hi Robert,

That Warning about first-order conditions is (somewhat not well-known) bug that crops up when using the operator form calling sequence of NLPSolve. The (Trackered) problem seems to relate to numeric computation of derivatives, ie. gradients of objectives and jacobians of constraints.

I have had some prior success with two types of workaround for this bug. One is to use expression form, rather than operator form. In this case, passing unquoted E(c0, c1, c2) to NLPSolve results in an unevaluated integral, which seems "good enough" to serve as an expression in this particular need.

The other type of workaround I've used before is to manually and explicitly supply procedures for the gradient (and jacobian if necessary for non-simple constraints). For that, I have previously used fdiff, because in such cases I expect that NLPSolve has already tried codegen:-Gradient (A.D.) and failed. I have previously mailed examples of such to Joe R., but I'd like to find time to post them here too.

Of course, I too replaced `int` with `Int` inside the `F` procedure body.

I also specified the variables, in my NLPSolve call. Ie,

> Optimization:-NLPSolve(E(c0, c1, c2),
>        variables=[c0,c1,c2],
>        initialpoint=[c0=.1035805,c1=1.001279,c2=-.561381]);
Error, (in Optimization:-NLPSolve) complex value encountered

I didn't yet investigate whether this last error is truthful, or whether or not some tiny imaginary component from the numeric integral result might not be just a negligible artefact of floating-point computation, or how to work around such.

acer

Shouldn't the call to Probability in proc f instead be to ST:-Probability ?

acer

Shouldn't the call to Probability in proc f instead be to ST:-Probability ?

acer

It seems plausible that there is some (other) example for which the "I" terms are not yet collected in the simplify@evalc result, while the direct solve result is not in the desired form.

Anyway, this is all verysimple tweaking, so arguing about it is mere quibbling. And beauty is in the eye of the beholder. So, yet another quibble: someone else might prefer the result of,

map((normal@Re)+(normal@Im)*I,simplify(evalc([sol])));

acer

It seems plausible that there is some (other) example for which the "I" terms are not yet collected in the simplify@evalc result, while the direct solve result is not in the desired form.

Anyway, this is all verysimple tweaking, so arguing about it is mere quibbling. And beauty is in the eye of the beholder. So, yet another quibble: someone else might prefer the result of,

map((normal@Re)+(normal@Im)*I,simplify(evalc([sol])));

acer

While the random 10x10 Matrix might likely have a full set of linearly independent eigenvectors, it could of course just be that the OP used RandomMatrix for illustration. It may be that for the OP's typical problems the chance of the input being a defective Matrix is more likely. (We don't yet know.)

Such decomposition approaches to computing the floating-point exponential generally fail for defective Matrices. See Method 14 of section 6 of the well-known paper by Moler and Van Loan. I doubt that this method gets "better" by applying it to the mixed float-symbolic problem, which is what we have here for an unassigned symbol `t`. (For a float Matrix, one might be able to compute the condition number quickly, as an initial test.)

Naturally, I fully expect that Robert is aware of all this. I only mention it for the possible benefit of others.

acer

While the random 10x10 Matrix might likely have a full set of linearly independent eigenvectors, it could of course just be that the OP used RandomMatrix for illustration. It may be that for the OP's typical problems the chance of the input being a defective Matrix is more likely. (We don't yet know.)

Such decomposition approaches to computing the floating-point exponential generally fail for defective Matrices. See Method 14 of section 6 of the well-known paper by Moler and Van Loan. I doubt that this method gets "better" by applying it to the mixed float-symbolic problem, which is what we have here for an unassigned symbol `t`. (For a float Matrix, one might be able to compute the condition number quickly, as an initial test.)

Naturally, I fully expect that Robert is aware of all this. I only mention it for the possible benefit of others.

acer

"MMA", sometimes given as "Mma", is a popular abbreviation for Mathematica,

acer

"MMA", sometimes given as "Mma", is a popular abbreviation for Mathematica,

acer

Are N.T.'s papers on these topics purely for numeric (and not symbolic) computation? Maple doesn't seem too slow for that. It seems to be the symbolic computation and interpolation (even though here using very quickly computed float eigenvalues) that gets bogged down.

acer

Are N.T.'s papers on these topics purely for numeric (and not symbolic) computation? Maple doesn't seem too slow for that. It seems to be the symbolic computation and interpolation (even though here using very quickly computed float eigenvalues) that gets bogged down.

acer

Note that you are computing a symbolic result, not a pure floating-point result. Actually, it's a mixed float-symbolic problem, which often is a special kind of animal. And the result will be very lengthy. Note, too, that pure floating-point Matrix exponentiation should be pretty fast.

So perhaps you should ask yourself: what do you expect to do with such a result? Often, such lengthy symbolic results don't give much insight on their own. One often ends up numerically instantiating (possibly many) values for the symbol (here, the `t`) in some other calculation or plot.

By my rough estimation, one can do about 1500 such 10x10 complex[8] pure floating-point Matrix exponentials in the same time that one can do the given symbolic example.

So rather than try for the general univariate, symbolic result (which might just end up being evaluated at many individual t values) perhaps it would be more efficient to alter one's methodology and compute the many individual purely floating-point Matrix exponentials directly.

acer

Try passing `plot` the operator F1, instead of the expression that is returned by F1(t). Ie,

plot(F1, -2..2);

I also saw quick success (even using the expression form) after simplifying the integrand (with `simplify`) following a conversion of the .2 to exact rational 2/10, or injecting a lower tolerance option into the evalf@Int call. The times when it was faster than your original call (by about a factor of 5) were pretty much each the same.

So, that can make plotting F1 faster, and fill those gaps. But it's still not very fast to plot F1 over a range of x. That's because computing F1 at each value of x value involves a separate infinite oscillatory numerical integral. There isn't a good general scheme for such in Maple yet.

acer

So, first you are doing this?

a:=readdata(`testdata.txt`,string,16):

If you've done that, then you could extract the contents of each list in `a` by replacing your lengthy command,

seq(op(i,op(1,a)),i=1..nops(op(1,a)));

with simply this,

op(a[1]);

As an alternative to using StringTools to remove the W characters, perhaps try something like this,

a:=seq(map(t->`if`(t="W",NULL,parse(t)),T),
       T in readdata(`testdata.txt`,string,16)):

That gives you a sequence of lists. Remember: you can produce a sequence from a list by simply applying op to it. Ie, at that point,

op(a[5]); # sequence, from 5th line/list in `a`

But I can't see why you need to convert the lists to sequences right away. A list is often just as easy to manipulate. I understand, that you intend to somehow plot the data.

acer

So, first you are doing this?

a:=readdata(`testdata.txt`,string,16):

If you've done that, then you could extract the contents of each list in `a` by replacing your lengthy command,

seq(op(i,op(1,a)),i=1..nops(op(1,a)));

with simply this,

op(a[1]);

As an alternative to using StringTools to remove the W characters, perhaps try something like this,

a:=seq(map(t->`if`(t="W",NULL,parse(t)),T),
       T in readdata(`testdata.txt`,string,16)):

That gives you a sequence of lists. Remember: you can produce a sequence from a list by simply applying op to it. Ie, at that point,

op(a[5]); # sequence, from 5th line/list in `a`

But I can't see why you need to convert the lists to sequences right away. A list is often just as easy to manipulate. I understand, that you intend to somehow plot the data.

acer