Джулия или R: создать симметричную матрицу по минимуму недиагональных элементов - PullRequest
Новая
       39

Джулия или R: создать симметричную матрицу по минимуму недиагональных элементов

0 голосов
/ 23 марта 2019

Вот сложная проблема.У меня есть произвольная квадратная матрица в R (может быть и в Юлии), например 

> set.seed(420)
> A <- matrix(runif(16),nrow = 4,byrow = T)
> A
          [,1]      [,2]      [,3]      [,4]
[1,] 0.6055390 0.9702770 0.1744545 0.4757888
[2,] 0.7244812 0.8761027 0.3775037 0.6409362
[3,] 0.6546772 0.5062158 0.3033477 0.7162497
[4,] 0.2905202 0.1962252 0.3225786 0.8404279

Я хочу преобразовать эту матрицу в симметричную матрицу, чтобы недиагональные элементы всегда были минимальнымииз 2 соответствующих недиагональных элементов матрицы A.В приведенном выше случае результат должен быть следующим: 

          [,1]      [,2]      [,3]      [,4]
[1,] 0.6055390 0.7244812 0.1744545 0.2905202
[2,] 0.7244812 0.8761027 0.3775037 0.1962252
[3,] 0.1744545 0.3775037 0.3033477 0.3225786
[4,] 0.2905202 0.1962252 0.3225786 0.8404279

Хотелось бы получить эффективный и не трудоемкий способ его программирования в Юлии и / или R. Должен быть применим к любому виду квадратной матрицы.

Ответы [ 2 ]

2 голосов
/ 23 марта 2019

В Юлии однострочник - 

B = [min(A[i,j], A[j,i]) for i in axes(A, 1), j in axes(A, 2)]

Однако, это в два раза больше работы, чем необходимо. Следующий подход более эффективен, работает на месте, изменяя исходную матрицу, и генерирует приятное сообщение об ошибке, если ваша матрица не квадратная: 

function symmetrize_min!(A::AbstractMatrix)
    ax1 = axes(A, 1)
    axes(A, 2) == ax1 || error("A must be square")
    for j in ax1
        for i in j+1:last(ax1)
            A[i,j] = A[j,i] = min(A[i,j], A[j,i])
        end
    end
    return A
end

В Julia есть соглашение, что функции, изменяющие свои аргументы, заканчиваются на !, как предупреждение для пользователей. Если вы не хотите, чтобы это изменило A, то вы можете copy(A) перед вызовом или внести небольшое изменение: 

function symmetrize_min(A::AbstractMatrix)
    ax1 = axes(A, 1)
    axes(A, 2) == ax1 || error("A must be square")
    B = similar(A)
    for j in ax1
        B[j,j] = A[j,j]
        for i in j+1:last(ax1)
            B[i,j] = B[j,i] = min(A[i,j], A[j,i])
        end
    end
    return B
end

В Юлии вы должны обнять петли - они быстрые. И поскольку их легко написать, казалось бы, сложные проблемы становятся простыми.

2 голосов
/ 23 марта 2019

Вероятно, лучшим способом было бы создать пару циклов for с Rcpp, но вы также можете попробовать: 

idx <- A <= t(A)
A * idx + t(A) * (1L - idx)

Редактировать

Решение Base R неэффективно как по времени, так и по памяти. Если вам когда-нибудь понадобятся скорость и / или память, вот функции Rcpp, которые делают матрицу симметричной. Один изменяет его на месте, а другой возвращает копию. 

library(Rcpp)

cppFunction('
void inplaceSymmetric(NumericMatrix A) {
  for (int i = 1; i < A.nrow(); ++i)
    for (int j = 0; j < i; ++j)
      if (A(i, j) < A(j, i)) 
        A(j, i) = A(i, j);
      else 
        A(i, j) = A(j, i);
}')

cppFunction('
NumericMatrix copySymmetric(NumericMatrix A) {
  NumericMatrix C = clone(A);
  for (int i = 1; i < A.nrow(); ++i)
    for (int j = 0; j < i; ++j)
      if (A(i, j) < A(j, i))
        C(j, i) = C(i, j);
      else
        C(i, j) = C(j, i);
  return C;
}')
...