matrix of higher order differences

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

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

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.

  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.