Randomize non-diagonal elements of symmetric matrix
I have a symmetric matrix that I want to randomly shuffle while keeping the diagonal elements unchanged. The rows all sum to 1 and should still sum to 1 after shuffling.
Toy example below:
A <- rbind(c(0.6,0.1,0.3),c(0.1,0.6,0.3),c(0.1,0.3,0.6))
A
# [,1] [,2] [,3]
# [1,] 0.6 0.1 0.3
# [2,] 0.1 0.6 0.3
# [3,] 0.1 0.3 0.6
I would like a matrix B with the same diagonal elements as A and still symmetric, but with the elements randomly shuffled to generate something like
B <- rbind(c(0.6,0.3,0.1), c(0.3,0.6,0.1), c(0.3,0.1,0.6))
B
# [,1] [,2] [,3]
# [1,] 0.6 0.3 0.1
# [2,] 0.3 0.6 0.1
# [3,] 0.3 0.1 0.6
My aim is to do that on a 24 *24 matrix, so the code can be messy and no need to have something with a low computational cost. So far, I have tried with a loop but the code quickly gets excessively complicated and I was wondering whether there was a more straightforward way to do it.
Get indices of non-diagonal elements. Subset values and row indices. Within each row, shuffle values and assign back.
i = row(A) != col(A)
A[i] = ave(A[i], row(A)[i], FUN = sample)
A
# [,1] [,2] [,3]
# [1,] 0.6 0.1 0.3
# [2,] 0.3 0.6 0.1
# [3,] 0.3 0.1 0.6
If you don't want to overwrite the original matrix, assign to a copy instead.
A = rbind(c(0.6,0.1,0.3), c(0.1,0.6,0.3), c(0.1,0.3,0.6))
i = row(A) != col(A)
A2 = A
set.seed(1)
A2[i] = ave(A[i], row(A)[i], FUN = sample)
A2
# [,1] [,2] [,3]
# [1,] 0.6 0.1 0.3
# [2,] 0.1 0.6 0.3
# [3,] 0.3 0.1 0.6
set.seed(12)
A2[i] = ave(A[i], row(A)[i], FUN = sample)
A2
# [,1] [,2] [,3]
# [1,] 0.6 0.3 0.1
# [2,] 0.3 0.6 0.1
# [3,] 0.1 0.3 0.6
Since you want to keep rowwise sum as 1 you can only shuffle each row elements excluding diagonal elements in each row.
set.seed(2021)
t(sapply(seq(nrow(A)), function(x) {
tmp <- A[x, ]
tmp[-x] <- sample(tmp[-x])
tmp
}))
# [,1] [,2] [,3]
#[1,] 0.6 0.1 0.3
#[2,] 0.3 0.6 0.1
#[3,] 0.1 0.3 0.6
Try the code below
t(mapply(
function(x, k) replace(x, k, sample(x[k])),
asplit(A, 1),
asplit(row(A) != col(A), 1)
))
which gives
[,1] [,2] [,3]
[1,] 0.6 0.1 0.3
[2,] 0.3 0.6 0.1
[3,] 0.1 0.3 0.6
One option could be:
set.seed(123)
t(mapply(function(x, y) {
ind <- which(seq_along(x) != y)
`[<-`(x, ind, sample(x[ind]))
},
x = asplit(A, 1),
y = 1:nrow(A)))
[,1] [,2] [,3]
[1,] 0.6 0.1 0.3
[2,] 0.1 0.6 0.3
[3,] 0.1 0.3 0.6