From: Ravi Varadhan <rvaradhan_at_jhmi.edu>

Date: Wed, 27 Apr 2011 17:41:22 -0400

Ravi Varadhan, Ph.D.

Assistant Professor,

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

*>
*

*> Hi.
*

*>
*

> The following avoids the inner loop and it was faster

*> for x of length 100 and 1000.
*

*>
*

*> diffmatrix2 <- function(x){
*

*> 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.
*

Date: Wed, 27 Apr 2011 17:41:22 -0400

Peter, I have indeed worked with Gregory-Newton and divided differences in my very first numerical analysis course a couple of decades ago! However, I am perplexed by the particular form of this matrix where the differences are stored along the diagonals. I know that this is not the *same* as the Wronskian, but was just wondering whether it is an established matrix that is some kind of an *ian* like Hermitian, Jacobian, Hessian, Wronskian, Laplacian, ...

Best,

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: peter dalgaard [mailto:pdalgd_at_gmail.com]
Sent: Wednesday, April 27, 2011 4:59 PM

To: Ravi Varadhan

Cc: R Help

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

On Apr 27, 2011, at 21:34 , Ravi Varadhan wrote:

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

*>
**> What might one do with a matrix of all order finite differences? It seems that such a matrix might be related to the Wronskian (its discrete analogue, perhaps).
**>
**> http://en.wikipedia.org/wiki/Wronskian
*

Not quite, I think. This is one function at different values of x, the Wronskian is about n different functions.

Tables of higher-order differences were used fundamentally for interpolation and error detection in tables of function values (remember those?), but rarely computed to the full extent - usually only until the effects of truncation set in and the differences start alternating in sign.

*>
*

> 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) >> } >> # ----

> The following avoids the inner loop and it was faster

-- Peter Dalgaard Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes_at_cbs.dk Priv: PDalgd_at_gmail.com ______________________________________________ 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 - 21:55:27 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 - 22:00:36 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.
*