Why am I getting different answer for same algorithm on R studio desktop vs R studio cloud?

Question

I run this simple code from class on my R studio desktop(M1 pro MacBook Pro) and the answer is different to the answer I get on R studio cloud. On the MacBook, I get 239.3093 vs 237.7821 on cloud. Am I missing something?

A = matrix(NA,nrow = 50,ncol = 50)
for (i in 1:50) {
  for (j in 1:50) {
    A[i,j] = sin(i) + cos(j)
  }
}

pivot.above = function(A, r, c) {
  A[r, ] = A[r, ]/A[r,c]
  for (i in (r-1):1) {
    A[i, ] = A[i, ]-A[i,c]*A[r, ]
  }
   A
}

for (i in 0:48) {
  j = 50-i
  A = pivot.above(A,j,j)
}

A[1, ] = A[1, ]/A[1,1]

print(sum(A))

Answer 1

This matrix has a very large condition number. From Wikipedia:

If the condition number is very large, then the matrix is said to be ill-conditioned. Practically, such a matrix is almost singular, and the computation of its inverse, or solution of a linear system of equations is prone to large numerical errors.

kappa(A)
[1] 1.352794e+19

This means that tiny numerical differences, such as the difference between compilers, system linear algebra libraries, operating systems, chips, etc., will lead to non-negligible differences in the results — as you found out.

As @DaveArmstrong commented, what the OP is doing here is not exactly either of the operations mentioned in the Wikipedia article (inversion or solution of a linear system); we're reducing the matrix to lower triangular form by Gaussian elimination. However, they're closely related (once we've reduced the matrix to triangular form, we've done most of the work of solving a linear system, because we can now use simple back-substitution). I don't know offhand why the sum of the first column of the reduced matrix should be a quantity that is especially sensitive, but it's not surprising to me that the condition number is related to this sensitivity. In these notes, Luke Tierney (an R-core member ...) describes why simple Gaussian elimination is unstable, and why a more sophisticated partial pivoting approach is more stable.

This is a perfectly reasonable question to ask, but is also something you're going to have to learn about/get used to as you move ahead with numerical computation. It would actually be a great thing to ask your instructor about.

Answer 2

Apologies that this isn't really an answer, but perhaps it's somewhere to start. In response to Ben's response to my question, I understand that the algorithm isn't just adding numbers, but so far as I can see, it is just vectorized scalar arithmetic (addition, subtraction, multiplication, division) - there are no matrix operations going on. The fact that the inputs and outputs are organized in a matrix is incidental. Consider what I think is an equivalent set of operations for data organized in a rectangular data frame rather than a square matrix.

B <- expand.grid(
  row = 1:50, 
  col = 1:50)

B <- B %>% 
  mutate(A = sin(row) + cos(col))

I rewrote the pivot.above function to work for this long-form data:

pivot.above.long <- function(A, rowind, colind){
  A$A[which(A$row == rowind)]  <- A$A[which(A$row == rowind)]/A$A[which(A$row == rowind & A$col == colind)]
  for(i in (rowind - 1):1){
    A$A[which(A$row == i)] = A$A[which(A$row == i)]-A$A[which(A$row == i & A$col == colind)]*A$A[which(A$row == rowind)]
  }
  A
}

Then, I applied that to the long-form data:

for (i in 0:48) {
  j = 50-i
  B = pivot.above.long(B,j,j)
}

B$A[which(B$row == 1)] <- B$A[which(B$row == 1)]/B$A[which(B$row == 1 & B$col == 1)]

And, I get exactly the same result as I got from the matrix operations - one of the two values returned in the original question.

print(sum(B$A))
# [1] 239.3093