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);
--
View this message in context: http://r.789695.n4.nabble.com/matrix-of-higher-order-differences-tp3477339p3477339.html
Sent from the R help mailing list archive at Nabble.com.
matrix of higher order differences
7 messages · Jeroen Ooms, Hans W Borchers, David Winsemius +3 more
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)
}
# ----
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
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)
}
# ----
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.
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
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)
}
# ----
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.
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)
}
# ----
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.
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
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)
}
# ----
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.
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