optimizing symbolic code
Schaefer Frank R
Does anybody know how one could possibly use the fact that some
variables are constant
in the optimization of a symbolic formula ?
I have a function that is called extremely often, whereby some variables
are
allways constant and others change. The idea is now to minimize the
computation time
by isolating terms with constant variables. So, constant expressions
are computed
only once, wherelse non-constants are computed in every function call.
Currently, I'm using maples "codegen.,optimize" with the
'tryhard'-option, which really
is a boost. I suppose, however, if there would be an algorithm that uses
the fact
that some expressions are constant for optimization, it would increase
efficiency by at least factor two.
Is there anyone out there who can help ?
Thanks,
Frank.
--
>| Frank Schaefer (Phone: 001-864-656-0482, FAX: 001-864-656-4435)
>| R. 236, Fluor Daniel Building, Mech. Eng.
>| Clemson University, SC 29634, USA
Paul Lutus
news:3A6F0CB7...@ces.clemson.edu...
> Dear Fellows,
>
> Does anybody know how one could possibly use the fact that some
> variables are constant
> in the optimization of a symbolic formula ?
Well, to begin with, "variables that are constant" is a bit confusing. They
are either constants or they are variables.
> I suppose, however, if there would be an algorithm that uses
> the fact
> that some expressions are constant for optimization, it would increase
> efficiency by at least factor two.
Well, obviously an expression consisting entirely of constants can be
converted into a single constant.
If I say --
y = x * (a+(b/c)-(c/d)^e)
-- and a through e are constants, the entire right-hand side can be
converted to a single constant. The basic idea is that constants can be
precomputed -- they don't need to be computed at runtime, because they do
not change at runtime.
--
Paul Lutus
www.arachnoid.com
Frank S.
Paul Lutus wrote:
> y = x * (a+(b/c)-(c/d)^e)
>
> -- and a through e are constants, the entire right-hand side can be
> converted to a single constant....
Ok, you got the basic idea. However, I obviously did not express
myself clearly.
Maple's 'codegen,optimize' program plays arround with different ways of
writing the same mathematical expression in order to optimize
computation
time. It strives to minimize the amount of multiplications, division,
additions,
substractions, subscripts, function calls and assignments.
The problem, now is that 'codegen,optimize' does not know that one
multiplication of two constants is in fact *no multiplication* because
it does only appear during the initialization. So, 'codegen,optimize'
might prefer a writing with 196 multiplications where 5 multiplications
are constant expressions (=> real 191 mults) over a writing with 250
multiplications where
150 multiplications are constant (=> real 100 mults). The final formula
is therefore not
optimal for computation.
So, is there a package like Maple's 'codegen,optimize' that can handle
this ?
Thanks,
Frank.
Paul Lutus
news:3A704775...@ces.clemson.edu...
Well, now that we've dispensed with the basic principles, I cannot add
anything specific because I am not familiar with Maple. In general, try to
arrange things so Maple knows the variables are predefined, that they have
fixed numerical values. This should lead to an elimination process based on
their fixed definitions.
FWIW, Mathematica knows how to do this intrinsically. When an expression is
entered and acted upon, it is converted to its simplest form by default.
--
Paul Lutus
www.arachnoid.com
Richard Fateman
each time through a loop because a and b don't change
within the loop, then a traditional compiler optimization
might notice this and set up something like
temp1:=a*b
....
You want Maple to do this even if the expression
is (perhaps) something like (a+x)*(b+y) which
could be rearranged to
a*b+ a*y +b*y +x*y,
and therefore you could pre-compute a*b.
First of all, note that this is a loser since the expanded
form takes more effort even if you have a*b free.
Secondly, if you are able to notice relatively constant
subexpressions, you can create your own temp computations
outside the loop,
and within the loop you make substitutions of temp for a*b,
etc.
I don't know how clever Maple is in codegen, however.
There are other ways of optimization, e.g. "chains of recurrences"
which may not be part of the standard system.
Good luck
RJF
Richard B. Kreckel
[...]
: FWIW, Mathematica knows how to do this intrinsically. When an expression is
: entered and acted upon, it is converted to its simplest form by default.
Oops, this statement, as it stands, is not correct. Putting aside the
obvious question whether there is such a thing as the "simplest form"
(there isn't) it of course depends on what you mean by "act upon". If
you replace a constant in an expression (this is what we have been
talking about in this thread), the result is not converted to anything
much simpler directly:
In[1]:= a/(a-b+x) - b/(a-b-x) /. x->0
a b
Out[1]= ----- - -----
a - b a - b
In[2]:= Together[%]
Out[2]= 1
Try variants of the above. Quite frankly, I wouldn't want to have a
system that converts my expression to a simpler form "by default" each
time I do something with the expression. I expect the system to do
only what *I* want it to do.
Regards
-richy.
--
Richard Kreckel
<Richard...@GiNaC.DE>
<http://www.ginac.de/~kreckel/>
Paul Lutus
news:94psbl$dp2$1...@bambi.zdv.Uni-Mainz.DE...
> Paul Lutus <nos...@nosite.com> wrote:
> [...]
> : FWIW, Mathematica knows how to do this intrinsically. When an expression
is
> : entered and acted upon, it is converted to its simplest form by default.
>
> Oops, this statement, as it stands, is not correct. Putting aside the
> obvious question whether there is such a thing as the "simplest form"
> (there isn't) it of course depends on what you mean by "act upon". If
> you replace a constant in an expression (this is what we have been
> talking about in this thread), the result is not converted to anything
> much simpler directly:
>
> In[1]:= a/(a-b+x) - b/(a-b-x) /. x->0
I simply meant in the case where symbolic names have been predefined as
constants (assume a,b,c,d,e are initially undefined):
In[10]:= (a-b)/(c+d)^e
Out[10]= (a-b)/(c+d)^e
In[5]:= a=b=c=d=e=5
Out[5]= 5
In[7]:= (a-b)/(c+d)^e
Out[7]= 0
And I should have said "converted to a simpler form by default," rather than
"simplest form." The simplest form achievable in a finite time takes a few
more steps.
> Try variants of the above. Quite frankly, I wouldn't want to have a
> system that converts my expression to a simpler form "by default" each
> time I do something with the expression. I expect the system to do
> only what *I* want it to do.
I can appreciate this, but Mathematica never robs expressions of their
mathematical meaning, and there are ways to prevent the default behavior.
--
Paul Lutus
www.arachnoid.com
Richard Fateman
Paul Lutus wrote:
> <snip>
>
> I can appreciate this, but Mathematica never robs expressions of their
> mathematical meaning, and there are ways to prevent the default behavior.
>
> --
> Paul Lutus
> www.arachnoid.com
I wouldn't go this far. An expression has a meaning to Mathematica.
That may be different from what ordinary humans think of as
the mathematical meaning of that same expression. There are many
ways of showing such flaws, but in particular Mathematica
gets its underwear tied in knots when doing mixed machine-float
and high-precision float arithmetic.
We generally believe that a=b and b=c means a=c. Mathematica 3.0
disagrees, since
1.00000000000000000000000-1/10^18 ==1.0 is True
1.0==1 is True
but
1.00000000000000000000000-1/10^18 ==1 is False
Now starting from this (or similar logical flaws), one
can prove many things that are false.
And we haven't even used any symbols.
Cheers.
RJF
Paul Lutus
news:3A708149...@cs.berkeley.edu...
>
>
> Paul Lutus wrote:
> > <snip>
>
> >
> > I can appreciate this, but Mathematica never robs expressions of their
> > mathematical meaning, and there are ways to prevent the default
behavior.
> >
> > --
> > Paul Lutus
> > www.arachnoid.com
>
>
> I wouldn't go this far. An expression has a meaning to Mathematica.
> That may be different from what ordinary humans think of as
> the mathematical meaning of that same expression. There are many
> ways of showing such flaws, but in particular Mathematica
> gets its underwear tied in knots when doing mixed machine-float
> and high-precision float arithmetic.
I agree completely. Without saying so I was limiting myself to symbolic
expressions, and in any case my original "never" might have been a bit
strong.
Mathematica (like all similar computer programs) is much more robust when
doing symbolic mathematics, and integer-only mathematics, than the
get-our-hands-dirty floating-point variety.
--
Paul Lutus
www.arachnoid.com
Richard Fateman
Paul Lutus wrote:
> Mathematica (like all similar computer programs) is much more robust when
> doing symbolic mathematics, and integer-only mathematics, than the
> get-our-hands-dirty floating-point variety.
>
> --
> Paul Lutus
> www.arachnoid.com
If you are doing polynomial arithmetic on symbols and integers, there
is not much opportunity to mess up. However,
Integrate[1/x,{x,a,b}]
gives an answer -Log[a]+Log[b]
which is clearly a problem for some values of a and b.
For a=-1, b=1, we get -I*Pi.
The implicit premise in Mathematica that it is OK to
give a "generic" answer which is wrong at specific
places represents a divergence from mathematics as
usually practiced.
RJF
Paul Lutus
I agree, and this is certainly a classic. I suppose if I look through the
documentation, I'll find a disclaimer about discontinuities that may appear
"in some cases." :)
The general equation solver generally avoids this sort of error, but it also
provides a straightforward result for y = 1/x without warning about
inappropriate values:
Solve[1/x==y,x]
x == 1/y
In the future, maybe a more sophisticated program will spell out the
caveats.
--
Paul Lutus
www.arachnoid.com
Daniel Lichtblau
>
> "Richard Fateman" <fat...@cs.berkeley.edu> wrote in message
> news:3A709362...@cs.berkeley.edu...
> >
> >
> > Paul Lutus wrote:
> >
> > > Mathematica (like all similar computer programs) is much more robust
> when
> > > doing symbolic mathematics, and integer-only mathematics, than the
> > > get-our-hands-dirty floating-point variety.
I'd be curious to see an example of this non-robustness. I do not claim
that no such examples exist, but rather that in the majority of cases
one sees expected consequences of reasonable design decisions rather
than actual flaws.
On this topic, there is a recent article that discusses Mathematica
numerics for statistics purposes, from an outsiders (that is, non-WRI)
point of view.
Computational Statistics
Volume 15 Issue 2 (2000) pp 279-299
The accurary of Mathematica 4 as a statistical package
B. D. McCullough
Federal Communications Commission, USA,
The article makes some mention of numerical accuracy and general
robustness issues.
> > > --
> > > Paul Lutus
> > > www.arachnoid.com
> >
> > If you are doing polynomial arithmetic on symbols and integers, there
> > is not much opportunity to mess up. However,
> > Integrate[1/x,{x,a,b}]
> >
> > gives an answer -Log[a]+Log[b]
> >
> > which is clearly a problem for some values of a and b.
> > For a=-1, b=1, we get -I*Pi.
Handling of branch cuts for symbolic integration is always problematic;
this example in essence proves that fact. There has been work within the
computer algebra community regarding placing assumptions on the
parameters, and some of this is reflected in our version currently under
development. But I doubt whether anyone has found a generally good way
to handle symbolic branch cut issues.
> > The implicit premise in Mathematica that it is OK to
> > give a "generic" answer which is wrong at specific
> > places represents a divergence from mathematics as
> > usually practiced.
>
> I agree, and this is certainly a classic. I suppose if I look through the
> documentation, I'll find a disclaimer about discontinuities that may appear
> "in some cases." :)
>
> The general equation solver generally avoids this sort of error, but it also
> provides a straightforward result for y = 1/x without warning about
> inappropriate values:
>
> Solve[1/x==y,x]
>
> x == 1/y
Documentation for Solve explicitly notes that solutions are always
generic. If you think a moment about the Solve output design (as
replacement rules, not in the equation form you indicate) you realize it
must be this way: how do you specify e.g. "NotEqual" as a replacement
rule?
In[9]:= Solve[1/x==y,x]
1
Out[9]= {{x -> -}}
y
For nongeneric results one can use Reduce.
In[10]:= Reduce[1/x==y, x]
1
Out[10]= x == - && y != 0
y
> In the future, maybe a more sophisticated program will spell out the
> caveats.
>
> --
> Paul Lutus
> www.arachnoid.com
Daniel Lichtblau
Wolfram Research
Richard Fateman
Daniel Lichtblau wrote:
>
> Paul Lutus wrote:
> >
> > "Richard Fateman" <fat...@cs.berkeley.edu> wrote in message
> > news:3A709362...@cs.berkeley.edu...
> > >
> > >
> > > Paul Lutus wrote:
> > >
> > > > Mathematica (like all similar computer programs) is much more robust
> > when
> > > > doing symbolic mathematics, and integer-only mathematics, than the
> > > > get-our-hands-dirty floating-point variety.
>
> I'd be curious to see an example of this non-robustness. I do not claim
> that no such examples exist, but rather that in the majority of cases
> one sees expected consequences of reasonable design decisions rather
> than actual flaws.
>
Hi Dan!
I'd like to invite any readers to see my (now 11+ years old) review of
Mathematica at
http://www.cs.berkeley.edu/~fateman/papers/mma.review.pdf
I have heard that people at WRI claim that "all the bugs have been fixed"
that were reported in that review. Since my intention in 1990 was NOT to
mention bugs that were easily fixed (and also so common), but to mention
design flaws, I looked forward to new releases of Mathematica.
While some of the details have been moved around, in my view, nearly all
the flaws remain.
Many of them are hardly expected consequences to the users of mathematica,
although I agree with Dan that one could say that someone truly familiar with the design
decisions might be able to anticipate the bizarre consequences. That's
how I (and W. Kahan) came up with these examples.
Here's an example from that paper.
p=314159265358979323;
q=314159265358979323.;
r=314159265358979323.00000000000000000000;
s=p+0.00000000000000000000;
{Tan[s],N[Tan[p]],Tan[q],Tan[r]}
The result is
{1.59981,1.59981,-0.,-1.1297926523089085443}
in the Mathematica I used for the review, I got ComplexInfinity instead of -0.
{p==q, q==r, r==s, r==p}
The result is
{True,True,True,True}
> On this topic, there is a recent article that discusses Mathematica
> numerics for statistics purposes, from an outsiders (that is, non-WRI)
> point of view.
>
> Computational Statistics
> Volume 15 Issue 2 (2000) pp 279-299
> The accurary of Mathematica 4 as a statistical package
> B. D. McCullough
> Federal Communications Commission, USA,
>
> The article makes some mention of numerical accuracy and general
> robustness issues.
I do not have access on-line to this paper, but it refers to another
paper on a similar topic (using EXCEL for statistics) by the same
author. It seems that EXCEL and many other packages do not use
the best known numerical techniques. This may be attributed to
ignorance (e.g. choosing a bad pseudo-random number generator)
age, speed tradeoffs, etc. I have no doubt that Mathematica's
choice of floating-point library routines could be more astute.
I also have no doubt that --given sufficient programming effort--
Mathematica's big-float capability can be used to compute high
accuracy results. I DO however, believe that naive users
-- who do not
make a study of the Mathematica design, but merely use it as
though it were doing mathematics as taught in numerical analysis
courses-- will be unpleasantly surprised on occasion. Or
if they are not surprised, they may just be getting wrong
answers.
I mostly agree with Dan's comments except that the idea that
it is OK to do something wrong if it is documented (i.e. in
the Solve command description) should be discouraged. It
is not ok to convert a bug into a feature by documenting it!
RJF
Paul Lutus
news:3A71D62F...@wolfram.com...
> Paul Lutus wrote:
> >
> > "Richard Fateman" <fat...@cs.berkeley.edu> wrote in message
> > news:3A709362...@cs.berkeley.edu...
> > >
> > >
> > > Paul Lutus wrote:
> > >
> > > > Mathematica (like all similar computer programs) is much more robust
> > when
> > > > doing symbolic mathematics, and integer-only mathematics, than the
> > > > get-our-hands-dirty floating-point variety.
>
> I'd be curious to see an example of this non-robustness. I do not claim
> that no such examples exist, but rather that in the majority of cases
> one sees expected consequences of reasonable design decisions rather
> than actual flaws.
I agree completely. I meant "robust" in a more basic sense, not in a
disparaging sense at all. The same reasoning could be applied to a skilled
mathematician confronted by a problem expressed entirely symbolically and/or
with integers vs. an "equivalent" floating-point expression. There are some
innate limitations in many FP expressions.
But, since we're on the topic, here is an example from another thread:
Compare the outcome for:
(-1)^(0.5)
vs.:
(-1)^(1/2)
I know *why* this happens, but it is disconcerting nevertheless.
--
Paul Lutus
www.arachnoid.com
Richard Fateman
Paul Lutus wrote:
>
>.....
>
> But, since we're on the topic, here is an example from another thread:
>
> Compare the outcome for:
>
> (-1)^(0.5)
>
> vs.:
>
> (-1)^(1/2)
>
> I know *why* this happens, but it is disconcerting nevertheless.
I think any time a computer algebra system surprises people it is
an indication of a potential problem. If the surprise cannot be
explained
by mathematics, but has to resort to "this is the way we did it",
then it is an actual problem. (A good explanation is
"the system does it this way because it is consistent with the
larger body of mathematics that you may not yet know about"... e.g.
why freshman calculus answers involving log(abs(...)) are not
produced... the system knows about complex numbers!)
There is nothing inevitable about Mathematica's bug on this problem.
(-1)^0.5 coming out 6.12..x10^-17+ 1. I.
In fact Macsyma returns this expression unchanged.
If you insist on numerically evaluating it, you get 1.0*I.
There is no necessity for 0.5, which is representable exactly in
binary floating point, to produce any error. It is certainly NOT
dictated by any mathematics.
There is an issue about whether the answer should include +I and -I,
in my view.
RJF
Richard B. Kreckel
[...]
: On this topic, there is a recent article that discusses Mathematica
: numerics for statistics purposes, from an outsiders (that is, non-WRI)
: point of view.
:
: Computational Statistics
: Volume 15 Issue 2 (2000) pp 279-299
: The accurary of Mathematica 4 as a statistical package
: B. D. McCullough
: Federal Communications Commission, USA,
:
: The article makes some mention of numerical accuracy and general
: robustness issues.
Now this makes me curious. I always felt Mathematica's precision
arithmetic is slightly ridiculous -- it tries to be more clever than
the programmer at points where it really should just compute a number
and let the user care about precision. The attempt is of course
well-meant but I always found it makes the product unuseful for my
purposes. :-(
Isn't there an online version of that article around that we can all
have a look at? I remember when you publish something at a Springer
Journal they give you an DOI (digital object identifier), basically a
link to your article in PDF. In this case, I can't find one.
Regards
-richy.
--
Richard Kreckel
<Richard...@Uni-Mainz.DE>
<http://wwwthep.physik.uni-mainz.de/~kreckel/>
Paul Lutus
> Daniel Lichtblau <da...@wolfram.com> wrote:
> [...]
> : On this topic, there is a recent article that discusses Mathematica
> : numerics for statistics purposes, from an outsiders (that is, non-WRI)
> : point of view.
> :
> : Computational Statistics
> : Volume 15 Issue 2 (2000) pp 279-299
> : The accurary of Mathematica 4 as a statistical package
> : B. D. McCullough
> : Federal Communications Commission, USA,
> :
> : The article makes some mention of numerical accuracy and general
> : robustness issues.
>
> Now this makes me curious. I always felt Mathematica's precision
> arithmetic is slightly ridiculous -- it tries to be more clever than
> the programmer at points where it really should just compute a number
> and let the user care about precision. The attempt is of course
> well-meant but I always found it makes the product unuseful for my
> purposes. :-(
I assume you mean Mathematica's convention of only providing as many decimal
places as the original numbers warrant? Like:
Sqrt[2]^2 = 2, but
1.7`2/3 = 0.57, and
1.7`8/3 = 0.56666667
I think this is a nice feature. It contradicts a human failing -- to assume
that all the digits in a result have meaning.
--
Paul Lutus
www.arachnoid.com
Daniel Lichtblau
Our overactive news server deleted alot of replies that were not even
three days old, so I am digging stuff out of the deja archives. I'll put
it all into one message, quoting previous posters as indicated. Opinions
herein are my own and do not necessarily reflect those of my employer.
--------------------
[Richard Fateman:]
I'd like to invite any readers to see my (now 11+ years old) review of
Mathematica at
http://www.cs.berkeley.edu/~fateman/papers/mma.review.pdf
I have heard that people at WRI claim that "all the bugs have been
fixed"
that were reported in that review. Since my
intention in 1990 was NOT to mention bugs that were easily fixed (and
also
so common), but to mention design flaws, I looked forward to new
releases
of Mathematica. While some of the details have been moved around, in
my
view, nearly all the flaws remain.
Many of them are hardly expected consequences to the users of
mathematica,
although I agree with Dan that one could say that someone truly
familiar
with the design decisions might be able to anticipate the bizarre
consequences.
That's how I (and W. Kahan) came up with these examples.
p = 314159265358979323;
q = 314159265358979323.;
r = 314159265358979323.00000000000000000000;
s = p+0.00000000000000000000;
{Tan[s],N[Tan[p]],Tan[q],Tan[r]}
[...]
------------------
[my response:]
Here are the results in version 4.1:
In[9]:= {Tan[s],N[Tan[p]],Tan[q],Tan[r]}
Out[9]= {1.60125, 1.60125, ComplexInfinity, -1.1297926523089085443}
Since s is a machine number, Tan[s] and N[Tan[p]] should
(and do) give the same result. The fact that that number has no correct
digits is by no means wrong. It is a consequence of a decision to use
machine arithmetic in evaluating elementary functions at machine
number arguments. I certainly believe this is the best thing to do.
As you know, there is a body of literature related to this, on validated
predicate testing ("at least give me the correct sign!") e.g. for
computational geometry. And often one must resort to hard labor to avoid
exactly the pitfalls of naive machine arithmetic. My point is that
Mathematica users (and implementors) are not magically immunized to the
harsh reality that this sort of thing is not foolproof.
We have r to high precision, so Tan works with it at that precision, and
gives
the "correct" result. We only know q to low precision as a Mathematica
"bignum". Hence after reducing modulo Pi we "think" we have zero. I
suppose
this may be a surprise to a naive user, but of course that does not make
it
wrong.
--------------------
[Richard Fateman, cont'd:]
I DO however, believe that naive users -- who do not
make a study of the Mathematica design, but merely use it as
though it were doing mathematics as taught in numerical
analysis courses-- will be unpleasantly surprised on occasion.
[...]
--------------------
[response, cont'd:]
Certainly. This even happens to sophisticated users from time to time.
Sometimes this indicates problematic design but far more often it
is simply misunderstanding about reasonable design decisions.
The fact that {p,q,r,s} are all deemed equal to one another is a
consequence of equality testing in Mathematica (in short: coerce to
lowest precision and test). I am not sure if it is under dispute in
this instance so I'll only note that I believe it is the most useful way
to do this. And I acknowledge that there may be room for disagreement.
---------------------
[Richard Fateman, cont'd:]
I mostly agree with Dan's comments except that the idea that
it is OK to do something wrong if it is documented (i.e. in
the Solve command description) should be discouraged. It
is not ok to convert a bug into a feature by documenting it!
---------------------
[response, cont'd:]
Not sure I followed this, but if you are claiming that Solve[1/y==x,y]
is "wrong" then you are wrong (and as every kid knows, the fact that
I used two wrongs makes me right). My opinion is that a generic Solve
is an important tool and should not be held hostage to examples with
measure zero "bad" sets. We have Reduce[] to handle those.
---------------------
[Paul Lutus:]
But, since we're on the topic, here is an example from another thread:
Compare the outcome for:
(-1)^(0.5)
vs.:
(-1)^(1/2)
I know *why* this happens, but it is disconcerting nevertheless.
--------------------
[Richard Fateman, in response:]
There is nothing inevitable about Mathematica's bug on this problem.
(-1)^0.5 coming out 6.12..x10^-17+ 1. I.
In fact Macsyma returns this expression unchanged.
If you insist on numerically evaluating it, you get 1.0*I.
There is no necessity for 0.5, which is representable exactly in
binary floating point, to produce any error. It is certainly
NOT dictated by any mathematics.
--------------------
[My response:]
This one involves some subtleties and I suspect there is no perfect
handling. I'll discuss some of the issues. For one, is the "correct"
result
still correct after the problem is perturbed? Specifically, if you get a
result of 1.0*I, what does that same program give for
(-1)^0.500000000001?
Or for (-1)^0.4999999999999? If you do not get real parts, then the
software that makes (-1)^(.5) purely imaginary is not behaving
gracefully.
I'll offer some opinions about (-1)^(.5). First, if one numerically
evaluates by finding exp(log(-1)*(.5)), a not unreasonable thing to do,
then the spurious real part is inevitable. While 0.5 is exactly
representable
as a binary float this is most certainly not the case for 0.5*Pi.
Suppose
instead one makes special cases for negatives raised to floats? Is this
a useful thing to do? From the point of view of avoiding a spurious real
part, no, not really. The reason is that this only helps for half
integer
powers. So these are indeed quite special.
Is it really so useful to avoid such spurious real parts? I think this
depends alot on specific needs. For example, it is far less important
than
recognizing spurious imaginary parts e.g. in numeric eigenvalues. (I
have
a bias here. I rely on numeric eigen-methods to do some root finding in
NSolve of polynomial systems, and people want real roots to be
explicitly
real.) On the other hand, avoiding spurious real parts can certainly
play a role in correct handling of branch cuts. This, however, is not
unlike the issues raised in the tangent evaluations at the top: near
discontinuities, getting evaluations "correct" is nontrivial and often
cannot be done with machine arithmetic. Moreover, this issue goes well
beyond getting square roots of reals to lie on the imaginary axes. For
example, how might you correctly handle cube roots?
On the topic of NSolve...
--------------------
[Richard Kreckel:]
Now this makes me curious. I always felt Mathematica's precision
arithmetic is slightly ridiculous -- it tries to be more clever than
the programmer at points where it really should just compute a number
and let the user care about precision.
[...]
--------------------
[My response:]
The technology used by NSolve often involves computation of an
approximate Groebner basis, which in turn uses Mathematica's
significance
arithmetic. While there are other ways to achieve NSolve functionality,
this approach has worked out fairly well. It would be MUCH harder
without
that precision control mechanism. This is but a particular example
of the utility of that arithmetic. Another would be in a topic I touched
upon above, validated predicate testing for computational geometry.
Again,
there may be other ways to achieve the same functionality, but that does
not diminish the importance of significance arithmetic as a tool for
such purposes.
--------------------
As a concluding remark, I will point out that the examples posted all
exhibit some form or another of pathology. Such examples can of course
demonstrate actual software bugs but more often serve to illuminate
design decisions and/or limitations intrinsic to numerical computation.
This is also quite important, I just thought I'd point out the
distinction.
Daniel Lichtblau
Wolfram Research
CINDY SMITH
Daniel Lichtblau <da...@wolfram.com> writes:
> [Paul Lutus:]
> But, since we're on the topic, here is an example from another thread:
>
> Compare the outcome for:
>
> (-1)^(0.5)
>
> vs.:
>
> (-1)^(1/2)
I did this in MathCad 8 and got the answer "i" in both cases.
Hmmmm.....don't tell me, it isn't a bug, it's a feature..... :-)
> I know *why* this happens, but it is disconcerting nevertheless.
> [Richard Fateman, in response:]
> Daniel Lichtblau
> Wolfram Research
--
Cindy Smith oooooooooo O oooooooooo O oooo
c...@dragon.com ooo O () o
c...@5sc.net O oooooooooo O oooooooooo O ooooo
|
A Real Live --+-- "Whoever misgoverns a house inherits the wind,
Catholic in | and the fool becomes slave to the wise."
Georgia! | NJB Proverbs 11:29
Spawn of a Jewish Carpenter ><> IXQYS <>< Go against the flow!
Paul Lutus
> > But, since we're on the topic, here is an example from another thread:
> >
> > Compare the outcome for:
> >
> > (-1)^(0.5)
> >
> > vs.:
> >
> > (-1)^(1/2)
>
> I did this in MathCad 8 and got the answer "i" in both cases.
>
> Hmmmm.....don't tell me, it isn't a bug, it's a feature..... :-)
Well, there HAD to be something MathCad does better than Mathematica. :)
--
Paul Lutus
www.arachnoid.com
Dima Pasechnik
[...]
> [Paul Lutus:]
> But, since we're on the topic, here is an example from another thread:
>
> Compare the outcome for:
>
> (-1)^(0.5)
>
> vs.:
>
> (-1)^(1/2)
>
> I know *why* this happens, but it is disconcerting nevertheless.
>
> --------------------
>
> [Richard Fateman, in response:]
> There is nothing inevitable about Mathematica's bug on this problem.
> (-1)^0.5 coming out 6.12..x10^-17+ 1. I.
> In fact Macsyma returns this expression unchanged.
> If you insist on numerically evaluating it, you get 1.0*I.
>
> There is no necessity for 0.5, which is representable exactly in
> binary floating point, to produce any error. It is certainly
> NOT dictated by any mathematics.
>
> --------------------
>
> [My response:]
>
> This one involves some subtleties and I suspect there is no perfect
> handling. I'll discuss some of the issues. For one, is the "correct"
> result
> still correct after the problem is perturbed? Specifically, if you get a
> result of 1.0*I, what does that same program give for
> (-1)^0.500000000001?
Please define "the same program".
Same source code on the same hardware?
Same source code on different hardware?
Does Mathematica hold the user hostage of the particular hardware
used at the moment ?
[...]
>
> On the topic of NSolve...
>
> --------------------
>
> [Richard Kreckel:]
>
>
> Now this makes me curious. I always felt Mathematica's precision
> arithmetic is slightly ridiculous -- it tries to be more clever than
> the programmer at points where it really should just compute a number
> and let the user care about precision.
>
> [...]
>
> --------------------
>
> [My response:]
>
> The technology used by NSolve often involves computation of an
> approximate Groebner basis, which in turn uses Mathematica's
> significance
> arithmetic. While there are other ways to achieve NSolve functionality,
> this approach has worked out fairly well. It would be MUCH harder
> without
> that precision control mechanism.
Hmm, what do the inner guts of NSolve have to do with the casual
use of the system?
Something that might work well for a developer
of the system need not be OK for a user of the system.
[...]
> Daniel Lichtblau
> Wolfram Research
Dima Pasechnik
http://ssor.twi.tudelft.nl/~dima/
Daniel Lichtblau
Good point. I had in mind same source, arbitrary hardware/OS. Now that
you mention it, that may have been naive. Mathematica uses OS (and, by
extension, hardware) dependent routines for machine arithmetic. So I
think the answer is that the outcome of such computation can be machine
dependent.
Regardless of whether or not (-1)^(.5) gives machine or OS dependent
results, there are machine arithmetic computations that will exhibit
such dependencies in Mathematica. I do not regard this as a bug or
design flaw (quite the contrary), but others might disagree. Does it
really amount to being "held hostage"?
My sense is that alot of OS-dependent anomalous behavior in mathematical
and related software has helped to push hardware/OS vendors toward IEEE
compliance. This of course has cut back on the more extreme examples of
hardware dependent results. But I doubt they are entirely gone.
> [...]
> >
> > On the topic of NSolve...
> >
> > --------------------
> >
> > [Richard Kreckel:]
> >
> >
> > Now this makes me curious. I always felt Mathematica's precision
> > arithmetic is slightly ridiculous -- it tries to be more clever than
> > the programmer at points where it really should just compute a number
> > and let the user care about precision.
> >
> > [...]
> >
> > --------------------
> >
> > [My response:]
> >
> > The technology used by NSolve often involves computation of an
> > approximate Groebner basis, which in turn uses Mathematica's
> > significance
> > arithmetic. While there are other ways to achieve NSolve functionality,
> > this approach has worked out fairly well. It would be MUCH harder
> > without
> > that precision control mechanism.
> Hmm, what do the inner guts of NSolve have to do with the casual
> use of the system?
> Something that might work well for a developer
> of the system need not be OK for a user of the system.
The question was, roughly, "Why might one believe Mathematica's
significance arithmetic is useful"? My answer was roughly, "It is quite
useful to NSolve". Indeed this is arcane, and most casual use will not
require this sort of functionality. But certainly some does. Moreover I
do not think the original question really had casual use in mind: the
phraseology "let the user care about precision" indicates to me a
far-from-casual user. Hence an example involving functionality that is
difficult to achieve without significance arithmetic is reasonable.
(Have you any idea how much harder it is to obtain reliable approximate
Groebner bases using fixed precision arithmetic?)
I confess to some confusion. I do not really understand what exactly is
the issue you wish to raise. Are you saying that significance is not
useful because most users do not directly use it? That would hardly be a
reason not to support it in a high-end program with numerous
sophisticated users in additional to the casual ones.
There is, moreover, indirect use, and this does affect casual users. Let
me give an answer that I probably should have given the first time. Low
precision significance arithmetic is used to handle numeric predicate
testing e.g. in
In[42]:= E^Pi < Pi^E
Out[42]= False
Certainly there are predicates of this sort that one cannot successfully
evaluate, and these simply return unevaluated. But many can be done this
way, and much more easily than when using fixed precision arithmetic.
Daniel Lichtblau
Wolfram Research
Richard Fateman
of an approximate Grobner basis to actual solutions? I thought
I was following the approximate GCD stuff, but maybe I've missed
this.
In any case, the idea that one can more easily compute something
using Mathematica's significance arithmetic vs. without it could
be supported by an attempt to program something without it vs.
with it.
Clearly the tools needed ... arbitrary-length "bigfloats" are
available in just about every computer algebra system... In fact, I know
Macsyma has a predicate testing program on constants and presumably
does NOT use significance arithmetic for it, though it will crank
up the precision as necessary (or, perhaps, give up without a decision
because the two sides seem to be the same but can't be proved to
be the same).
RJF
Richard Fateman
Maybe another example would help.
Consider
test= N[10^30,20]
this is 10^30 to 20 decimal places in Mathematica.
compute logtest= Log[test]. You get 69.0775527898213705205
Mathematica thinks it is accurate to 20 decimal places.
Now forceably convert it to 35 places :
SetPrecision[logtest]
69.07755278982137052053974364053093
All these digits are correct, too. So Mathematica's
intuition here is wrong.
As most of us know, Log changes rather slowly. Mathematica
doesn't know this, even though it could in principle figure
this out by taking a derivative. So a carefully
written program would likely have to do this:
(a) compute the proper error bounds, perhaps based on
other formulas, derivatives, etc.
(b) disabling/resetting the precision and/or accuracy
of some (all?) default numerical results from mathematica
to the correct ones.
Will Mathematica tell you that more digits are correct
than are justified? I haven't fiddled enough with the
latest Mathematica to find any examples, and it could
be that it is now sufficiently conservative that this
phenomenon is gone.
RJF
Daniel Lichtblau
>
> I'm not sure what NSolve is doing; what is the relationship
> of an approximate Grobner basis to actual solutions? I thought
> I was following the approximate GCD stuff, but maybe I've missed
> this.
An approximate GB can be used to cough up an endomorphism matrix
representing e.g. multiplication by x1 in the zero dimensional algebra
C[x1,...,xn]/<polynomial system> (assuming this algebra is in fact zero
dimensional). This all follows from a 1988 Auzinger & Stetter article
and independent work of others shortly thereafter. Eigenvalues of this
matrix give values of x1 at all roots of that system, and more generally
there are eigen-methods to get all solutions in a reasonably efficient
and numerically accurate manner. There is a body of literature on this
technology though I think generally either exact Groebner bases or other
methods entirely are used to get that matrix.
Gauging the accuracy of the roots may be tricky, but really no more so
than is the case for a Newton-type root finder.
> In any case, the idea that one can more easily compute something
> using Mathematica's significance arithmetic vs. without it could
> be supported by an attempt to program something without it vs.
> with it.
> Clearly the tools needed ... arbitrary-length "bigfloats" are
> available in just about every computer algebra system... In fact, I know
> Macsyma has a predicate testing program on constants and presumably
> does NOT use significance arithmetic for it, though it will crank
> up the precision as necessary (or, perhaps, give up without a decision
> because the two sides seem to be the same but can't be proved to
> be the same).
> RJF
> [...]
How does it handle catastrophic cancellation? If it ignores this issue,
it is likely to give (occasionally) erroneous results when something
equivalent to zero appears to be nonzero and is deemed as such. If it
pays careful attention to this issue then it is most likely emulating
significance arithmetic at least to a modest extent.
Daniel Lichtblau
Wolfram Research