Julien Langou wrote: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.)
But the default value is 2^-52=2.22e-16 on PC right?
* 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.
Julien Langou wrote: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.
Julien Langou wrote:* 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.)
Users browsing this forum: No registered users and 1 guest