@herclau You don't have to explain to me the use of SVD for Least Squares solving.

Look:

> y.U.Ss.V;

         [0.00520590831446885, 1.00159598094724, 0.999421478032667]

> V^%T.Ss.U^%T.y;

                            [0.00520590831446883]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

Those two above are just transpositions of each other. Earlier, I suggested you use y.U.Ss.V. But of course you can use V^%T.Ss.U^%T.y as well, because they are just transpositions of each other! I simply chose the version which required no additional cost of actually doing transpositon.

You constructed Ss as the inverse of DiagonalMatrix(S[1 .. 3]) so of course Ss is symmetric. And so when multiplying y by V^%T.Ss.U^%T you are (least squares) solving M.x=y.

In your P2 you actually used V.U^%T.y.Ss which is no good. I don't know where you get that formula from.

Remember, SingularValues returns Vt which is the transpose of V. Forgetting that is not enough to generate your strange formula used in your P2, I think.

For fun (but perhaps not optimally using the 'thin' option, I didn't check), since Maple uses SVD to compute the non-square float pseudo-inverse,

> MatrixInverse(M).y;

                            [0.00520590831446895]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

@herclau You don't have to explain to me the use of SVD for Least Squares solving.

Look:

> y.U.Ss.V;

         [0.00520590831446885, 1.00159598094724, 0.999421478032667]

> V^%T.Ss.U^%T.y;

                            [0.00520590831446883]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

Those two above are just transpositions of each other. Earlier, I suggested you use y.U.Ss.V. But of course you can use V^%T.Ss.U^%T.y as well, because they are just transpositions of each other! I simply chose the version which required no additional cost of actually doing transpositon.

You constructed Ss as the inverse of DiagonalMatrix(S[1 .. 3]) so of course Ss is symmetric. And so when multiplying y by V^%T.Ss.U^%T you are (least squares) solving M.x=y.

In your P2 you actually used V.U^%T.y.Ss which is no good. I don't know where you get that formula from.

Remember, SingularValues returns Vt which is the transpose of V. Forgetting that is not enough to generate your strange formula used in your P2, I think.

For fun (but perhaps not optimally using the 'thin' option, I didn't check), since Maple uses SVD to compute the non-square float pseudo-inverse,

> MatrixInverse(M).y;

                            [0.00520590831446895]
                            [                   ]
                            [   1.00159598094724]
                            [                   ]
                            [  0.999421478032667]

I wonder whether this is due to some inefficiency related to (lack of decent) buffering. I guess that I mean: it doesn't sound so good to read from disk a byte at a time (even if not importing integer[1] and even if endianness is a concern). I wonder what the performance would be from just using something like Linux fread (buffering according to cache size, say)?

acer

I believe that this has been mentioned before. But a reminder serves its purpose.

p.s.! Welcome back, Alec, whether here for long or not.

acer

@herclau So first one must correct the cos term your P2 and P3, by multiplying the call to cos() by a factor of 0.5 as mentioned. After doing so then one can get a plot that matches your followup plot on the restricted range 4.0 to 5.7, as above. OK. And then you are unsatisfied with these results of your by-hand use of SVD to do the least-squares fitting, demonstrated by how the red curve is not matching well in that restricted plot.

Did you want,

> y.U.Ss.V;

         [0.00520590828213008, 1.00159598094467, 0.999421478049389]

instead of your,

> V.U^%T.y.Ss;

          [0.0103501902515756, 0.967599719533608, 1.05732003209726]

?

Using y.U.Ss.V gets the overlaid plot below, where the red and green curves coincide so much that only the green curve (on top) is easily visible,

acer

@herclau So first one must correct the cos term your P2 and P3, by multiplying the call to cos() by a factor of 0.5 as mentioned. After doing so then one can get a plot that matches your followup plot on the restricted range 4.0 to 5.7, as above. OK. And then you are unsatisfied with these results of your by-hand use of SVD to do the least-squares fitting, demonstrated by how the red curve is not matching well in that restricted plot.

Did you want,

> y.U.Ss.V;

         [0.00520590828213008, 1.00159598094467, 0.999421478049389]

instead of your,

> V.U^%T.y.Ss;

          [0.0103501902515756, 0.967599719533608, 1.05732003209726]

?

Using y.U.Ss.V gets the overlaid plot below, where the red and green curves coincide so much that only the green curve (on top) is easily visible,

acer

@hirnyk ...that doesn't really explain how the smooth overlaid curve in the parent post was obtained in Mathematica. Which is presumably why the Answer above already linked to that help page you snapped. What value of Mma's ListPlot's InterpolationOrder would fill in such a good a match for sin+cos/2 by mere splines, in the tight corners where there are few data points?

While on topic, the original post had another funny thing, at least to my eye. The 'errors' introduced by the addition of that final random array, with points between 0.0 and 0.1, look a little small on the plot. it doesn't look like the points vary from the sin-cos/2 function by that much (ie, not by 0.05 in y direction on average). Maybe my eyes are going. (Note, the Answer above uses 0.0..0.01 for those 'errors'. But it is the original Mma's 0.0..0.1 range for them that seems under-represented in the plotted data.)

Or maybe the original plot was made with some (other) command variant(s).

@hirnyk ...that doesn't really explain how the smooth overlaid curve in the parent post was obtained in Mathematica. Which is presumably why the Answer above already linked to that help page you snapped. What value of Mma's ListPlot's InterpolationOrder would fill in such a good a match for sin+cos/2 by mere splines, in the tight corners where there are few data points?

While on topic, the original post had another funny thing, at least to my eye. The 'errors' introduced by the addition of that final random array, with points between 0.0 and 0.1, look a little small on the plot. it doesn't look like the points vary from the sin-cos/2 function by that much (ie, not by 0.05 in y direction on average). Maybe my eyes are going. (Note, the Answer above uses 0.0..0.01 for those 'errors'. But it is the original Mma's 0.0..0.1 range for them that seems under-represented in the plotted data.)

Or maybe the original plot was made with some (other) command variant(s).

@Maple_Maple I am not sure what you mean by "related at T". The solutions I showed each specified what T (along with all the other variables) had to be.

Do you mean that you think that the only valid answers are ones in which the other variables are defined in terms of T? If that's so, then run,

solve({fff},AllSolutions);

yourself and see what from the whole solution set might satisfy your need.

You could try your luck with something like,

eliminate({fff}, indets({fff},name) minus {constants} minus {T});

but I recall having a problem with that -- some eventual murky error message issued from `eliminate/recursive`.

The command solve({fff},AllSolutions) eventually returns, after about an hour for me in Maple 14, giving even more solutions (involving RootOf, some of them long) but also including those I showed above.

note to self: there are some duplicate entries in the solutions I showed above (differing only in the particular name of the assumed _Z parameters).

@Maple_Maple I am not sure what you mean by "related at T". The solutions I showed each specified what T (along with all the other variables) had to be.

Do you mean that you think that the only valid answers are ones in which the other variables are defined in terms of T? If that's so, then run,

solve({fff},AllSolutions);

yourself and see what from the whole solution set might satisfy your need.

You could try your luck with something like,

eliminate({fff}, indets({fff},name) minus {constants} minus {T});

but I recall having a problem with that -- some eventual murky error message issued from `eliminate/recursive`.

The command solve({fff},AllSolutions) eventually returns, after about an hour for me in Maple 14, giving even more solutions (involving RootOf, some of them long) but also including those I showed above.

note to self: there are some duplicate entries in the solutions I showed above (differing only in the particular name of the assumed _Z parameters).

@epostma Erik, you are a gentleman and a scholar.

An entertaining post!

That financial summary at the end is the best part.

acer

@Alejandro Jakubi Yes, all that information is provided by CodeTools:-Usage if I invoke it like,

CodeTools:-Usage(...something...,quiet,output=[output,realtime,cputime,bytesused,bytesalloc]);

But that's more effort and typing than I care for. That's why I like the idea of `tools/bench`, or something shorter still and universally avalable which has the default bahaviour that I want. (You write yourself that you like `tools/bench`, but of course you too could just use CodeTools:-Usage with the appropriate options.) It's just nicer to have a short alias. Of course I could stick a short alias for it in all my personal initialization files, and merely suffer a little when using other nonnetworked machines. I guess that it's human nature to imagine that everyone else wants the same defaults as oneself -- where I, you, John, etc might be naturally prone to that just like most others are.

Which reminds me of something I've been wondering: if Maple can share worksheets on the Cloud then why not also allow "login" via Maple's GUI to load GUI preferences (there are so many!) and even .mapleinit material?

Since `tools/bench` has been mentioned here, I'd like to add that I usually want the more pieces of information, including change in time, change in memory "used", and change in memory allocated. Eg,

(st,kbu,kba):=time(),kernelopts(:-bytesused),kernelopts(:-bytesalloc):
 ...some computation...
time()-st,kernelopts(:-bytesused)-kbu,kernelopts(:-bytesalloc)-kba:

My reason is that the change in bytesused shows more how much memory is being "managed", while I am often more interested in how much extra memory allocation on Maple's behalf is brought about by the computation. (I do realize that there is interplay between the concepts, and that once the gcfreq limit is hit Maple might collect memory before it allocates more. So care is needed in interpreting bytesused and bytesalloc if done in mid-session.)

I also find that, w.r.t. measuring Thread'd applications, I'm also interested in the change in wall-clock time. Now, the time() command measures total cycles used by all Threads, and adds them up. So in order to effectively estimate true elapsed computation time I have to resort to time[real] (and hope that other processes on my machine aren't slowing it down and affecting the wall-clock difference). So these days I find I do,

(swt,st,kbu,kba):=time[:-real](),time(),kernelopts(:-bytesused),kernelopts(:-bytesalloc):
 ...some computation...
time[:-real]()-swt,time()-st,kernelopts(:-bytesused)-kbu,kernelopts(:-bytesalloc)-kba:

I have another quibble. I often want to benchmark a computation and see both the performance stats as well as the actual computational result. I find the invocation time(...some command...) disappointing in that it doesn't provide me with the actual result. And I find that `tools/bench` suffers in the same way. I'd prefer it if `tools/bench` also included its actual result eval(e) as part of the returned expression sequence. I don't want to have to make my computational argument to `tools/bench` be wrapped in `assign`.

As for the ::uneval topic, I'll mention this use (which I'd only mentioned to Jacques privately before I saw that this is a popular Post). It's not something most people will ever want to do I think, but basically one can uneval top-level commands and then manipulate them as one wishes. Of course this too does not prevent automatic simplification, as Claire has mentioned.

I enjoyed this series of related posts very much.

They do remind me of a related topic which was discussed (a few different times) on this site some years ago. it was the issue of how to efficiently generate many different medium-sized samples of a given distribution. At that time, the focus was mostly on reducing overhead within the Statistics package, or reducing the overhead of garbage collecting (a.k.a. memory management) of Array/Vector containers for such samples. The judgement then (IIRC) was in preference of generating a single sample up front -- as large as one can bear in terms of memory allocation.

Now it seems to me that you have mostly been discussing user-defined distributions, and that presumably there are some better/dedicated methods for sample well-known stock distributions. Please correct me if that is wrong. So my related question is: for these relevant distributions how can the overhead of duplicated recomputation of all these subdivisions (that you have shown above) be avoided in the case that one wishes to resample a distribution many different times? Eg, if I wish to generate finite-sized samples of 10,000 entries, repeated say 100,000 times.

acer