How to use recursive formula in loop and if statement

I may be overcomplicating this problem.

The problem is as such

Variables

a = lower triangular matrix/dataframe 20x20 
b = a 1x20 matrix/vector 
c = the previous row result of the formula (recursive bit)

I want to have a lower triangular matrix (I,j) where it is recursively defined that

pseudocode

if(I<j) = 0

else if (I==j) = 1

else (I>j) = sum(a * b * c )/ b[I] - b[J] (where I is current row position and J is current column position)

Example of how I hope the formula/ R’s vectorisation will work shown on a smaller matrix to make life easier.

I = row

j = column

b(I,j)

b(1,2) refers to position in matrix

Pseudo code of example matrix

    Matrix C
         1                                                                       |  2                                                                    |3  |4  |5  |6  |
    1|i=j=1                                                                      |i<j=0                                                                  |...|...|...|...| 
    2|i>j=a(2,1)*b(1)*c(1,1)/b(2)-b(1)                                           |i=j=1                                                                  |...|...|...|...|
    3|i>j= (a(3,1)*b(1)*c(1,1) + a(3,2)*b(2)*c(2,1))/(b(3)-b(1)                  |i>j= a(3,2)*b(2)*c(2,2)/(b(3)-b(2)                                     |...|...|...|...|
    4|i>j = a(4,1)*b(1)*c(1,1)+a(4,2)*b(2)*c(2,1)+a(4,3)*b(3)*c(3,1))/b(4)-b(1)  |i>j= a(4,2)*b(2)*c(2,2) + a(4,3)*b(3)*c(3,2)/b(4)-b(2)                 |...|...|...|...|
    5|...                                                                        |...                                                                    |...|...|...|...|
    6|...                                                                        |...                                                                    |...|...|...|...|

Onto the code I have so far,

Create variable first as shown below:

    my_int <- 20
    
    nr <- as.integer(my_int)
    
    #create a n x n matrix with zeroes 
    
    a <- matrix(0, nr, nr)
    
    # For each row and for each column, assign values based on position
    
    # These values are the product of two indexes
    
    for(i in 1:dim(a)[1]) {
      for(j in 1:dim(a)[2]) {
        a[i,j] = if(i<j) {
          0
        }else if(i==j) {
          1
        }else {
          3
        }
        }
        }
    
    # make into dataframe
    
    mymat <- data.frame(mymat)
    
    # create b variable
    
    z <- rep(1:20)
    
    # made it into lower diagonal matrix to make it easier to work with
    
    b <- matrix(0, length(z), length(z))
    b[lower.tri(b, diag = TRUE)] <- z[sequence(length(z):1)]
    b
    
    # create matrix for formula to operate with 
    
    # create variable "c"
    
    
    c <- matrix(0, nr, nr)
    
    # For each row and for each column, assign values based on position
    # These values are the product of two indexes
    for(i in 1:dim(c)[1]) {
      for(j in 1:dim(c)[2]) {
        mymat2[i,j] = if(i<j) {
          0
        }else if(i==j) {
          1
        }else {
          5 # place holder for now
        }
        }
        }

Formula for calculating results, I think works thanks to R’s vectorisation

    sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,])

My issue is then how to insert this in the if statements for the creation of the recursively defined matrix as shown below

    # calculate recursively defined lower triangular matrix
    
    c <- matrix(0, nr, nr)
    
    # For each row and for each column, assign values based on position
    # These values are the product of two indexes
    for(i in 1:dim(c)[1]) {
      for(j in 1:dim(c)[2]) {
        mymat2[i,j] = if(i<j) {
          0
        }else if(i==j) {
          1
        }else {
          sum(a*b*lag(c), na.rm = TRUE)/(b[,j]-b[I,]) # formula for calculation of values for lower triangular matrix
        }
        }
        }

This gives error

Error in mymat2[i, j] <- if (i < j) { : number of items to replace is not a multiple of replacement length

I can link an excel spreadsheet where this formula works if it would help. It is only achievable in excel by a lot of manual inputs etc.

Example of expected results using real data

a = 5x5 lower triangle matrix

0|0   |0   |0   |0
5|0   |0   |0   |0
5|0.56|0   |0   |0
5|0.20|0.61|0   |0
5|0.06|0.16|0.61|0


b = 1x5 matrix/vector

0.27917|0.499|0.83|1.191|1.48

c = recursive matrix results

1   |0   |0   |0   |0
6.36|1   |0   |0   |0
5.77|0.84|1   |0   |0
5.50|0.77|1.43|1   |0
5.3 |0.72|1.80|2.46|1

Example of how I calculated "c" in excel

1   |0   |0   |0   |0

=(INDEX(a,2,1)*INDEX(b,1)*INDEX(c,1,1))/(INDEX(b,2)-INDEX(b,1))|1   |0   |0   |0

=(INDEX(a,3,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,1))/(INDEX(b,3)-INDEX(b,1))|=(INDEX(a,3,2)*INDEX(b,2)*INDEX(c,2,2))/(INDEX(b,3)-INDEX(b,2))|1   |0   |0

=(INDEX(a,4,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,1))/(INDEX(b,4)-INDEX(b,1))|=(INDEX(a,4,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,2))/(INDEX(b,4)-INDEX(b,2))|=(INDEX(a,4,3)*INDEX(b,3)*INDEX(c,3,3))/(INDEX(b,4)-INDEX(b,3))|1   |0

=(INDEX(a,5,1)*INDEX(b,1)*INDEX(c,1,1)+INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,1)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,1)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,1))/(INDEX(b,5)-INDEX(b,1)) |=(INDEX(a,5,2)*INDEX(b,2)*INDEX(c,2,2)+INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,2)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,2))/(INDEX(b,5)-INDEX(b,2))|=(INDEX(a,5,3)*INDEX(b,3)*INDEX(c,3,3)+INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,3))/(INDEX(b,5)-INDEX(b,3))|=(INDEX(a,5,4)*INDEX(b,4)*INDEX(c,4,4))/(INDEX(b,5)-INDEX(b,4))|1

Above code shows how it can be calculated in excel via index and manually setting a large amount of cell positions.

# FUNCTIONS
vpad <- function(vec) { 
    extend <- 20 - length(vec)
    c(vec, rep(0, extend))
}

clip <- function(x) { x[x >= 0 & x <= 20] }

pad_a <- function(a, iter) {
    mat <- cbind(a[, iter:20], matrix(0, 20, iter-1))
    mat * lower.tri(mat)
}

pad_b <- function(b, iter) {
    tmp <- vpad(b[iter:20])
    mat <- matrix(rep(tmp, each=20), 20, 20)
    mat * lower.tri(mat)
}

pad_c <- function(c, iter) {
    L <- 21 - iter
    i <- clip(seq(L) + iter - 1)
    j <- seq(i)
    v <- vpad(mapply(function(i, j) c[i, j], i, j))
    mat <- matrix(rep(v, each=20), 20, 20)
    mat <- mat * (lower.tri(mat) + diag(20))
    mat

    end <- 21 - iter
    mat1 <- rbind(matrix(0, iter - 1, 20), mat[1:end, ])
    mat1 * lower.tri(mat1)
}

# MAKE FAKE DATA
set.seed(1)
m <- matrix(sample(400)/400, 20, 20)
a <- m * lower.tri(m)
b <- sample(20)/21

# INITIALIZE VARS
add_prev <- matrix(rep(0, 200), 20, 20)
init_mat <- lower.tri(m) + diag(20)
overwrite_mat <- init_mat

for (i in 1:20) {
    overwrite_mat <- (pad_a(a, i) * pad_b(b, i) * pad_c(overwrite_mat, i)) + add_prev
    add_prev <- overwrite_mat
}


denom1 <- matrix(rep(b, times=20), 20, 20)
denom1 <- denom1 * lower.tri(denom1) + diag(20)
denom2 <- matrix(rep(b, each=20), 20, 20)
denom2 <- denom2 * lower.tri(denom2) + diag(20)

final <- ((overwrite_mat / denom1) - denom2) + diag(20)
final <- replace(final, is.nan(final), 0)

 [,1]         [,2]        [,3]        [,4]        [,5]       
 [1,]   0.0000000   0.00000000  0.00000000  0.00000000   
 [2,]  -0.5129073   0.00000000  0.00000000  0.00000000   
 [3,]   0.6738060  -0.05926190  0.00000000  0.00000000   
 [4,]   9.2856557   5.97660268 -0.15059524  0.00000000   
 [5,]   0.1971807   0.01591345 -0.14775033 -0.08718254