You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@sis.apache.org by Martin Desruisseaux <ma...@geomatys.fr> on 2013/09/11 23:36:40 UTC
Report on matrix work
Hello all
The javadoc package documentation below tries to explain the rational
behind the SIS matrix package (still work in progress):
https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/package-summary.html
The idea is strongly inspired by "vecmath", which was part of Java3D.
"Vecmath" is actually hard to replace by other matrix package (JBLAS,
JAMA or other) because the design goals are not the same. There is
thousands of possible operations for matrix, and the operations of
interest are not the same for resolving linear equations than for
performing coordinate transformations. Because "vecmath" was designed
for coordinate conversions in 3D spaces, it was a pretty good fit to
geospatial needs.
There is a few methods in org.apache.sis.referencing.operation.matrix
that are more closely related to coordinate transformations rather than
general matrix operations: isAffine(), toAffineTransform(),
normalizeColumns() - more will come soon. I have not yet connected the
SIS package to JAMA. But I came to think that the connection (if any)
should not be visible in the public API. For example when inverting a
non-square matrix, JAMA uses a least squares approach. This is probably
a good choice for resolving linear equations, but sis-referencing needs
a more conservative approach. For this reason, I don't think that we
should expose the JAMA Matrix.inverse() method directly.
I will post more when the package will be closer to completion, unless
someone would like clarifications or changes?
Martin
Re: Report on matrix work (continuing)
Posted by Adam Estrada <es...@gmail.com>.
Outstanding! Thanks Martin!
A
On Fri, Oct 4, 2013 at 4:31 PM, Martin Desruisseaux <
martin.desruisseaux@geomatys.fr> wrote:
> Hello all
>
> I just completed the work on extended floating point precision in matrix
> multiplication and inversion. This is new work compared to what we had in
> Geotk. I added in the package javadoc [1] a small section, namely "Extended
> floating point precision", which try to explain why we care about the
> accuracy of those operations.
>
> Now the last work needed is to port existing Geotk code for handling of
> NaN values and some special kind of non-square matrix, and I should be done
> with matrix. The next work will be on geodetic datum.
>
> Martin
>
>
> [1] https://builds.apache.org/job/**sis-jdk7/site/apidocs/org/**
> apache/sis/referencing/**operation/matrix/package-**summary.html<https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/package-summary.html>
>
>
Report on matrix work (continuing)
Posted by Martin Desruisseaux <ma...@geomatys.fr>.
Hello all
I just completed the work on extended floating point precision in matrix
multiplication and inversion. This is new work compared to what we had
in Geotk. I added in the package javadoc [1] a small section, namely
"Extended floating point precision", which try to explain why we care
about the accuracy of those operations.
Now the last work needed is to port existing Geotk code for handling of
NaN values and some special kind of non-square matrix, and I should be
done with matrix. The next work will be on geodetic datum.
Martin
[1]
https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/package-summary.html
Re: Report on matrix work
Posted by Martin Desruisseaux <ma...@geomatys.fr>.
Thanks Adam :-)
Le 20/09/13 15:56, Adam Estrada a écrit :
> Thanks Martin! I am interested in what you find out in your research
> regarding floating point precisions. I (we) have been noticing some
> precision inconsistencies between 64bit Python and 64bit Java and I am
> curious to learn what others may encounter.
It depends on which operation the inconsistencies are found:
*
On the basic +, -, *, / operations, it depends on whether Phython
makes use of the Intel 80 bits registers or not. Java explicitly
forbid that (actually: mandates to truncate to 64 bits after every
operations, except for the exponent part if the "strictfp" keyword
is not present). So if the Phython language does not impose such
restriction (which I do not know), then Phython is likely to be more
accurate than Java on Intel processors, at the expense of slightly
different results on different architectures.
*
For trigonometric operations (sin, cos, tan, etc.), it depends on
whether Phython checks for large angles or not. The x87 Intel
processors are accurate for angles between -45° to +45°. For angles
greater than that, the processor uses an inaccurate approximation of
Pi. The error is small for small angles, but become greater as the
angle increase. Makers can build more accurate processors today, but
Intel can not do that on Pentium for compatibility reasons (I have
been told that AMD tried many years ago, but had to revert back to
the inaccurate algorithm because producing the correct answer was
breaking too many softwares). See the "Evaluation" section of Bug ID
4857011 [1] for more information (very interesting in my opinion).
Note that I do not know is using SSE unit instead of x87, as
proposed in [2], would affect accuracy.
About SIS matrix, the situation is as below:
*
Matrix multiplications and inversions are usually performed only
once when preparing a map projection. Then the result is used for
transforming a large bunch of coordinates. Consequently performance
of (matrix*matrix) and (matrix ^ -1) operations do not really matter
for Apache SIS. Performance of (matrix*vector) do matter a lot, but
this will be the job of an other class (MathTransform - to be
committed later).
*
In Apache SIS, the numerical values in the matrices will often be
specified by the CRS definition. The values may be conversion factor
from feet to metres (0.3048), from metre to kilometres (0.001), from
gradians to degrees (0.9), or map projection parameters defined by
the user (false easting/northing, etc.). It is quite common that a
chain of operations go from CRS A to CRS B to CRS C than back to CRS
A. After such circle, we often see 0.3047999999 instead of 0.3048,
etc. It is this kind of errors that I would like to reduce. While
apparently minor, we found the annoyance to have surprising large
implications.
In summary, for SIS matrices, accuracy is more important than
performance. Yesterday I was considering to use java.math.BigDecimal for
internal computation (not to be visible to users). I was concerned by
the high cost of such objects, but I though that the fact that
BigDecimal works in base 10 could actually be good for SIS matrices
because many factors are defined in base 10 (e.g. conversion from feet
to metres is defined as exactly 0.3048, which does not have an exact
representation as a Java double). However BigDecimal does not support
NaN and infinities, and we really needs the SIS matrices package to be
able to work with NaN/infinities.
Today I'm exploring the use of "double-double" arithmetic as an
alternative [3]. It would be more compact, more efficient, support
NaN/infinities and would hopefully have enough precision for our needs
regarding 0.3048 and similar factors.
Martin
[1] http://bugs.sun.com/bugdatabase/view_bug.do?bug_id=4857011
[2] http://bugs.sun.com/bugdatabase/view_bug.do?bug_id=7175279
[3]
http://en.wikipedia.org/wiki/Double-double_%28arithmetic%29#Double-double_arithmetic
Re: Report on matrix work
Posted by Adam Estrada <es...@gmail.com>.
This is really amazing work too, BTW! I am totally blown away :)
Adam
On Fri, Sep 20, 2013 at 9:56 AM, Adam Estrada <es...@gmail.com>wrote:
> Thanks Martin! I am interested in what you find out in your research
> regarding floating point precisions. I (we) have been noticing some
> precision inconsistencies between 64bit Python and 64bit Java and I am
> curious to learn what others may encounter.
>
> Adam
>
>
> On Thu, Sep 19, 2013 at 12:51 PM, Martin Desruisseaux <
> martin.desruisseaux@geomatys.fr> wrote:
>
>> Hello all
>>
>> The static utility methods in Matrices class [1] should now be completed.
>> There is not so many of them (I tried to trim down compared to Geotk), but
>> those which are there are important. In particular, 'createTransform(**AxisDirection[],
>> AxisDirection[])' will be involved in about every coordinate
>> transformations, so a bug in this method would have serious consequences
>> for SIS (thankfully, this method is relatively simple). I tried to improve
>> the javadoc with examples.
>>
>> If some peoples would like to look at the javadoc of the following
>> Matrices methods, I would like to fix any inconsistency or unclear
>> statements that may be identified:
>>
>> * createTransform
>> * createDimensionSelect
>> * createPassThrough
>>
>> The work that remain to be done is to port the special handling of NaN
>> values, inversion of non-square matrices, and maybe try to perform
>> multiplications and inversions using extended floating point precisions
>> (could be the subject an other email).
>>
>> Martin
>>
>>
>> [1] https://builds.apache.org/job/**sis-jdk7/site/apidocs/org/**
>> apache/sis/referencing/**operation/matrix/Matrices.html<https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/Matrices.html>
>>
>>
>
Re: Report on matrix work
Posted by Adam Estrada <es...@gmail.com>.
Thanks Martin! I am interested in what you find out in your research
regarding floating point precisions. I (we) have been noticing some
precision inconsistencies between 64bit Python and 64bit Java and I am
curious to learn what others may encounter.
Adam
On Thu, Sep 19, 2013 at 12:51 PM, Martin Desruisseaux <
martin.desruisseaux@geomatys.fr> wrote:
> Hello all
>
> The static utility methods in Matrices class [1] should now be completed.
> There is not so many of them (I tried to trim down compared to Geotk), but
> those which are there are important. In particular, 'createTransform(**AxisDirection[],
> AxisDirection[])' will be involved in about every coordinate
> transformations, so a bug in this method would have serious consequences
> for SIS (thankfully, this method is relatively simple). I tried to improve
> the javadoc with examples.
>
> If some peoples would like to look at the javadoc of the following
> Matrices methods, I would like to fix any inconsistency or unclear
> statements that may be identified:
>
> * createTransform
> * createDimensionSelect
> * createPassThrough
>
> The work that remain to be done is to port the special handling of NaN
> values, inversion of non-square matrices, and maybe try to perform
> multiplications and inversions using extended floating point precisions
> (could be the subject an other email).
>
> Martin
>
>
> [1] https://builds.apache.org/job/**sis-jdk7/site/apidocs/org/**
> apache/sis/referencing/**operation/matrix/Matrices.html<https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/Matrices.html>
>
>
Re: Report on matrix work
Posted by Martin Desruisseaux <ma...@geomatys.fr>.
Hello all
The static utility methods in Matrices class [1] should now be
completed. There is not so many of them (I tried to trim down compared
to Geotk), but those which are there are important. In particular,
'createTransform(AxisDirection[], AxisDirection[])' will be involved in
about every coordinate transformations, so a bug in this method would
have serious consequences for SIS (thankfully, this method is relatively
simple). I tried to improve the javadoc with examples.
If some peoples would like to look at the javadoc of the following
Matrices methods, I would like to fix any inconsistency or unclear
statements that may be identified:
* createTransform
* createDimensionSelect
* createPassThrough
The work that remain to be done is to port the special handling of NaN
values, inversion of non-square matrices, and maybe try to perform
multiplications and inversions using extended floating point precisions
(could be the subject an other email).
Martin
[1]
https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referencing/operation/matrix/Matrices.html
Re: Report on matrix work
Posted by "Mattmann, Chris A (398J)" <ch...@jpl.nasa.gov>.
This sounds great to me, excellent work, Martin.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Chris Mattmann, Ph.D.
Senior Computer Scientist
NASA Jet Propulsion Laboratory Pasadena, CA 91109 USA
Office: 171-266B, Mailstop: 171-246
Email: chris.a.mattmann@nasa.gov
WWW: http://sunset.usc.edu/~mattmann/
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Adjunct Assistant Professor, Computer Science Department
University of Southern California, Los Angeles, CA 90089 USA
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-----Original Message-----
From: Martin Desruisseaux <ma...@geomatys.fr>
Organization: Geomatys
Reply-To: "dev@sis.apache.org" <de...@sis.apache.org>
Date: Wednesday, September 11, 2013 4:36 PM
To: Apache SIS <de...@sis.apache.org>
Subject: Report on matrix work
>Hello all
>
>The javadoc package documentation below tries to explain the rational
>behind the SIS matrix package (still work in progress):
>
>https://builds.apache.org/job/sis-jdk7/site/apidocs/org/apache/sis/referen
>cing/operation/matrix/package-summary.html
>
>The idea is strongly inspired by "vecmath", which was part of Java3D.
>"Vecmath" is actually hard to replace by other matrix package (JBLAS,
>JAMA or other) because the design goals are not the same. There is
>thousands of possible operations for matrix, and the operations of
>interest are not the same for resolving linear equations than for
>performing coordinate transformations. Because "vecmath" was designed
>for coordinate conversions in 3D spaces, it was a pretty good fit to
>geospatial needs.
>
>There is a few methods in org.apache.sis.referencing.operation.matrix
>that are more closely related to coordinate transformations rather than
>general matrix operations: isAffine(), toAffineTransform(),
>normalizeColumns() - more will come soon. I have not yet connected the
>SIS package to JAMA. But I came to think that the connection (if any)
>should not be visible in the public API. For example when inverting a
>non-square matrix, JAMA uses a least squares approach. This is probably
>a good choice for resolving linear equations, but sis-referencing needs
>a more conservative approach. For this reason, I don't think that we
>should expose the JAMA Matrix.inverse() method directly.
>
>I will post more when the package will be closer to completion, unless
>someone would like clarifications or changes?
>
> Martin
>