0

In R I'm interested in the general case to generate a matrix from a formula such as:

X = some other matrix
Y(i, j) = X(i, j) + Y(i - 1, j - 1)

Unfortunately I can't find how to account for the matrix self-referencing.

Obviously order of execution and bounds checking are factors here, but I imagine these could be accounted for by the matrix orientation and formula respetively.

Thanks.

c z
  • 7,726
  • 3
  • 46
  • 59

2 Answers2

0

Well, you can always use a for loop:

Y <- matrix(0, ncol=3, nrow=3)
#boundary values:
Y[1,] <- 1
Y[,1] <- 2

X <- matrix(1:9, ncol=3)

for (i in 2:nrow(Y)) {
  for (j in 2:ncol(Y)) {
    Y[i, j] <- X[i, j] + Y[i-1, j-1]
  }
}

If that is too slow you can translate it to C++ (using Rcpp) easily.

Roland
  • 127,288
  • 10
  • 191
  • 288
0

This solution assumes that you want Y[1,n] == X[1,n] and Y[n,1] == X[n,1]. If not, you can apply the same solution on the sub-matrix X[-1,-1] to fill in the values of Y[-1,-1]. It also assumes that the input matrix is square.

We use the fact that Y[N,N] = X[N,N] + X[N-1, N-1] + ... + X[1,1] plus similar relations for off-diagonal elements. Note that off-diagonal elements are a diagonal of a particular sub-matrix.

# Example input
X <- matrix(1:16, ncol=4)

Y <- matrix(0, ncol=ncol(X), nrow=nrow(X))

diag(Y) <- cumsum(diag(X))

Y[1,ncol(X)] <- X[1,ncol(X)]
Y[nrow(X),1] <- X[nrow(X),1]

for (i in 1:(nrow(X)-2)) {
  ind <- seq(i)
  diag(Y[-ind,]) <- cumsum(diag(X[-ind,])) # lower triangle
  diag(Y[,-ind]) <- cumsum(diag(X[,-ind])) # upper triangle
}
Matthew Lundberg
  • 42,009
  • 6
  • 90
  • 112