From bug-octave-request at che dot utexas dot edu  Tue Oct 26 13:40:24 1993
Subject: Solving Ax=b
From: Tony Mullins <jamull at che dot utexas dot edu>
To: bug-octave
Date: Tue, 26 Oct 93 13:40:20 EDT


I believe that there is an inconsistency in the reporting of solutions
to Ax=b when A lacks full rank.  For non-square A, Octave returns a
minimum norm, least squares solution.  Why not continue that behavior
for square A?  Consider the following
example:

octave:1> A = [1,2;2,4]
A =

  1  2
  2  4

octave:2> b = [3;5]
b =

  3
  5

octave:3> x =A\b 
error: matrix singular to machine precision, rcond = 0
error: evaluating binary operator `\' near line 3, column 5
error: evaluating assignment expression near line 3, column 3


Here is what I believe the correct behavior should be:

octave:4> [V,S,U] = svd(A)
V =

  -0.44721  -0.89443
  -0.89443   0.44721

S =

  5.00000  0.00000
  0.00000  0.00000

U =

  -0.44721  -0.89443
  -0.89443   0.44721

octave:6> r = rank(A)
r = 1

octave:10> [m,n] = size(A)
m = 2

n = 2

octave:11> S_plus = zeros(m,n)
S_plus =

  0  0
  0  0

octave:12> S_plus(1:r,1:r) = inv(S(1:r,1:r))
S_plus =

  0.20000  0.00000
  0.00000  0.00000

octave:13> A_plus = U*S_plus*V'
A_plus =

  0.040000  0.080000
  0.080000  0.160000

octave:14> x = A_plus*b
x =

  0.52000
  1.04000

This x is the minimum norm solution to Ax=p where p is the projection
of b onto the range of A.
I don't know what the behavior of Matlab is, but this behavior seems
consistent with the effect of \ for non-square systems.  I realize
that this complicates the analysis of singularity for square systems.
But if singularity is detected for square systems, a solution still
could be returned with the message that the solution is a minimum norm
least squares solution to Ax=b.

Tony Mullins
jamull at che dot utexas dot edu
Department of Chemical Engineering
The University of Texas at Austin