[R] matrix of higher order differences

[Hans W Borchers]

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