CVS difference for ai05s/ai05-0047-1.txt

--- ai05s/ai05-0047-1.txt	2007/07/26 02:58:06	1.2
+++ ai05s/ai05-0047-1.txt	2007/08/03 04:14:09	1.3
@@ -1,4 +1,4 @@
-!standard G.3.1                              07-06-21    AI05-0047-1/02
+!standard G.3.1                              07-08-02    AI05-0047-1/03
 !standard G.3.2
 !class binding interpretation 07-04-04
 !status work item 07-04-04
@@ -73,6 +73,25 @@
     X * Conjugate (Y) which mimics the mathematical notation where conjugation is always
     made explicit.  At any rate, a clarification would be useful.
 
+7 - The eigenvectors returned by Eigensystem in Generic_Real_Arrays
+    and Generic_Complex_Arrays are not well specified, and this has
+    unfortunate effects. One problem is that the result of these subprograms
+    is hard to test. A more serious problem is that different implementations
+    may return different eigenvectors, and this may cause portability
+    problems. Also, even in a single implementation, changes to the internal
+    algorithms, or to the code generated for the annex G generics, may cause
+    the result of Eigensystem to change. This non-determinism is unpleasant.
+
+    Consider first the real case. Given a set of orthonormal eigenvectors,
+    changing the signs of any subset of these vectors also result is a set of
+    orthonormal eigenvectors (this is easily seen geometrically: it amount to
+    changing the direction of the vectors in an orthonormal basis). It would
+    be nice to be more specific.
+
+    Consider now the complex case. Things are more nasty here because there
+    is much more freedom: multiplying each vector by a complex number of
+    modulus 1 leaves the set of eigenvectors orthonormal. Again, it would be
+    nice to lift the ambiguity.
 
 !recommendation
 
@@ -87,5 +106,52 @@
 
 
 !appendix
+
+[Summary of private mail of July 2007; used by demand, er, permission.]
+
+Pascal Leroy:
+
+It appears that the eigenvectors returned by Eigensystem in
+Generic_Real_Arrays and Generic_Complex_Arrays are not well specified, and
+this has unfortunate effects.  One problem is that the result of these
+subprograms is hard to test.  A more serious problem is that different
+implementations may return different eigenvectors, and this may cause
+portability problems.  Also, even in a single implementation, changes to
+the internal algorithms, or to the code generated for the annex G
+generics, may cause the result of Eigensystem to change.  This
+non-determinism is unpleasant.
+
+Consider first the real case.  If V1, V2, ..., Vn are eigenvectors that
+are mutually orthonormal, then changing the signs of any subset of these
+vectors also result is a set of eigenvectors that are mutually orthonormal
+(this is easily seen geometrically: it amount to changing the direction of
+the vectors in an orthonormal basis).  It would have been nice to be more
+specific, for instance by requiring the first nonzero component of each
+vector to be positive.
+
+Consider now the complex case.  Things are more nasty here because you
+have much more freedom: you can multiply each vector by a complex number
+of modulus 1, and still have an orthonormal set.  One way to lift the
+ambiguity would be to require that the first nonzero component of each
+vector be a positive real.  But this would require a lot of complex
+divisions, which would degrade the accuracy somewhat.  Not sure what to do
+here.
+
+John Barnes:
+
+Another source of difficulty is if two or more eigenvalues are the same.
+This results in a whole subspace of  vectors and any orthogonal set in that
+subspace will do. I don't see how one could overcome that easily either.
+
+Pascal Leroy:
+
+True, but it's hard to be excited about this.  Even if the matrix has two
+mathematically identical eigenvalues, it's hard to believe that the
+numerical algorithm used will find the floating-point values of the
+eigenvalues be exactly identical, so the eigenvectors won't form a whole
+subspace after all.  If, conversely, it turns out that the matrix didn't
+have mathematically identical eigenvalues, but the numerical algorithm
+produces two identical floating-point values, then the whole thing is
+probably ill-conditioned, so you get what you get.
 
 ****************************************************************