Complex left-divide and SVD give ireproducible results

[David Bateman]

Dirk Laurie wrote:
> --------
> Bug report for Octave 2.1.73 configured for i486-pc-linux-gnu    
>
> Description:
> -----------
>
> The following operations give irreproducible results:
>     z=A\b  when A is a nonsquare complex matrix
>     [U,S,V]=svd(A) when A is a complex matrix
>
> By 'irreproducible' is meant that the results are not identical
> to the last bit when the operation is run more than once on
> identical input, although the results in both cases are numerically
> correct in the sense that b-A*z and A-U*S*V respectively are at
> roundoff error level.  The two bugs are reported as one because 
> it is thought probable that they share a common cause.
>
> A possible cause might be that an uninitialized value is used in 
> a comparison of which the result would normally not critically 
> affect the output, e.g. whether columns in a full-rank orthogonal 
> factorization are to be interchanged. 
>   

Note that 2.1.73 is an older release and you should probably be using
3.0.1, but that doesn't change this discussion. In fact due to the
change in the use of the lapack xGELSS function to use xGELSY instead in
recent versions of Octave, with your test case I can no longer exposes
the issue. There has been many e-mails about issues similar to this to
the Octave lists., and its related to many many different operators and
functions in Octave of floating point values.

Basically, I think you just have to live with it for the latest
compilers and if you use an optimized blas. My take on the reason for
such issues is that the x86 processors have 80 bit internal registers
and optimized blas implementations and the 4.0.x release of the gnu tool
chain and later make use of these additional bits of precision. That is
between operations on the x86 registers, the floating-point values don't
have to be stored back to the 64-bit value that contains it saving the
cost of a memory copy from and back to the registers. This can give
significant improvement in the performance and a marginal (and
unpredictable) improvement in the accuracy.

However, during a context switch on the x86 processors, all of the
registers are swapped back to memory and those 64-bit values and because
you can't control when the context switches might happen when the matrix
sizes start to get large you can get unpredictable differences in the
numerical solutions.

So what can you do if you absolutely have to always have the same
numerical values

* Don't use optimized version of the blas/lapack
* Use a gnu toolchain earlier than 4.0.x or use the -ffloat-store option
to g++, gcc and gfortran.

Though the cost in terms of speed will be horrendous. Also, if you
algorithm is so sensitivity to those last couple of bits of accuracy you
probably need to reconsider your algorithm. In any case I don't consider
this to be a bug in Octave as such.

Cheers
David


-- 
David Bateman                                David.Bateman at motorola.com
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary