matrix of higher order differences

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

matrix of higher order differences

Jeroen Ooms
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);

Reply | Threaded
Open this post in threaded view
|

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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: 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

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: matrix of higher order differences

Petr Savicky-2
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 <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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: matrix of higher order differences

Ravi Varadhan
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: [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 <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.

______________________________________________
[hidden email] 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.

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: matrix of higher order differences

Peter Dalgaard-2

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: [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 <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.
>
> ______________________________________________
> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

Re: matrix of higher order differences

Ravi Varadhan
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: [hidden email]


-----Original Message-----
From: peter dalgaard [mailto:[hidden email]]
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: [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 <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.
>
> ______________________________________________
> [hidden email] 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.
>
> ______________________________________________
> [hidden email] 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: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] 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.