pseudo inverse for ill-conditioned matrix

[fkengun]

Now I am using CLAPACK, I want to calculate the pseudo inverse for ill conditioned matrix by calling dgelss and setting B to a identity matrix. The matrix is overdetermined(100*64) and not full rank. The condition number is 2.2e+16. The result I got is different from that of MATLAB. There is no huge number in the result so it seems right, but actually it's wrong(if MATLAB is right). But for a small matrix, which is 2*2 with a condition number 2.2e+8, the results are both right.

I'm wonderring is there any way that I can get a right result? Or should I call dgesvd and do the rest calculate myself? Or should I treat the matrix first to eliminate the ill condition before doing the dgelss?

For deeper understanding, what's the different between the algorithms of pinv in MATLAB and dgelss in CLAPACK? Is there any additional ill condition elimination treatment in MATLAB?


I appreciate any comment.

Best regards,

[Julien Langou]

The pseudoinverse is not a continuous function of the matrix entries, the discontinuity occurs in particular at the interface of singular and nonsingular matrices. Being near a discontinuity is really a bad place to be when doing numerical computation.

If the entries in the pseudoinverse do not blow up using Matlab's pinv() and the ones in LAPACK DGELSS do, it is probable that Matlab is using some kind of filtering to remove the small singular values in your initial matrix.
Setting RCOND to 1.0D-10 in DGELSS may help in your case. (You can also try RCOND < 0.)

Julien.

[fkengun]

The pseudoinverse is not a continuous function of the matrix entries, the discontinuity occurs in particular at the interface of singular and nonsingular matrices. Being near a discontinuity is really a bad place to be when doing numerical computation.

If the entries in the pseudoinverse do not blow up using Matlab's pinv() and the ones in LAPACK DGELSS do, it is probable that Matlab is using some kind of filtering to remove the small singular values in your initial matrix.
Setting RCOND to 1.0D-10 in DGELSS may help in your case. (You can also try RCOND < 0.)

Julien.

Thanks for your reply first.
That is what I think, MATLAB uses some kind of processing to prevent itself from blowing up. And I want to know how it is done.
I will try to configure the RCOND. But the default value is 2^-52=2.22e-16 on PC right? I thought it is the precision value which is the higher the better, am I right? Anyway, I will try it first.

Thanks again.
Best Regards

[Julien Langou]

But the default value is 2^-52=2.22e-16 on PC right?

Well ... what do you mean by default ? ...

Code: Select all

*  RCOND   (input) DOUBLE PRECISION
*          RCOND is used to determine the effective rank of A.
*          Singular values S(i) <= RCOND*S(1) are treated as zero.
*          If RCOND < 0, machine precision is used instead.


The subroutine checks RCOND. If RCOND is set greater than 0, (RCOND > 0,) then the subroutine uses RCOND as threshold to cut the singular values ( S( i ) ) smaller than THR = RCOND * S( 1 ) where S( 1 ) is the largest singular value. If RCOND is set to zero, (RCOND = 0,) then the subroutines does not filter any singular values except the exactly zero ones if any. (Same logic as RCOND > 0.) If RCOND is set less than 0, (RCOND < 0,) then the subtoutine uses EPS (the machine precision) to cut the singular values ( S( i ) ) smaller than THR = EPS * S( 1 ). So RCOND < 0 is the same as RCOND = EPS.

So there is in some sense "no default". You do need to set RCOND to a value. The default one might be a negative number, but you need to take action.
If you do nothing and have RCOND uninitialized in entry, this is a pretty bad idea.
You may filter, not filter, etc .

Two comments:

1) You probably want to use DGELSD, it is the same functionality as DGELSS, only faster. DGELSY might also be interesting to give a shot. It is not the same algorithm but should give a good answer (i.e., filtered as you ask for) in most of the cases. DGELSY would be really faster than bother DGELSS and DGELSD.

2) Calling this subroutine with ( RCOND > 1 ) makes no sense whatsoever - no matter what the matrix is. (The subroutine will behave at it specifies. Pseudoinverse would be identity. Rank would be 0. Fun.)

Cheers,
Julien.

[fkengun]

Thanks for your huge reply.

The default value in my words is exactly the EPS. Usually I set RCOND to -1, which means use EPS to cut singular value. But when I called DGELSS to get pseudo inverse with RCOND=-1, the result blew up. If I set RCOND to 1e-15, it worked. My question is why it works better if I set the precision to a low value.

About the DGELSD and DGELSY, I will check them out today. Thanks for your comments

Regards

[sven]

Your matrix probably has a singular value close to EPS, in which case you are likely to get a big difference depending upon whether or not that singular value is treated as zero.

I believe that the MATLAB pinv function uses a tolerance of EPS*max(m,n)*max(singular value), so if you wish to compare LAPACK and MATLAB you should take RCOND = EPS*max(m,n) in DGELSS. pinv and DGELSS both use the SVD.

It might be worth trying your original suggestion of computing the SVD and then forming A^+ = V S^+ U^H, where X^+ denotes the pseudo-inverse. S^+ should be the diagonal matrix whose nonzero entries are the reciprocals of the diagonal entries of those element of S that you regard as nonzero. That is, set any diagonal elements of S that you regard as negligible to zero, and reciprocate the rest. Of course, the number of nonzero entries is the numerical rank of A. By looking at the diagonal entries of S you will be able to see the distribution of the singular values.

Best wishes,

Sven Hammarling.

[fkengun]

1) You probably want to use DGELSD, it is the same functionality as DGELSS, only faster. DGELSY might also be interesting to give a shot. It is not the same algorithm but should give a good answer (i.e., filtered as you ask for) in most of the cases. DGELSY would be really faster than bother DGELSS and DGELSD.

I tried the DGELSD and DGELSY. The DGELSD does run faster than DGELSS(65us against 75us), and the result is right. But the DGELSY didn't return right answer to me, no matter how I configured the parameter. I set the JPVT to all zero to disable the column permutation and all the rest parameter are almost the same as DGELSS. I don't think there is anything wrong. So I'm confused.

[Julien Langou]

* Your matrices are fairly small so the difference DGESLD is DGELSS is not that important. For larger matrices, the improvement is much more significant.

* Regarding DGELSY, this is not exactly the same algorithm, the answer should have satisfied you, if not, forget it. It was good to try though. (I think the RCOND needs to be set a little higher for DGELSY. But there is not a one-to-one match between DGELSS/DGELSD and DGESLY.)

Julien

[fkengun]

* Your matrices are fairly small so the difference DGESLD is DGELSS is not that important. For larger matrices, the improvement is much more significant.

* Regarding DGELSY, this is not exactly the same algorithm, the answer should have satisfied you, if not, forget it. It was good to try though. (I think the RCOND needs to be set a little higher for DGELSY. But there is not a one-to-one match between DGELSS/DGELSD and DGESLY.)

Julien

My matrix is 100*64. And that will be the standard size in my program. DGELSD, as a faster function, will definitely be my choice.

I still have some confusion. The parameter RCOND in all the functions should be the precision during the processing. But why sometime we should set it a lower value than the machine precision 2.22e-16? I thought the precision should be the higher the better. Could you please explain a little bit detail of that to me?

Regards

[Julien Langou]

If you need to set RCOND lower than EPS to get "good results", (whatever this means,) this is probably because your matrix has some particular structure. (I am guessing.) For example, some columns are scaled way smaller than others.This is not the general case. In the general case, you would expect the condition number to be at most the reciprocal of EPS. (This is because round off errors naturally brings everything to that level.) RCOND is there and you can tune it to adapt to your case. If setting it smaller than EPS is "better" for you, please do so. Julien.