Quick summary: Maybe there needs to be a numerics-sig before any
of these changes are adopted...
And Huaiyu Zhu writes:
-
- (4) Similarity with Matlab (but not really compatibility).
Take care here. Matlab has made many arbitrary choices that kill
consistency. Their treatment of 1x1 `matrices' has been discussed
to death in the NA-digest. And the lack of algorithm specification
for matrix /, \, etc. has caused more than a few surprises.
- Notation: To avoid confusion, I'll adopt matlab notation and invent
- a set of unique names in the following discussions:
I'll trade all the fuss over notation for fuss over _definitions_.
I'm accustomed to using xGEMM, xGESVX, etc. in any language. Why
those rather than whatever the language throws at me? They're
well-defined. I know their error bounds and behavior in the face
of floating-point oddities. Give me those definitions, and I'll
take almost any notation. Hell, I've been putting up with the
BLAS and LAPACK naming scheme.
Before anyone adds new notations for operations, you need to be
sure you have definitions for those operations. And good definitions
will exacerbate some of the numeric issues in Python. For example,
Python doesn't provide reasonable ways to test for Inf or NaN. These
are needed in new eigenvalue routines (the only known optimal ones).
Yes, I know, there are some non-IEEE architectures left, but they
shouldn't be encouraged. (Keep in mind that I'm at Berkeley; expect
a strong IEEE-754 bias.)
And then there are issues with complex. -0.0j = -0j, but
complex(0.0,-0.0) = 0j. And there seems to be no way to construct
1-0j at all. This is painful when you think about branch cuts
(http://mathworld.wolfram.com/BranchCut.html) combined with
conditional tests. The typical graphical example shows an airfoil
with flow going in a completely wrong path, but I can't find that
on-line right now.
There are standard definitions for most matrix and vector operations
in the BLAS and LAPACK. The next BLAS standard is being drafted now,
see
http://www.netlib.org/cgi-bin/checkout/blast/blast.pl
The authors are trying to simplify the current text as much as
possible, so it may be a while before it's fully released. The LAPACK
and BLAS definitions are much more consistent than the Matlab ones.
Some extra reading on numeric issues in some languages:
Some Java-related:
http://www.validgh.com/java/
(Of course, I recommend Dr. Kahan's slides, available there)
Some Matlab-related:
http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
http://www.cs.berkeley.edu/~wkahan/Cantilever.pdf
http://www.netlib.org/na-digest-html/99/v99n43.html#2
(and previous issues of na-digest)
Great IEEE754 feature rationale:
http://www.cs.berkeley.edu/~wkahan/ieee754status/ieee754.ps
- - How to deal with left division?
How to compute it? Iterative or direct? With what pivoting
strategy if direct, or preconditioning if iterative? For which
matrix sizes? What about symmetric matrices? Treating them as
unsymmetric adds a _lot_ of error.
- - How to deal with power? It is already two characters
How to define power? In the matrix^scalar case, are you going to
use the eigenvalues or calculate it explicitly? There are huge
differences in both numeric behavior and performance. Which to
use depends on matrix size...
Matrix^matrix does exist, but it's terribly difficult to calculate
and almost never used. (I think. Could be wrong.)
- - Do we want a copy operator := which is different from =? It allows
- efficient code. The following code does not change A
- B := A
- B += C
Neither does B = A + C. And := is already a horribly overloaded
symbol in the literature. Sometimes it means definition (aliasing,
referencing), sometimes it means assignment (copying). Blah.
- This definitely will be usable for other classes. Call it __copy__?
Deep or shallow? Letting it up to the class involved makes it
unpredictable, so it shouldn't be an operator.
- These are the things that come to mind. Now, opinions, please.
Please don't rush to a syntax before you have a semantics. Sorry for
the length of this, but all I've seen to far is discussion over what
syntax to use for various undefined operations. That bothers me.
It's how Matlab started, and now they're trying to drag it back
out of its craziness. Getting hit by an odd numerical bug doesn't
happen often, but when it does things blow up or get recalled in
huge numbers... Seriously, debuging numerical problems is painful
even for most numerical analysts. It's better to try to prevent
them.
There's the possibility of Python being _better_ than Matlab, Octave,
SciLab, etc. and not just being equal to them.
Jason