Re: Behavior of mldivide in Octave relative to Matlab

[Jaroslav Hajek]

On Tue, Apr 8, 2008 at 12:57 PM, David Bateman
<address@hidden> wrote:
> I was recently contacted by Marco Caliari in particular concerning the
>  behavior of the mldivide operator in Octave and Matlab for
>  underdetermined systems. He pointed out the thread
>
>  http://icl.cs.utk.edu/lapack-forum/viewtopic.php?p=497&;
>
>  where the behavior of the Matlab mldivide operator is analyzed. It seems
>  that Matlab is not using the LAPACK xGEDLD function like Octave does,
>  but rather a formulation based on the QR factorization. The formulation
>  can be written as
>
>  m=10; n=20; A = randn(m,n); b=rand(m,1); [Q,R,E]=qr(A);
>  [A\b E(:,1:m)*(R(:,1:m)\(Q'*b))]
>
>  Result returned by this is not the minimum 2-norm solution, like in
>  Octave but is equally valid. It also appears that the matlab method is
>  faster and uses less memory. Marco points out that values like m=50 and
>  n=20000 in the above will result in an exhaustion of the memory of the
>  machine, whereas Matlab is still capable of returning a result.
>

While certainly faster, this approach seems to be inferior in other ways.

A numerical question arises: how big can the null space component become,
relative to the minimum-norm solution? Can it be nicely bounded, or can it be
arbitrarily big? Consider this example:

m = 10; n = 10000; A = ones(m,n)+1e-6*randn(m,n); b =
ones(m,1)+1e-6*randn(m,1); norm(A\b)

while Octave's minimum-norm values are around 3e-2, Matlab's results
are 50-times larger.

For another issue, try this code:

 m = 5; n = 100; j = floor(m*rand(1,n))+1; b = ones(m,1); A =
zeros(m,n); A(sub2ind(size(A),j,1:n)) = 1; x = A\b; [dummy,p] =
sort(rand(1,n)); y = A(:,p)\b; norm(x(p)-y)

It shows that unlike in Octave, mldivide in Matlab is not invariant
w.r.t. column permutation
if there are multiple columns of the same norm - i.e. permuting
columns of the matrix gets you different result than ermuting the
solution vector. I bet that many Matlab users would be surprised by
this behaviour.

In my opinion, since \ is often part of a more complex expression,
(where there is no room to react to warnings or flags) it should
prefer intelligence (robustness) to speed, and that speaks for
Octave's implementation.

>  Another issue is what is the result that is returned by Matlab and
>  Octave for singular matrices. Consider a case like
>
>  A = [eps, 0, -eps; 0, eps, -eps; -eps, 0, 10]; x = A \ ones(3,1)
>
>  Octave falls back to using xGELSD to solve the matrix in this case
>  whereas as Matlab returns the result based on LU factorization with a
>  warning of its inaccuracy.
>

This is again the intelligence vs. speed.


>  Should we duplicate Matlab's behavior for the two cases above?
>

I'd suggest making it an option, and keep Octave's behavior the default.

cheers,


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