

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:(n1)){
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);


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:(n1)){
> 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 1000by1000 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:(n1)){
x < diff(x) # use 'diff' in a loop
for(j in 1:(ni)){ # 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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:(n1)){
>> 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 1000by1000
> 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:(n1)){
> x < diff(x) # use 'diff' in a loop
> for(j in 1:(ni)){ # 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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:(n1)){
> > 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 1000by1000 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:(n1)){
> x < diff(x) # use 'diff' in a loop
> for(j in 1:(ni)){ # 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 ith column by (i1) 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


My apologies in advance for being a bit offtopic, 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) 5022619
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:(n1)){
> > 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 1000by1000 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:(n1)){
> x < diff(x) # use 'diff' in a loop
> for(j in 1:(ni)){ # 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 ith column by (i1) 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


On Apr 27, 2011, at 21:34 , Ravi Varadhan wrote:
> My apologies in advance for being a bit offtopic, 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/WronskianNot quite, I think. This is one function at different values of x, the Wronskian is about n different functions.
Tables of higherorder 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) 5022619
> 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:(n1)){
>>> 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 1000by1000 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:(n1)){
>> x < diff(x) # use 'diff' in a loop
>> for(j in 1:(ni)){ # 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 ith column by (i1) 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/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Peter, I have indeed worked with GregoryNewton 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) 5022619
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 offtopic, 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/WronskianNot quite, I think. This is one function at different values of x, the Wronskian is about n different functions.
Tables of higherorder 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) 5022619
> 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:(n1)){
>>> 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 1000by1000 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:(n1)){
>> x < diff(x) # use 'diff' in a loop
>> for(j in 1:(ni)){ # 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 ith column by (i1) 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/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, 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/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

