Re: Behavior of mldivide in Octave relative to Matlab

[Jaroslav Hajek]

On Wed, Apr 9, 2008 at 9:18 AM, Marco Caliari <address@hidden> wrote:
>
> >
> > >  Only a final note: the case m=50 and n=20000 works in Matlab's A\b
> whereas
> > > it fails in Matlab's [Q,R,E]=qr(A) (same problem of exhaustion of the
> memory
> > > as in Octave).
> > >
> >
> > I think you can do `[Q,R,E] = qr (A,0)' instead, to get E as
> > permutation vector rather than matrix.
> >
>
>  To summarize: in order to imitate the Matlab's behaviour of A\b, you can do
>
>  [Q,R,E] = qr(A,0);
>  x = (R(:,1:m)\(Q'*b));
>  x(n) = 0;
>  x(E) = x;
>
>  Thanks for the hint, Jaroslav. Corcerning the underdetermined case, I'm
> completely satisfied now (yes, David, I have a code depending on sparse
> solutions of underdetermined systems), because I can reproduce Matlab with
> the four clear lines above.

If you think it is useful, I suggest you contribute it to Octave or
Octave-Forge.


>  Concerning near-singular matrices, for the script
>
>  for i = 12:16
>   A = [10^(-i), 0, -10^(-i); 0, 10^(-i), -10^(-i); -10^(-i), 0, 10];
>   b = ones(3,1);
>   x = A \ b;
>   norm(A*x-b)
>  end
>
>  I prefer Matlab's result. So, I would appreciate a different function or an
> optional argument for mldivide. Not sure what should be the default
> behaviour.
>

OK, that's a question of preference. Consider again a counter-example:

# solve, nullify a tiny element of A, solve again, get relative difference
eps = 1e-17; A = [eps,0,-eps;0,eps,-eps;-eps,0,10]; b = ones(3,1); x1
= A\b; A(1,1) = 0; x2 = A\b; norm(x1-x2)/max(norm(x1),norm(x2))

Octave's result is on the level of machine epsilon, while in Matlab x1
and x2 are off by 100%.
The question, of course, is, whether the change A(1,1) = 0 is actually tiny.
>From the point of view of the whole matrix, yes (matrix norm is ~ 10).
>From the point of view of its row or column, no. In practice, the
answer will depend on how the matrix was created.

Matlab seems to always use LUP factorization in the square case, so
that's where its behaviour comes from. Again, to make it more
counter-intuitive, try instead:

# add one null equation (zero row and zero rhs) and watch what happens.
eps = 1e-17; A = [eps,0,-eps;0,eps,-eps;-eps,0,10]; b = ones(3,1); x1
= A\b; x2 = [A;zeros(1,3)]\[b;0]; norm(x1-x2)/max(norm(x1),norm(x2))

again, I can easily imagine Matlab users being surprised here.

Both approaches are valid and useful (though in different situations).
I think that Octave's approach is somewhat more consistent, in the
sense that it does not treat the square case much specially.
Again, a special function can be provided to get Matlab-like
behaviour. A fairly simple and robust approach would be just wrap the
LAPACK's xGESVX solvers (which balance the matrix, perform iterative
refinement and compute error estimates). That would involve a DLD
function (in C++) - if you don't feel up to it, you can leave this one
to me, I'll do it in some near future.

cheers,

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz