## 1394 Reputation

13 years, 134 days
Maplesoft

## Social Networks and Content at Maplesoft.com

I am the manager of the Mathematical Software Group, working mostly on the Maple library. I have been working at Maplesoft since 2007, mostly on the Statistics and Units packages and on contract work with industry users. My background is in abstract algebra, in which I completed a PhD at Eindhoven University of Technology. During my studies I was always searching for interesting questions at the crossroads of math and computer science. When I came to Canada in 2007, I had nothing but a work permit, some contacts at Maplesoft (whom I had met at a computer algebra conference a year earlier), and a plan to travel around beautiful Canada for a few months. Since Maplesoft is a company that solves more interesting math and computer science related questions than just about anywhere else, I was quite eager to join them, and after about three months, I could start.

## Yes...

This is intentional behaviour; it is sometimes used inside the library to select between different versions of the same command, or to select certain options. For example, see ?evalf or ?map or ?Statistics/StandardError. map can even take multiple indices.

The way to access this information is by examining ?procname. That is the name by which the current procedure was invoked. In the example above, procname will be equal to a[b]. To check whether a specific index was given, you could do something like this:

`a := proc(\$)  if type(procname, 'indexed') and [op(procname)] = ['b'] then    printf("invoked with index b\n");  elif not type(procname, 'indexed') then    printf("invoked without indices\n");  else    error "invalid input: expected %1's index to be 'b' or absent, but received %2", op(0, procname), [op(procname)];  end if;end proc;`

Note that op(procname) is typically enclosed in brackets to make it a list, because you don't know a priori whether someone might try to invoke your procedure with multiple indices.

Hope this helps,

Erik Postma
Maplesoft.

## See this help page for investigation hin...

Hi Manohar,

The default method of integration for the stiff case is the Rosenbrock algorithm. The help page for  ?dsolve/rosenbrock lists some details that might be helpful. As a first step, I would set infolevel[dsolve] pretty high, to 5 or even 10, in order to get an idea of where it is that dsolve is spending all its time. The very first clue is, when does it take so much time? Is it during the initial call to dsolve, or is it later when you ask it to evaluate the solution at some point? (I'm assuming you're not using the calling sequence that includes the 'range' option.)

If the time is taken during the initial call, the issue might be computing the Jacobian; with a bit of work and use of the 'jacobian' option, you might be able to get around that. If the time is taken during actual evaluation of the solution, you could try playing with the tolerance options ('abserr', 'relerr', and 'initstep') if you have some flexibility in that.

Finally, the Rosenbrock method is most suitable for the case where you're not looking for extremely high accuracy - if you are already setting the tolerance options in order to get very high accuracy, then you might be better off using method=lsode[x] for some value of x. Or is that what you meant where you were talking about the "islode solver"?

Hope this helps,

Erik Postma
Maplesoft.

## Are you sure this is what you want to co...

Hi alex_01,

One thing I found somewhat surprising (from a mathematical / statistical point of view - I can see how it could seem natural from an experimentalist point of view) is that after taking the subsets of data that are sufficiently far away from the mean, you take just the correlation coefficient of these data sets - which means: where X and Y are the subsets of data. I would think you might want to consider using the mean of the whole data set together, instead of Mean(X) and Mean(Y). Anyway - maybe it doesn't make much of a difference.

Either way, the condition on which the correlation is conditional, is that and . (If N goes to infinity, then the sample mean approaches the population mean arbitrarily well (for the normal distribution), so we'll ignore that difference.)

Now we can define the bivariate normal PDF (see e.g. Wikipedia), set mu to 0 and sigma to 1 (any different situation can be obtained from that by translation and scaling), restrict its domain to X and Y both greater than w, and start writing down integrals. I did this for a while and got the following horrendous expression for the covariance:

`Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*(x-(Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*x/(Int(Limit(Pi^(1/2)*exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2))-erf((w-w*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), y = w .. infinity)), x = w .. infinity), y = w .. infinity)))*(y-(Int(Int(exp(-(1/2)*(x^2+y^2-2*rho*x*y)*(1-rho^2))*x/(Int(Limit(Pi^(1/2)*exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2))-erf((w-w*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), y = w .. infinity)), x = w .. infinity), y = w .. infinity)))/(Int(Limit(Pi^(1/2)*exp(-(1/2)*(-1+rho^2)^2*y^2)*(erf((x-x*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2))-erf((w-w*rho^2+rho^3*y-rho*y)/(2-2*rho^2)^(1/2)))/(2-2*rho^2)^(1/2), x = infinity), y = w .. infinity)), x = w .. infinity), y = w .. infinity)`

I'm not very optimistic that we are going to be able to figure out the behaviour of this expression as w grows for various values of rho. It's not unlikely that there are other, better, explanations, though.

Erik.

## Use it as a component...

@vicaut : You can use Kamel's idea as a component to building a solution that includes all the roots. In fact, I would suggest using the ?RootFinding/NextZero command; it has somewhat better numerical behaviour and its behaviour is a bit better specified than fsolve for this type of question (where you want to "avoid" many roots). Define the following (Maple 15) procedure:

`allRoots := proc(f :: procedure, rng :: range, \$)local result, root;  result := Vector(0, 'datatype = float');  root := RootFinding:-NextZero(f, lhs(rng), 'maxdistance' = rhs(rng) - lhs(rng));  while root <> FAIL and root <= rhs(rng) do    result(numelems(result) + 1) := root;    root := RootFinding:-NextZero(f, root, 'maxdistance' = rhs(rng) - root);  end do;  return result;end proc;`

Now allRoots(rhs(ODE_4), 10 .. 20) should give you a Vector with all the (100 or so) roots that you were talking about.

Hope this helps,

Erik Postma
Maplesoft.

## Here is a light weight approach...

Here is an alternative approach that is fully in-memory. It won't give you a very beautiful plot, but the idea is there. My example (in Maple 15) uses the ODE with parameter a.

Actually, in production code, I'd use option cache rather than remember because it is less likely to blow up with memory use. However, here option remember will use less memory than the plot structures involved, so it will do for this example.

`sys := [a*diff(f(x), x, x) + diff(f(x),x) + f(x) = 0, f(0)=0, D(f)(0) = 1];solver := dsolve(sys, numeric, parameters=[a]):solvedf := x -> eval(f(:-x), solver(x));zerolist := proc(aa)global solver, solvedf, parameters, a;option remember;local x0, result;  solver(parameters=[a=aa]);  x0 := RootFinding:-NextZero(solvedf, 0);  result := Vector(0);  while x0 < 20 do    result(numelems(result) + 1) := x0;    x0 := RootFinding:-NextZero(solvedf, x0);  end do;  return result;end proc;dispatcher := i -> proc(aa)local result;  result := zerolist(aa);  if numelems(result) >= i then    return result[i];  else    return NULL;  end if;end proc;plot(dispatcher(1), 1/5 .. 5):# This does adaptive plotting, populating the remember table of zerolist in an appropriate way.rememberedXValues := sort([indices(op(4, eval(zerolist)), 'nolist')]);maxlength := max(map(numelems, [entries(op(4, eval(zerolist)), 'nolist')]));plot([seq(dispatcher(i), i = 1 .. maxlength)], 1/5 .. 5, 'sample' = rememberedXValues,      'adaptive = false', 'view = [DEFAULT, 0..15]');`

This results in the following graph: Hope this helps,

Erik Postma
Maplesoft.

## Use Statistics...

Hi zbay,

The easiest solution to achieve what you ask for is to simply run the Statistics:-Mean and Statistics:-StandardDeviation procedures. So you load the Statistics package (by typing with(Statistics): ), enter your values into a  ?Vector V, then type Mean(V), StandardDeviation(V); to get the two values.

If you want to use the procedure you have written already, then I'd adapt it as follows:

`myStandardDeviation := proc(S)local x, mean;  mean := Statistics:-Mean(S);  # or if you prefer, mean := add(x, x in S) / numelems(S);  return sqrt(add((x - mean)^2, x in S) / (numelems(S) - 1));end proc;`

Note the following things:

1. In cases such as this, where you're adding a concrete set of numbers instead of constructing a symbolically described sum, ?add is more appropriate than  ?sum .
2. You forgot to specify which the running variable is in the sum.
3. I use Maple 15's command  ?numelems , which makes this automatically work for lists, Vectors, sets, and essentially all data structures you can think of.

Hope this helps,

Erik Postma
Maplesoft.

## The number of variables is not going to ...

Hi mutmurat,

It's highly unlikely that the number of variables is going to be the limiting factor in your computation. You can easily, say, differentiate simple expressions in thousands of variables. However, if the computation requires a lot of computations of large, complicated expressions, which are subsequently used in building ever larger expressions, there are circumstances where this will quickly grind to a halt. (A notorious, innocent-seeming example is inverting a symbolic matrix. Some of these matrices look quite innocent but symbolic inverses cannot fit into memory of a current machine.)

I'm aware that this does not really answer your underlying question. Your best bet is probably simply trying your computation. Be sure to ask here again if you get stuck, and provide lots of details about what you're doing!

Hope this helps,

Erik Postma
Maplesoft.

## RealDomain...

Try the RealDomain package.

`RealDomain:-solve(x^4=1); # returns -1, 1`

It's not perfect, but it's pretty good.

Hope this helps,

Erik Postma
Maplesoft.

## Normalizer and Testzero...

One way that you can sometimes attack these problems is by reducing the amount of normalization that Maple does. For matrix inversion, for example, every potential pivot has to be checked for being equal to 0, which is by default done by converting to a factored normal form. The method of normalization used in LinearAlgebra is configurable, however, as is the method of zero testing, and that can sometimes lead to a solution in a case like this. However, the result, since it is not normalized as much as you would like, is generally not very useful.

In this case, we define C as above and then do:

`Normalizer := x -> x; # No normalization at all.Testzero := testeq; # Fast heuristic zero testing - still extremely reliable.inv := LinearAlgebra:-MatrixInverse(C):`

The computation is now actually pretty fast. But note that I put a colon after that last statement, not a semicolon. That's very much intentional, since Maple can't possibly print any of the entries of the matrix - they are much too big!  For example,

`length(inv[1, 1]);`

returns 1611432722. That is, writing out the [1, 1] entry would require around 10^10 words of memory! The reason Maple can deal with such expressions at all is that it uses a DAG structure for representing expressions, in which repeating subexpressions use memory only once.

So can you use inv at all? Yes, you can. You can, for example, evaluate inv at a point in (x1, ..., y4)-space by substituting numbers in. However, it's not certain that this is actually faster than substituting the values into C and then inverting the matrix numerically. I only did a preliminary test, as follows:

`variables := convert(indets(C), list):mypoint := map(v -> v = RandomTools:-Generate(float(range=0..1, method=uniform)), variables);CodeTools:-Usage(eval(inv, mypoint), iterations=1000):# memory used=0.54MiB, alloc change=0 bytes, cpu time=3.88ms, real time=3.89msCodeTools:-Usage(eval(C, mypoint)^(-1), iterations=1000):# memory used=63.66KiB, alloc change=0 bytes, cpu time=890.00us, real time=892.00us`

This test is unfair, since the matrix inversion happens in a BLAS (in optimized external code) whereas evaluating inv happens in interpreted Maple code. So maybe the evaluation of inv could be improved by setting things up so that you can use evalhf or even Compiler:-Compile, but that would take some work.

All in all, I think in this case you're better off substituting the values first.

Hope this helps,

Erik Postma
Maplesoft.

## Best solution is by transforming to a li...

Hi Andreas,

I entered the data in your example; this looks like data that would be good to fit with an exponential function. (The Excel notation is a little puzzling, but it seems that that's what happens there too.) An exponential function is of course nonlinear, but the fitting process itself is linear: you take the logarithms of all the observations, do a linear fit, and then apply exp to the result. Maple supports this using the ?Statistics/ExponentialFit command:

`with(Statistics):data := <0.07, 2,0.16, 0.24, 0.36, 0.5, 0.8, 1.2, 1.8, 8, 4, 6, 9, 15, 20, 30, 43, 66.9, 100, 150>;PointPlot(data, axis=[mode=log]); # to show that an exponential fit will workresult := ExponentialFit(<seq(1 .. 20)>, data, x);` This assigns result := 0.0866208931860853*exp(0.365097831283158*x). In order to find the R^2 value for this (linear!) fit, we need the residual sum of squares (easy to obtain with ExponentialFit) and the so-called total corrected sum of squares. For this total corrected sum of squares, we can't use ExponentialFit and we'll need to do the transformation by hand. We can obtain the total sum of squares by seeing that the second central moment is the average of the squared deviations from the mean, whereas the total sum of squares is the sum of the squared deviations from the mean - thus they differ by the factor given by the number of data points. Putting this together, we obtain:

`ssRes := Statistics:-ExponentialFit(<seq(1..20)>, data, x, 'output = residualsumofsquares');logdata := log~(data);ssTotal := numelems(logdata) * CentralMoment(logdata, 2);Rsquared := 1 - ssRes / ssTotal;`

This results in a value of 0.9167199079, which is what Excel also seems to get.

If you have a more esoteric fitting function and you need to use ?Statistics/NonlinearFit or, better, ?Statistics/Fit, you can use the same approach to find the pseudo-R^2, except you would take the central moment of data - since there is not necessarily an equivalent of logdata that makes sense in that case.

Hope this helps,

Erik Postma
Maplesoft.

## LinearAlgebra:-NullSpace...

Hi hthomas10,

Do you mean for a matrix θ? If so, you'll want to use ?LinearAlgebra/NullSpace.

Hope this helps,

Erik Postma
Maplesoft.

## Two answers tangential to the rest of th...

Hi @jimmyinhmb,

A very interesting question indeed, and leading to a fascinating discussion. I noticed that one of your questions has not been answered at all, and one only implicitly. First your question 1. To do an in-place sort of an rtable, you can call sort['inplace']. Example:

`restart():v := Statistics:-Sample(Normal(0, 1), 10^7):CodeTools:-Usage(sort['inplace'](v));`

tells me:

`memory used=0.52KiB, alloc change=0 bytes, cpu time=2.05s, real time=2.05s`

Question 3 about ArrayTools:-Alias was answered implicitly already, by acer's use of it: that warning is there because if you really do your best, you can crash Maple's kernel with it. But in order to do that, you need to:

1. manually clear the attributes of the returned rtable, so that it doesn't reference the original rtable anymore;
2. make sure there is no reference to the original rtable anywhere else;
3. cause a garbage collection to occur.

If you're not actually trying to do that, then use of ArrayTools:-Alias is perfectly safe.

Hope this helps,

Erik Postma
Maplesoft.

## Alternative that's efficient for huge ma...

Robert's solution is excellent and probably exactly what the asker wanted. However, I can't help but suggest a version that's substantially more efficient for very big matrices (and, unavoidably, a few lines longer):

`listshuffle := (mat :: Matrix) -> Matrix(op(1, mat), combinat['randperm'](convert(mat, 'list')));statsshuffle := proc(mat :: Matrix)local v :: Vector;  v := ArrayTools:-Alias(mat, [`*`(op(1, mat))]);  v := Statistics:-Shuffle(v);  return ArrayTools:-Alias(v, [op(1, mat)]);end proc;test := Matrix(1000, 1000, symbol=a);CodeTools:-Usage(listshuffle(test));CodeTools:-Usage(statsshuffle(test));`

On my machine, the listshuffle takes 2.89 seconds and 235 MB memory, whereas the statsshuffle takes 0.11 seconds and only 7.7 MB. This is because in statsshuffle:
1) there is only one occasion where the data is copied, and no conversions (ArrayTools:-Alias does not actually convert the data in the rtable that you pass to it, it just puts a different face on it),
2) that one copy occurs in Statistics:-Shuffle, which does its work in fast external code, and works with rtables only.
listshuffle, on the other hand, first converts the big matrix to a list (memory-intensive), then computes a random permutation of the integers in external code (fast), and then permutes the list according to that random permutation (slow, involving a lot of random memory access). It then recreates the Matrix (memory-intensive again).

To wrap this up, if you need to do something like this in a program, consider using tools that do not require you to convert your data type - and see if you can subvert ArrayTools:-Alias to your purposes!

Hope this helps,

Erik Postma
Maplesoft.

PS: I just redid the test with a datatype=float Matrix, and while the timings and memory usage are no different, statsshuffle has the advantage that the Matrix returned automatically has datatype=float. For listshuffle we would need to explicitly tell Maple to reuse the same rtable attributes.

## Wikipedia...

I can't believe I didn't check wikipedia: ﻿http://en.wikipedia.org/wiki/Higher-order_singular_value_decomposition explains that there are two notions of SVD for higher-order tensors, and the one I need is called candecomp-PARAFAC. I guess I'll be implementing that some time next week...

Erik.

## Approximation?...

What a short bike ride can't do for you (from work to home). Here is at least an approximation algorithm, valid for general kth order tensors v.

Let u be the best rank-one approximation to v. If v = u, stop, otherwise approximate v - u recursively.

Now the problem has been reduced to finding the best rank one approximation to a general kth order tensor v. I think we can do this as follows:

Add the entries along one dimension of the tensor, resulting in a (k - 1)th order tensor w. Take the best rank-one tensor approximating w. Find a vector t so that t tensor w is the best possible approximation to v, by considering the components of t individually.

Again a reduction, this time to the same problem for k - 1 (which we solve by induction with SVD as the base case) and to finding a scalar ti  such that ti * w approximates vi best. This sounds like a problem that is manageable; I think it's probably convex.

However, this is certainly only an approximation algorithm: for the tensor tensor minus tensor every resulting approximation is the zero matrix. I guess I will implement this and see how well it performs for images...

Erik Postma
Maplesoft.

﻿