## About DGELSD subroutine and it's rank estimation

Open discussion regarding features, bugs, issues, vendors, etc.

### About DGELSD subroutine and it's rank estimation

Hello,
As the help text of xGELSD subroutine says, the rank of the A matrix is estimated via svd decomnposition as "the number of singular values which are greater than RCOND*S(1)".
In GNU Octave software, xGELSD is used to solve a system of equations in the left division operator '\'. In GNU Octave a function called rank() exists too in order to determine the numerical rank of a matrix via svd decomposition. The help text of rank() function is:
Code: Select all
` -- Function File:  rank (A, TOL)     Compute the rank of A, using the singular value decomposition.     The rank is taken to be the number of singular values of A that     are greater than the specified tolerance TOL.  If the second     argument is omitted, it is taken to be          tol = max (size (A)) * sigma(1) * eps;     where `eps' is machine precision and `sigma(1)' is the largest     singular value of A.`

As we can see, rank() uses as criteria for consider a singular value as zero the expression max (size (A)) * sigma(1) * eps (in Octave, eps is the precision, equivalent to DLAMCH('P')).

The '\' operator in Octave uses RCOND=-1 in xGELSD (http://hg.savannah.gnu.org/hgweb/octave ... x.cc#l2430), so xGELSD uses internally the machine precision (help of http://www.netlib.org/lapack/double/dgelsd.f).

I have problems about rank computation with DGELSD (via \ operator) vs rank() in GNU Octave (my installation uses lapack 3.2.2 from Debian repositories and GNU Octave 3.2.4).

For example, we can generate a random matrix with rank deficiency and try to solve a system with '\' operator and try to calculate the rank of the matrix with rank() function as:
Code: Select all
`a=rand(1000,200);a(:,196:200)=a(:,1:5);b=a\a;rank(a)warning: dgelsd: rank deficient 1000x200 matrix, rank = 196ans =  195`

As we can see, by construction the matrix has a rank of 195 because the 5 last columns are the same than the first 5 columns. The operator rank() computes exactly a rank of 195, but the function DGELSD (via \) computes a rank of 196.

If we change the function rand() by randn() in the generation of the matrix, DGELSD says that the matrix has complete rank (no warning about rank deficiency appears), but rank() function again computes the correct result
Code: Select all
`a=randn(1000,200);a(:,196:200)=a(:,1:5);b=a\a;rank(a)ans =  195`

My knowledge about numerical algebra is limited but apparently the criteria max(size(A))*sigma(1)*RCOND (RCOND==eps in Octave) is better than sigma(1)*RCOND in rank revealing. Are I'm wrong?

Thanks
jgpallero

Posts: 24
Joined: Thu Jul 29, 2010 2:29 pm