Julia or R: create symmetric matrix by minimum of off-diagonal elements

Question

Here is a complex problem. I have an arbitrary square matrix in R (can be in Julia as well), for example

> 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

I want to transform this matrix in a symmetric matrix, so that the off-diagonal elements are always the minimum of the 2 corresponding off-diagonal elements of the A matrix. In the above case, the result must be:

          [,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

Would like an efficient and non-time consuming way of programming it in Julia and/or R. Must be applicable to any kind of square matrix.

josemz · Answer 1 · 2019-04-25T15:47:20.467

3

Probably the best way would be to make a couple of for loops with Rcpp, but you could give this a try as well:

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

Edit

The base R solution is inefficient, both time and memorywise. If you ever need the speed and/or memory, here are the Rcpp functions that make the matrix symmetric. One modifies it in-place and the other returns a copy.

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;
}')

edited Apr 25 '19 at 15:47

answered Mar 22 '19 at 23:02

josemz

1,283
7
15

Will try the second option because I'm not a C++ user (and therefore a RCpp user) ;-) Which language is the second option? RCpp? – coolsv Mar 23 '19 at 04:18
I've edited my answer to include the `Rcpp` functions. – josemz Mar 26 '19 at 16:11
The R base solution is brilliant and worked like a charm with my data – Javier Fajardo Aug 31 '20 at 16:39

score 2 · Accepted Answer · answered Mar 23 '19 at 07:02

In Julia, a one-liner is

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

However, this does twice as much work as needed. The following approach is more efficient, works in-place modifying the original matrix, and generates a nice error message if your matrix isn't square:

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

In Julia it's a convention that functions modifying their arguments end with !, as a warning to users. If you don't want this to modify A, then either you can copy(A) before calling it or make a small change:

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

In Julia you should embrace loops---they are fast. And because they are easy to write, seemingly complex problems become simple.

Another one-liner is `min.(A, transpose(A))`, which, as you note does does (a bit more than) twice the necessary work. Oddly, the in-place version: `A .= min.(A, transpose(A))` allocates as much as `B = min.(A, transpose(A))`. — DNF, Mar 23 '19 at 07:13
I never understood the `!` option in Julia, but with your explanation, it seems more clear ;-) — coolsv, Mar 23 '19 at 16:29

Julia or R: create symmetric matrix by minimum of off-diagonal elements

2 Answers2