From: Martin Maechler <maechler_at_stat.math.ethz.ch>

Date: Sat 07 Apr 2007 - 14:57:27 GMT

stopifnot(length(d) == 2, length(singValA) == min(d),

Ravi> Assistant Professor, The Center on Aging and Health .......

Ravi> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Apr 08 01:07:45 2007

Date: Sat 07 Apr 2007 - 14:57:27 GMT

>>>>> "Ravi" == Ravi Varadhan <rvaradhan@jhmi.edu>

>>>>> on Fri, 6 Apr 2007 12:44:33 -0400 writes:

Ravi> Hi, qr(A)$rank will work, but just be wary of the Ravi> tolerance parameter (default is 1.e-07), since the Ravi> rank computation could be sensitive to the tolerance Ravi> chosen.

rankMat <- function(A, tol = NULL, singValA = svd(A, 0,0)$d) {

## Purpose: rank of a matrix ``as Matlab'' ## ---------------------------------------------------------------------- ## Arguments: A: a numerical matrix, maybe non-square ## tol: numerical tolerance (compared to singular values) ## singValA: vector of non-increasing singular values of A ## (pass as argument if already known) ## ---------------------------------------------------------------------- ## Author: Martin Maechler, Date: 7 Apr 2007, 16:16d <- dim(A)

stopifnot(length(d) == 2, length(singValA) == min(d),

diff(singValA) < 0) # must be sorted decreasingly if(is.null(tol))

tol <- max(d) * .Machine$double.eps * abs(singValA[1])
else stopifnot(is.numeric(tol), tol >= 0)
sum(singValA >= tol)

}

A small scale simulation with random matrices, i.e., things like

## ranks of random matrices; here will have 5 all the time: table(replicate(1000, rankMat(matrix(rnorm(5*12),5, 12) )))# < 1 sec.

indicates that qr(.)$rank gives the same typically, where I assume one should really use

qr(., tol = .Machine$double.eps, LAPACK = TRUE)$rank

to be closer to Matlab's default tolerance.

Ok, who has time to investigate further? Research exercise:

1>> Is there a fixed number, say t0 <- 1e-15 1>> for which qr(A, tol = t0, LAPACK=TRUE)$rank is 1>> ``optimally close'' to rankMat(A) ?

2>> how easily do you get cases showing svd(.) to more reliable 2>> than qr(., LAPACK=TRUE)?

To solve this in an interesting way, you should probably
investigate classes of "almost rank-deficient" matrices,
and I'd also be interested if you "randomly" ever get matrices A
with rank(A) < min(dim(A)) - 1

(unless you construct some columns/rows exactly from earlier
ones or use all-0 ones)

Ravi> Ravi.

Ravi> -----------------------------------------------------------------Ravi> Ravi Varadhan, Ph.D.

Ravi> Assistant Professor, The Center on Aging and Health .......

Ravi> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

Ravi> -------------------------------------------------------------------- >> How about >> qr(A)$rank >> or perhaps >> qr(A, LAPACK=TRUE)$rank >> Cheers, >> Andy >> Hi! Maybe this is a silly question, but I need the >> column rank (http://en.wikipedia.org/wiki/Rank_matrix) >> of a matrix and R function 'rank()' only gives me the >> ordering of the elements of my matrix. How can I >> compute the column rank of a matrix? Is there not an R >> equivalent to Matlab's 'rank()'? I've been browsing>> for a time now and I can't find anything, so any help >> will be greatly appreciated. Best regards!

>> -- -- Jose Luis Aznarte M. >> http://decsai.ugr.es/~jlaznarte Department of Computer >> Science and Artificial Intelligence Universidad de >> Granada Tel. +34 - 958 - 24 04 67 GRANADA (Spain) Fax: >> +34 - 958 - 24 00 79 ______________________________________________R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Apr 08 01:07:45 2007

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sat 07 Apr 2007 - 23:31:17 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*