[R] matrix of higher order differences

[David Winsemius]

On Apr 27, 2011, at 7:25 AM, 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.

Every call to row() or col() creates a matrix of the same size as M.  
It might speed up if you created them outside the 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)
> }
> # ----

David Winsemius, MD
West Hartford, CT