[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Behavior of mldivide in Octave relative to Matlab
From: |
Jaroslav Hajek |
Subject: |
Re: Behavior of mldivide in Octave relative to Matlab |
Date: |
Tue, 8 Apr 2008 14:50:17 +0200 |
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