Как получить абсолютную разницу между значениями в двух столбцах матрицы

У меня есть матрица, подобная следующей

      i j value
 [1,] 3 6 0.194201129
 [2,] 3 5 0.164547043
 [3,] 3 4 0.107149279
 [4,] 4 3 0.004927017
 [5,] 3 1 0.080454448
 [6,] 1 2 0.003220612
 [7,] 2 6 0.162313646
 [8,] 3 3 0.114992628
 [9,] 4 1 0.015337253
[10,] 1 6 0.026550051
[11,] 3 2 0.057004116
[12,] 4 2 0.006441224
[13,] 4 5 0.025641026
[14,] 2 4 0.004885993
[15,] 1 1 0.036552785
[16,] 1 5 0.048249186
[17,] 1 4 0.006053565
[18,] 1 3 0.004970296

Как видите, для некоторых пар i, j существует обратная пара. Например, для i = 3, j = 1 есть пара с i = 1, j = 3.

Вот чего я хочу добиться.

Для каждой пары i, j нужно вычесть обратное значение и получить абсолютное значение вычитания. Для тех пар, у которых нет обратной пары, из них вычитается 0.

Вот несколько примеров:

Для i = 3, j = 5 обратной пары нет (i = 5, j = 3), поэтому расчет становится следующим:

abs(0.164547043 - 0)

Для i = 3, j = 1 в матрице есть обратная пара с i = 1, j = 3, поэтому расчет будет следующим:

abs(0.004970296 - 0.080454448)

Я подошел к этому, написав кучу кода (65 строк), полного циклов for, и его трудно читать и редактировать.

Поэтому мне было интересно, есть ли другой более эффективный способ сделать что-то подобное, используя более компактные функции.

Мотивированный предыдущим сообщением, где его ответ был довольно простым (с использованием функции агрегата()) и поиском этих функций в Интернете, я пытаюсь использовать здесь mapply(), но правда в том, что я не могу справиться с обратным пары.

РЕДАКТИРОВАТЬ:

dput()
    memMatrix <- structure(c(3, 3, 3, 4, 3, 1, 2, 3, 4, 1, 3, 4, 4, 2, 1, 1, 1, 
        1, 6, 5, 4, 3, 1, 2, 6, 3, 1, 6, 2, 2, 5, 4, 1, 5, 4, 3, 0.194201128983738, 
        0.164547043451226, 0.107149278958536, 0.00492701677834917, 0.0804544476798398, 
        0.00322061191626409, 0.162313646044361, 0.114992627755601, 0.0153372534398016, 
        0.0265500506171091, 0.0570041160347523, 0.00644122383252818, 
        0.0256410256410256, 0.00488599348534202, 0.0365527853282693, 
        0.0482491856677524, 0.0060535654765406, 0.00497029586494912), .Dim = c(18L, 
        3L), .Dimnames = list(NULL, c("i", "j", "value")))

Также вот код, который пока работает, но он намного сложнее

Где memMatrix — матрица, приведенная вверху поста. И здесь вы можете увидеть небольшую разницу в том, что я умножаю абсолютное значение на переменную с именем probability_distribution, но это не имеет большого значения. Я удалил это (умножение) из исходного поста, чтобы сделать его более простым.

subFunc <- function( memMatrix , probability_distribution )
{

  # Node specific edge relevance matrix
  node_edgeRelm <- matrix(ncol = 3)
  colnames(node_edgeRelm) <- c("i","j","rel")
  node_edgeRelm <- na.omit(node_edgeRelm)

  for ( row in 1:nrow( memMatrix ) )
  {
    pair_i <- memMatrix[row,"i"]
    pair_j <- memMatrix[row,"j"]

    # If already this pair of i and j has been calculated continue with the next pair
    # At the end of a new calculation, we store the i,j (verse) values in order from lower to higher
    # and then we check here for the inverse j,i values (if exists).
    if( pair_i < pair_j )
      if( any(node_edgeRelm[,"i"] == pair_i & node_edgeRelm[,"j"] == pair_j) ) next
    if( pair_j < pair_i )
      if( any(node_edgeRelm[,"i"] == pair_j & node_edgeRelm[,"j"] == pair_i) ) next

    # Verse i,j
    mepm_ij <- as.numeric( memMatrix[which( memMatrix[,"i"] == pair_i & memMatrix[,"j"] == pair_j ), "mep"] )
    if( length(mepm_ij) == 0 )
      mepm_ij <- 0
    # Inverse j,i
    mepm_ji <- as.numeric( memMatrix[which( memMatrix[,"i"] == pair_j & memMatrix[,"j"] == pair_i ), "mep"] )
    if( length(mepm_ji) == 0 )
      mepm_ji <- 0

    # Calculate the edge relevance for that specific initial node x and pair i,j
    edge_relevance <- probability_distribution * abs( mepm_ij - mepm_ji )

    # Store that specific edge relevance with an order from lower to higher node
    if ( pair_i < pair_j)
      node_edgeRelm <- rbind( node_edgeRelm, c( as.numeric(pair_i), as.numeric(pair_j), as.numeric(edge_relevance) ) )
    else
      node_edgeRelm <- rbind( node_edgeRelm, c( as.numeric(pair_j), as.numeric(pair_i), as.numeric(edge_relevance) ) )
  }

  na.omit(node_edgeRelm)
}

вы можете запустить его как subFunc(memMatrix, 1/3)