From: Ravi Varadhan <rvaradhan_at_jhmi.edu>

Date: Wed, 27 Apr 2011 15:34:09 -0400

Ravi Varadhan, Ph.D.

Assistant Professor,

Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University

}

R-help_at_r-project.org 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.

R-help_at_r-project.org 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 Wed 27 Apr 2011 - 19:37:48 GMT

Date: Wed, 27 Apr 2011 15:34:09 -0400

My apologies in advance for being a bit off-topic, but I could not quell my curiosity.

http://en.wikipedia.org/wiki/Wronskian

Ravi.

Ravi Varadhan, Ph.D.

Assistant Professor,

Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University

Ph. (410) 502-2619

email: rvaradhan_at_jhmi.edu

-----Original Message-----

From: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] On Behalf Of Petr Savicky
Sent: Wednesday, April 27, 2011 11:01 AM
To: r-help_at_r-project.org

Subject: Re: [R] matrix of higher order differences

On Wed, Apr 27, 2011 at 11:25:42AM +0000, Hans W Borchers wrote:

> Jeroen Ooms <jeroenooms <at> gmail.com> writes:

*>
**> >
**> > Is there an easy way to turn a vector of length n into an n by n matrix, in
**> > which the diagonal equals the vector, the first off diagonal equals the
**> > first order differences, the second... etc. I.e. to do this more
**> > efficiently:
**> >
**> > diffmatrix <- function(x){
**> > n <- length(x);
**> > M <- diag(x);
**> > for(i in 1:(n-1)){
**> > differences <- diff(x, dif=i);
**> > for(j in 1:length(differences)){
**> > M[j,i+j] <- differences[j]
**> > }
**> > }
**> > M[lower.tri(M)] <- t(M)[lower.tri(M)];
**> > return(M);
**> > }
**> >
**> > x <- c(1,2,3,5,7,11,13,17,19);
**> > diffmatrix(x);
**> >
**>
**> I do not know whether you will call the appended version more elegant,
**> but at least it is much faster -- up to ten times for length(x) = 1000,
**> i.e. less than 2 secs for generating and filling a 1000-by-1000 matrix.
**> I also considered col(), row() indexing:
**>
**> M[col(M) == row(M) + k] <- x
**>
**> Surprisingly (for me), this makes it even slower than your version with
**> a double 'for' loop.
**>
**> -- Hans Werner
**>
**> # ----
**> diffmatrix <- function(x){
**> n <- length(x)
**> if (n == 1) return(x)
**>
**> M <- diag(x)
**> for(i in 1:(n-1)){
**> x <- diff(x) # use 'diff' in a loop
**> for(j in 1:(n-i)){ # length is known
**> M[j, i+j] <- x[j] # and reuse x
**> }
**> }
**> M[lower.tri(M)] <- t(M)[lower.tri(M)]
**> return(M)
**> }
**> # ----
*

Hi.

The following avoids the inner loop and it was faster for x of length 100 and 1000.

n <- length(x) if (n == 1) return(x) A <- matrix(nrow=n+1, ncol=n) for(i in 1:n){ A[i, seq.int(along=x)] <- x x <- diff(x) } M <- matrix(A, nrow=n, ncol=n) M[upper.tri(M)] <- t(M)[upper.tri(M)] return(M)

}

Reorganizing an (n+1) x n matrix into an n x n matrix shifts i-th column by (i-1) downwards. In particular, the first row becomes the main diagonal. The initial part of each of the remaining rows becomes a diagonal starting at the first component of the original row.

Petr Savicky.

R-help_at_r-project.org 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.

R-help_at_r-project.org 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 Wed 27 Apr 2011 - 19:37:48 GMT

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

Archive generated by hypermail 2.2.0, at Wed 27 Apr 2011 - 21:10:34 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.
*