function R_opt = si_opt_r(T)
    arguments
        T (:,:) double
    end
    
    n = size(T, 2);
    R_opt = ones(n);
    
    for j = 1:n
        for k = 1:n
            R_opt(j, k) = ~any(all(T(:,[k, j]) == [0 1], 2));
        end
    end; clear j k;
    
    % Restore Symmetry
    for j = 1:n
        for k = 1:n
            R_opt(j, k) = R_opt(j, k) & R_opt(k, j);
        end
    end; clear j k;
    
    % Note: Could be done in one pass if processing upper / lower triag separately
end