Re: x = Z\z

[David Bateman]

oort wrote:
> Hello.
>
> If x = Z\z is the solution of Zx=z and only square systems have solution
> then why the operation of non-square matrices Z gives "numerical values"?
> Don't you think that it should give some error message? 
>
> For instance:
>
> octave:1> A = [3, 2, 6; 2, -2, 1; -1, 0.5, 3]
> A =
>
>    3.00000   2.00000   6.00000
>    2.00000  -2.00000   1.00000
>   -1.00000   0.50000   3.00000
>
> octave:2> a = [1; 2; 3]
> a =
>
>    1
>    2
>    3
>
> octave:3> A\a
> ans =
>
>   -0.74684
>   -1.26582
>    0.96203
>
> OK... "A" is a 3x3 matrice and "a" is a 3x1 matrice 
>
> But:
>
> octave:4> B = [3, 2; 2, -2; -1, 0.5]
> B =
>
>    3.00000   2.00000
>    2.00000  -2.00000
>   -1.00000   0.50000
>
> octave:5> b = [1; 2; 3]
> b =
>
>    1
>    2
>    3
>
> octave:6> B\b
> ans =
>
>    0.29801
>   -0.11479
>
> Now we have a system of 3 equations and 2 variables. It's a overdefined
> system. Curiously B*x does not give equal to "b"...
>
>
> And if we use a underdefined system we aldo reach a "numeric" result.
>
> octave:7> C = [3, 2, 6; 2, -2, 1]
> C =
>
>    3   2   6
>    2  -2   1
>
> octave:8> c = [1; 2]
> c =
>
>    1
>    2
>
> octave:9> C\c
> ans =
>
>    0.42175
>   -0.51459
>    0.12732
>
> How it is possible to have a "result" from a underdifined system?
> Curiously  if we put a third row with zeros in C and if we calculate C*x we
> obtain "c". However in octave x = [0.42175; -0.51459; 0.12732] and in MATLAB
> x = [0; -0.7857; 0.4286] 
>   

Although overdefined systems have no solution, a minimum norm supplies
the solution that is closest. That is Octave return the value of x that
minimizes norm(A*x - b)

Similarly for underdefined systems there are infinitely many solutions
and Octave return the unique solution that minimizes the norm of X.

Matlab does things differently for these cases, and in our minds
incorrectly. Copying from the draft version of the FAQ


Matlab's solvers as used by the operators mldivide (\) and mrdivide (/),
use a different approach than Octave's in the case of singular, under-,
or over-determined matrices. In the case of a singular matrix, Matlab
returns the result given by the LU decomposition, even though the underlying
solver has flagged the result as erroneous. Octave has made the choice
of falling back to a minimum norm solution of matrices that have been
flagged as singular which arguably is a better result for these cases.

In the case of under- or over-determined matrices, Octave continues to
use a minimum norm solution, whereas Matlab uses an approach that is
equivalent to

@example
@group
function x = mldivide (A, b)
  [Q, R, E] = qr(A);
  x = [A \ b, E(:, 1:m) * (R(:, 1:m) \ (Q' * b))]
end
@end group
@end example

@noindent
While this approach is certainly faster and uses less memory than
Octave's minimum norm approach, 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:

@example
@group
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)
@end group
@end example

@noindent
while Octave's minimum-norm values are around 3e-2, Matlab's results
are 50-times larger. For another issue, try this code:

@example
@group
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)
@end group
@end example

@noindent
It shows that unlike in Octave, mldivide in Matlab is not invariant
with respect to column permutations. If there are multiple columns of
the same norm, permuting columns of the matrix gets you different
result than permuting the solution vector. This will surprise many
users.

Since the mldivide (\) and mrdivide (/) operators are 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 so the Octave
developers
are firmly of the opinion that Octave's approach for singular, under- and
over-determined matrices is a better choice that Matlab's