# matrix of higher order differences

7 messages
Open this post in threaded view
|

## matrix of higher order differences

 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);
Open this post in threaded view
|

## Re: matrix of higher order differences

 Jeroen Ooms 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) } # ---- ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: matrix of higher order differences

 On Apr 27, 2011, at 7:25 AM, Hans W Borchers wrote: > Jeroen Ooms 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: matrix of higher order differences

 In reply to this post by Hans W Borchers On Wed, Apr 27, 2011 at 11:25:42AM +0000, Hans W Borchers wrote: > Jeroen Ooms 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. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: 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/WronskianRavi. ------------------------------------------------------- Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: [hidden email] -----Original Message----- From: [hidden email] [mailto:[hidden email]] On Behalf Of Petr Savicky Sent: Wednesday, April 27, 2011 11:01 AM To: [hidden email] Subject: Re: [R] matrix of higher order differences On Wed, Apr 27, 2011 at 11:25:42AM +0000, Hans W Borchers wrote: > Jeroen Ooms 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. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code. ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.