A simpler fix: return rank=NA when Lapack is used. Keith
> On Jan 23, 2018, at 5:36 AM, Serguei Sokol <so...@insa-toulouse.fr> wrote: > > Le 23/01/2018 à 08:47, Martin Maechler a écrit : >>>>>>> Serguei Sokol <so...@insa-toulouse.fr> >>>>>>> on Mon, 22 Jan 2018 17:57:47 +0100 writes: >> > Le 22/01/2018 à 17:40, Keith O'Hara a écrit : >> >> This behavior is noted in the qr documentation, no? >> >> >> >> rank - the rank of x as computed by the decomposition(*): always full >> rank in the LAPACK case. >> > For a me a "full rank matrix" is a matrix the rank of which is indeed >> min(nrow(A), ncol(A)) >> > but here the meaning of "always is full rank" is somewhat confusing. >> Does it mean >> > that only full rank matrices must be submitted to qr() when >> LAPACK=TRUE? >> > May be there is a jargon where "full rank" is a synonym of >> min(nrow(A), ncol(A)) for any matrix >> > but the fix to stick with commonly admitted rank definition (i.e. the >> number of linearly independent >> > columns in A) is so easy. Why to discard lapack case from it (even >> properly documented)? >> >> Because 99.5% of caller to qr() never look at '$rank', >> so why should we compute it every time qr() is called? > 1. Because R already does it for linpack so it would be consistent to do so > for lapack too. > 2. Because R pretends that it is a part of a returned qr class. > 3. Because its calculation is a negligible fraction of QR itself. > >> >> ==> Matrix :: rankMatrix() does use "qr" as one of its several methods. >> >> -------------- >> >> As wiser people than me have said (I'm paraphrasing, don't find a nice >> citation): >> >> While the rank of a matrix is a very well defined concept in >> mathematics (theory), its practical computation on a finite >> precision computer is much more challenging. > True. It is indeed depending of round-off errors during QR calculations and > tolerance > setting but putting it just as min(nrow(A), ncol(A)) and still calling it > rank of "full rank" > is by far the most misleading choice to my mind. > > Once again, if we are already calculating it for linpack let do it in most > consistent > way for lapack too. I can propose a patch if you will. > > Serguei. > >> >> The ?rankMatrix help page (package Matrix, part of your R) >> https://stat.ethz.ch/R-manual/R-devel/library/Matrix/html/rankMatrix.html >> starts with the following 'Description' >> >> __ Compute ‘the’ matrix rank, a well-defined functional in theory(*), >> somewhat ambigous in practice. We provide several methods, the default >> corresponding to Matlab's definition. >> >> __ (*) The rank of a n x m matrix A, rk(A) is the maximal number of linearly >> independent columns (or rows); hence rk(A) <= min(n,m). >> >> >> >>> On Jan 22, 2018, at 11:21 AM, Serguei Sokol <so...@insa-toulouse.fr> >> wrote: >> >>> >> >>> Hi, >> >>> >> >>> I have noticed different rank values calculated by qr() depending on >> >>> LAPACK parameter. When it is FALSE (default) a true rank is >> estimated and returned. >> >>> Unfortunately, when LAPACK is set to TRUE, the min(nrow(A), ncol(A)) >> is returned >> >>> which is only occasionally a true rank. >> >>> >> >>> Would not it be more consistent to replace the rank in the latter >> case by something >> >>> based on the following pseudo code ? >> >>> >> >>> d=abs(diag(qr)) >> >>> rank=sum(d >= d[1]*tol) >> >>> >> >>> Here, we rely on the fact column pivoting is activated in the called >> lapack routine (dgeqp3) >> >>> and diagonal term in qr matrix are put in decreasing order >> (according to their absolute values). >> >>> >> >>> Serguei. >> >>> >> >>> How to reproduce: >> >>> >> >>> a=diag(2) >> >>> a[2,2]=0 >> >>> qaf=qr(a, LAPACK=FALSE) >> >>> qaf$rank # shows 1. OK it's the true rank value >> >>> qat=qr(a, LAPACK=TRUE) >> >>> qat$rank #shows 2. Bad, it's not the expected value. >> >>> >> >> > -- >> > Serguei Sokol >> > Ingenieur de recherche INRA >> >> > Cellule mathématique >> > LISBP, INSA/INRA UMR 792, INSA/CNRS UMR 5504 >> > 135 Avenue de Rangueil >> > 31077 Toulouse Cedex 04 >> >> > tel: +33 5 6155 9849 >> > email: so...@insa-toulouse.fr >> > http://www.lisbp.fr [[alternative HTML version deleted]] ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
