Octave left division not robust enough compared to Matlab

[David Bateman]

Alexander Mamonov wrote:
> Consider the following example:
> for some vector v compute
> x = vander(v) \ v
> The answer is obvious, x is all zeros, except for x(end-1), which is one.
> Since the Vandermonde matrix is known to be exponentially ill-conditioned,
> computing x is a challenging task for a linear solver.
> If we take v in the form
> v = linspace(1, 2, N)';
> then for N=14 in Octave we already have
> 
>> x = vander(v) \ v
> warning: matrix singular to machine precision, rcond = 4.50748e-18
> warning: attempting to find minimum norm solution
> warning: dgelsd: rank deficient 14x14 matrix, rank = 13
> x =
> 
>   -4.5700e-007
>   8.9146e-006
>   -7.9956e-005
>   4.3658e-004
>   -1.6191e-003
>   4.3063e-003
>   -8.4507e-003
>   1.2388e-002
>   -1.3565e-002
>   1.0959e-002
>   -6.3481e-003
>   2.4969e-003
>   9.9940e-001
>   6.5752e-005
> 
> which is nowhere close to the correct solution.
> On the other hand Matlab produces exact answer for N up to 64
> (RCOND = 1.519193e-036).
> For the above case N=14 Matlab produces the following output
> 
>>> x = vander(v) \ v
> Warning: Matrix is close to singular or badly scaled.
>          Results may be inaccurate. RCOND = 4.490508e-018.
> 
> x =
> 
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      0
>      1
>      0
> 
> The comparison was made with
> Octave 3.0.2 Windows mingw build, and Matlab Windows R2006a.
> Similar behavior was observed for Octave 3.0.1 and 3.0.2 on Linux (Debian).
> 
> While not a language incompatibility, this issue is a "numerical
> incompatibility"
> of Octave with Matlab, i.e. the problems that are easily solved on
> Matlab produce
> unexpected behavior in Octave.
> 
> The above example may seem artificial, but I stumbled upon similar behavior
> in my actual research, when I tried to solve a particular
> ill-conditioned inverse problem.
> 
> Since I'm not familiar with Octave internal workings, it's difficult
> for me to judge
> what could be the source of the problem. My guess would be that single
> presicion is used somewhere in the left divide computation, which increases
> the sensitivity to the conditioning of the problem significantly.
> Another possible
> reason is that Octave switches to minimum-norm solver too early.

Read the thread,

http://www.nabble.com/Behavior-of-mldivide-in-Octave-relative-to-Matlab-to16561235.html#a16561235

in short the matlab behavior is wrong, as it uses an LUP factorization 
for a square matrix even if it is flagged as singular and ends up with a 
factorization that is not invariant with column permutations.

The thread above discusses how you can get the same behavior if you 
really want it, even though Octave's solution is numerically superior to 
the Matlab one.

D.
-- 
David Bateman                                dbateman at dbateman.org
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)