I don't have better code at the moment, but I can make a few comments.

The way that Matlab used to compute the Rank (and perhaps it still does) is to look at the ratio of the singular values to the largest singular value. The rank then being the number of sing. values for which the ratio satisfies a condition involving a (working- or double- precision) machine constant.

Maple also does the above, by first reducing to bidiagonal form (clapack dgebrd or nag f08kef), and then computing all singular values (clapack dbdsqr or nag f08mef), and then testing the ratios.

How could it be improved? For the both stages, it might be possible to  use a faster implementation of those named routines, from, say MKL. That is not something that could be tried with recompiling part of Maple, as the bundled MKL does not include the static .lib's as per usual Intel licensing of that product. Another possibility is the Eigen package, which is a quite separate implementation of floating-point linear algebra routines.

Another possibility might be to use a so-called rank-revealing QR computation.

Or, if not using MKL then the "optimal" blocksize as returned by (clapack) ilaenv (or nag equivalent) for the bidiagonal form reduction stage's performance could be checked and possibly tweaked.

Or, for the second stage it might be possible to use another routine, or another approach altogether. There are some notes which suggest that the newer divide and conquer routines are faster in the case that singular vectors are computed, but that otherwise those drivers use the same `dqds` algorithm as before. Are there other ways? At the expense of memory copies is there a possibility to compute only "selected" singular values for the bidiagonal form, so that the second stage might be be run in parallel (by splitting that second task in half, etc)? I don't know.

I have in fact been investigating this very topic, as I had noticed the timings posted at the URL you cited and have long wanted to look at the newer routines for computing singular vectors. (I might have seen it on usenet a short while back.) If I come up with an improvement on the Rank timing then I'll let you know.

I have a recollection that at one time in the recent past the URL you cited showed that Matlab was by far the slowest of the three M's on the rank calculation. Either I remember it wrongly, or they have improved their approach in the near past. It is very interesting (and likely no coincidence I suspect) how closely the Mma and Matlab timings you've cited now match each other across problem size. Also interesting is the performance ratio as compared to Maple 17, as the problem size gets large.

On a related note, I do have a simple proc that does wrapperless external call to the faster divide & conquer (and slower but more accurate jacobi method) drivers for the more expensive job of computing either all or only min(m,n) singular vectors. Perhaps I can find time for a post soon-ish.

acer

@marc005 That restriction on the Matrix brower appears to be present in 17. But it's not difficult to write custom procedures which do similar things.

For examination of the entries by value, you could just use a Data (table) Component. If the Matrix already is assigned to a name in the session, then inserting such an embedded component from the Components palette will very nicely prompt whether you wish to use it.

The attached worksheet took just a few moments to write. A little longer and it could be spiffier still.

quandlfun.mw

@marc005 The message "Unable to display structured data" arises from the Matrix browser (also available via right-click context-menu), when either dimension exceeds a certain value. It's not a consequence of that `GetCSV` procedure, per se.

The limiting value for either dimension appears to be 10001.

The Matrix browser has a few quirks. For example, the `image` view does not seem to work when either row or column dimension is just 1.

@Carl Love Yes, I get agreement (except on that one mentioned subexample). When an assumption is placed, including as additional qualification, a new assumed local is generated and assigned to the global name. It's always been problematic (a FAQ of sorts) that people get into trouble if they make assignments using some intermediate stage assumptions and then combine expressions involving the mix of assumed names. I don't know whether the word "intended" connotes to everyone, but that is the longstanding implementation and, presumably, the early design.

I find it less error prone (and easier to keep straight) to use `assuming` rather than `assume`, and this kind of thing is one reason.

@Carl Love Yes, I get agreement (except on that one mentioned subexample). When an assumption is placed, including as additional qualification, a new assumed local is generated and assigned to the global name. It's always been problematic (a FAQ of sorts) that people get into trouble if they make assignments using some intermediate stage assumptions and then combine expressions involving the mix of assumed names. I don't know whether the word "intended" connotes to everyone, but that is the longstanding implementation and, presumably, the early design.

I find it less error prone (and easier to keep straight) to use `assuming` rather than `assume`, and this kind of thing is one reason.

So your y1,...yn are floats? If so then what happens if you raise Digits in Maple? (Ie. to 15, and then to 30?) Numeric roundoff error can be different for different (equivalent) forms of an expression, even without the kind of algebraic cancellation that Carl described. Applying `simplify` will not necessarily reduce the numeric roundoff effects for subsequent floating-point substituion and evaluation.

Why use `subs` instead of `eval`, btw?

acer

So your y1,...yn are floats? If so then what happens if you raise Digits in Maple? (Ie. to 15, and then to 30?) Numeric roundoff error can be different for different (equivalent) forms of an expression, even without the kind of algebraic cancellation that Carl described. Applying `simplify` will not necessarily reduce the numeric roundoff effects for subsequent floating-point substituion and evaluation.

Why use `subs` instead of `eval`, btw?

acer

@Carl Love Using your code above I get,

  > about(x1);

Originally x, renamed x~:
  is assumed to be: RealRange(Open(1),infinity)

even after the additional qualification `x>2` is made.

I get that in Maple 11 through 17, on Windows.

Above your results show,

  > about(x1);

Originally x, renamed x~:
  is assumed to be: RealRange(Open(2),infinity)

which I don't obtain.

@Carl Love Using your code above I get,

  > about(x1);

Originally x, renamed x~:
  is assumed to be: RealRange(Open(1),infinity)

even after the additional qualification `x>2` is made.

I get that in Maple 11 through 17, on Windows.

Above your results show,

  > about(x1);

Originally x, renamed x~:
  is assumed to be: RealRange(Open(2),infinity)

which I don't obtain.

Is your acquaintance running 32bit or 64bit Maple on Windows?

acer

@Carl Love Carl, you wrote in a Comment above that, "Most interesting is that all copies end up with the same assumptions."

In which version did you obtain those results from the `about` command, may I ask?

@Carl Love Carl, you wrote in a Comment above that, "Most interesting is that all copies end up with the same assumptions."

In which version did you obtain those results from the `about` command, may I ask?

@emma hassan Well, you already have a table, as V;

eval(V);

If you want you can turn V into a Matrix (and then right-click and browse as one way to view its entries), or call plots:-surfdata on that.

matdata:=Matrix(M,N,V);

plots:-surfdata(matdata,axes=box);

@emma hassan Well, you already have a table, as V;

eval(V);

If you want you can turn V into a Matrix (and then right-click and browse as one way to view its entries), or call plots:-surfdata on that.

matdata:=Matrix(M,N,V);

plots:-surfdata(matdata,axes=box);

Are a,b,c,... each something numeric? Are F1 and F2 numeric? Are they floating-point?

acer