Simple question but not an obvious answer! The solution I suggest follows a slightly different strategy than yours and does not involve any for loops. First, I developed a function (i.e. area_cover_()) that produces a cross table with only those polygons that have at least one intersection. Then, in a second step, I developed another function (i.e. add_isolated_poly()) that adds the polygons with no intersection at the end of the cross table produced in the first step. This makes it easier to read the final table if you have many polygons without intersection. So, please find the reprex below.

NB: The input data for the reprex corresponds to your sf object s.f.2 with area column

Reprex

1. First step: create a cross table including only the polygons that have at least one intersection (not including the polygons without intersection makes the reading of the cross table more efficient). To do this, I developed the function area_cover()

• Code of the area_cover() function

library(sf)
library(dplyr)
library(tidyr)


area_cover <- function(x) {
  x %>% 
  st_intersection() %>% 
  filter(n.overlaps>1) %>%
  mutate(area_inter = st_area(.)) %>%  
  unnest(., cols = c(origins, geometry)) %>% 
  left_join(., as.data.frame(x), by = c("origins" = "id")) %>% 
  mutate(cover_percent = area_inter/area.y*100) %>% 
  select(.,origins, area.y, area_inter, cover_percent) %>% 
  rename("id" = "origins", "area" = "area.y") %>% 
  st_drop_geometry() %>%  
  group_by(area_inter) %>% 
  mutate(poly_X_id = rev(id)) %>% 
  relocate(poly_X_id, .before = area_inter) %>% 
  xtabs(cover_percent ~ id + poly_X_id, data = .) %>%  
  replace(.== 0, 100) %>% 
  round(., digits = 1)
  }

• Output of the area_cover() function

Results <- area_cover(s.f.2)

Results      # Only polygons with at least one intersection are present in this cross table
#>    poly_X_id
#> id      1     2
#>   1 100.0  31.6
#>   2  27.0 100.0  

class(Results) 
#> [1] "xtabs" "table"  # you can convert 'Results' into matrix with 'as.matrix()' if needed.

2. Second step (optional): add isolated polygons (i.e. with no intersection) at the end of the cross table 'Results' (i.e. result of the previous step). To do this, I developed the function add_isolated_poly() which creates a dataframe with n columns corresponding to the n id's of the isolated polygons and filled with 0s

• Code of the add_isolated_poly() function

add_isolated_poly <- function(y, z){ # 'y arg.' = s.f.2 and 'z arg.' = result of the function area_cover()

id_isolated_poly <- setdiff(y$id, colnames(z))

df_isolated_poly <- y %>% 
  filter(.,id %in% id_isolated_poly) %>% 
  st_drop_geometry() %>% 
  select(., id) %>% 
  t() %>% 
  as.data.frame() %>% 
  `colnames<-`(., id_isolated_poly) %>% 
  rbind(.,rep(list(rep(0, nrow(y))), length(id_isolated_poly))) %>% 
  slice(-c(1))

cbind.fill <- function(...){
  nm <- list(...)
  nm <- lapply(nm, as.matrix)
  n <- max(sapply(nm, nrow))
  do.call(cbind, lapply(nm, function (x)
    rbind(x, matrix(0, n-nrow(x), ncol(x)))))
}

Results %>% 
  cbind.fill(., df_isolated_poly) %>%  
  replace(., col(.) == row(.), 100) %>%  
  `rownames<-`(., c(colnames(z), id_isolated_poly))
}

• Output of the add_isolated_poly() function

Results_2 <- add_isolated_poly(s.f.2, Results)

Results_2
#>     1     2   3
#> 1 100  31.6   0
#> 2  27 100.0   0
#> 3   0   0.0 100

class(Results_2)
#> [1] "matrix" "array

^{Created on 2021-11-19 by the reprex package (v2.0.1)}

IMPORTANT EDIT

Although they produce the right result with the proposed minimal example, the two functions I proposed above are not generalizable and produce wrong results. After a lot of trial and error, here is a simple, very fast function... and, this time, right!!! So, please find the reprex below.

Solution

• Code of the area_cover() function

library(sf)
library(dplyr)

area_cover <- function(x){
  Results <- x %>% 
    st_intersection(.,.) %>% 
    mutate(area_inter = st_area(.),
           cover = area_inter/area.1*100) %>% 
    st_drop_geometry() %>% 
    xtabs(cover ~ id.1 + id, data = ., sparse = TRUE) %>%  
    round(., digits = 1) %>% 
    as.matrix(.) 
  
  names(dimnames(Results)) <- NULL
  
  return(Results)
}

• Output of the area_cover() function

area_cover(s.f.2)

#>     1     2   3
#> 1 100  31.6   0
#> 2  27 100.0   0
#> 3   0   0.0 100

^{Created on 2021-11-22 by the reprex package (v2.0.1)}

Benchmarking

I compared my function with the validated solution of @wibeasley based on the new dataset provided by @M. Beausoleil (see comments below).

For the comparison to be valid, I have slightly modified the @wibeasley's function so that the output is a matrix containing percentages (i.e. same output as my function)

• Data

set.seed(1234)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 100
l = vector("list", n)
for (i in 1:n)
  l[[i]] = p + 7 * runif(2)
s = st_sfc(l)
s.f = st_sf(s)
s.f$id = c(1,1,2,2,3,4,5,5,6,7,7,8,9,8,7,3,4,5,5,6)
s.f.2 = s.f %>% group_by(id) %>%   summarise(geometry = sf::st_union(s)) # %>% summarise(area = st_area(s))
s.f.2$area = st_area(s.f.2)

• Code to compare the two solutions

library(sf)
library(dplyr)
library(tidyr)
library(bench)

bench_results <- bench::mark(
  "wibeasley_validated_answer" = {  #!!!!!  NB: lighlty modified function for the comparison to be valid!!!!!
    wibeasley <- sf::st_intersection(s.f.2, s.f.2) %>% 
      dplyr::mutate(
        area       = sf::st_area(.),
        proportion = area / area.1 * 100
      ) %>%
      tibble::as_tibble() %>%
      dplyr::select(
        id_1 = id,
        id_2 = id.1,
        proportion,
      ) %>% 
      # tidyr::complete(id_1, id_2, fill = list(proportion = 0))
      tidyr::pivot_wider(
        names_from = id_1,
        values_from = proportion,
        values_fill = 0
      ) %>% 
      as.matrix(.,rownames.force = TRUE) %>% 
      `<-`(., .[,-c(1)]) %>%  
      round(.,1)
    },  
  "lovalery_answer" = {
    lovalery <- area_cover(s.f.2)
  },
  min_iterations = 1000,
  relative = TRUE, 
  check = TRUE)

• Results of the benchmarking (relative values)

bench_results
#> # A tibble: 2 x 6
#>   expression                   min median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>                 <dbl>  <dbl>     <dbl>     <dbl>    <dbl>
#> 1 wibeasley_validated_answer  1.13   1.10      1          1       1.16
#> 2 lovalery_answer             1      1         1.10      24.1     1

• Results of the benchmarking (absolute values)

bench_results
#> # A tibble: 2 x 6
#>   expression                      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>                 <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 wibeasley_validated_answer   46.4ms   49.1ms      20.2    2.54MB    1.11 
#> 2 lovalery_answer              41.5ms   43.7ms      22.7   61.09MB    0.969

• Final checking that the two functions produce the same result

all.equal(wibeasley, lovalery)
#> [1] TRUE

wibeasley
#>       1     2     3     4     5     6     7     8     9
#> 1 100.0  19.0  13.3   4.6  22.6  21.2  12.8   3.7  11.6
#> 2  18.3 100.0  28.7  31.9  33.0  14.3  32.2  25.1   5.1
#> 3  14.6  32.7 100.0  23.3  35.5   7.2  28.8   7.0  13.7
#> 4   5.1  36.4  23.3 100.0  20.3  23.7  26.2  26.8  14.7
#> 5  13.0  19.8  18.7  10.7 100.0  10.3  27.0  15.6  12.4
#> 6  21.3  15.0   6.6  21.8  18.0 100.0  26.9  24.4  15.8
#> 7   9.5  24.8  19.5  17.7  34.7  19.9 100.0  11.8  17.9
#> 8   3.8  26.8   6.6  25.1  27.7  24.9  16.3 100.0   1.3
#> 9  22.1  10.0  23.7  25.5  41.0  30.0  46.0   2.3 100.0

lovalery
#>       1     2     3     4     5     6     7     8     9
#> 1 100.0  19.0  13.3   4.6  22.6  21.2  12.8   3.7  11.6
#> 2  18.3 100.0  28.7  31.9  33.0  14.3  32.2  25.1   5.1
#> 3  14.6  32.7 100.0  23.3  35.5   7.2  28.8   7.0  13.7
#> 4   5.1  36.4  23.3 100.0  20.3  23.7  26.2  26.8  14.7
#> 5  13.0  19.8  18.7  10.7 100.0  10.3  27.0  15.6  12.4
#> 6  21.3  15.0   6.6  21.8  18.0 100.0  26.9  24.4  15.8
#> 7   9.5  24.8  19.5  17.7  34.7  19.9 100.0  11.8  17.9
#> 8   3.8  26.8   6.6  25.1  27.7  24.9  16.3 100.0   1.3
#> 9  22.1  10.0  23.7  25.5  41.0  30.0  46.0   2.3 100.0

^{Created on 2021-11-22 by the reprex package (v2.0.1)}

There are a few things I'd optimize. I can't really evaluate the optimization with only three polygons, and I'm guessing the intersection calculation is the only real expensive part, so I'd start there.

You'll get an instant ~50% reduction if you don't calculate both (1) polygon A & B, and then (2) polygon B & A. In a sense, calculate only the upper triangle, and reuse/reflect the values to the lower triangle.

# I think you wanted to create 5 empty columns.  Use `numeric(0)` instead of `NULL`
df.fill=data.frame(id1=numeric(0), id2=numeric(0), area=numeric(0), over1=numeric(0), over2=numeric(0))

for (k in seq_along(id.sort)) {
  for (op in seq(from = k, to = length(id.sort), by = 1)) { # Avoid the lower triangle
    int.out = st_intersection(
      data.sp[data.sp$id==id.sort[k],], 
      data.sp[data.sp$id==id.sort[op],]
    )
    
    if(nrow(int.out) != 0) {
      area.tmp = st_area(int.out)#/10000
      over1 = area.tmp/int.out$area.ha
      over2 = area.tmp/int.out$area.ha.1 
    } else {area.tmp = 0;over1=0;over2=0}
    
    df.fill.tmp.upper = data.frame(id1 = id.sort[k], id2=id.sort[op], 
                                   area = area.tmp,
                                   over1 = over1,
                                   over2 = over2)
    df.fill.tmp.lower = data.frame(id1 = id.sort[op], id2=id.sort[k], 
                                   area = area.tmp,
                                   over1 = over2,
                                   over2 = over1)
    df.fill <- 
      if (k == op) rbind(df.fill, df.fill.tmp.upper)
      else         rbind(df.fill, df.fill.tmp.upper, df.fill.tmp.lower)
  }
}
df.fill %>% 
  dplyr::mutate(
    over1 = as.numeric(over1) * 100
    over2 = as.numeric(over2) * 100
  ) %>%
  select(-area, -over2) %>% 
  pivot_wider(
    names_from = id2,
    values_from = over1, 
    values_fill = 0
  )

Output:

# A tibble: 3 x 4
    id1   `1`   `2`   `3`
  <dbl> <dbl> <dbl> <dbl>
1     1 100    31.6     0
2     2  27.0 100       0
3     3   0     0     100

Without a real example to benchmark, I'm not sure it's faster than your solution. But it's simpler and easier to understand (at least for my brain).

sf::st_intersection() is vectorized. So it will find & return all the intersections of the first & second argument for you. In this case, the two arguments are the same set of polygons.

sf::st_intersection(s.f.2, s.f.2) %>% 
  dplyr::mutate(
    area       = sf::st_area(.),
    proportion = area / area.1
  ) %>%
  tibble::as_tibble() %>%
  dplyr::select(
    id_1 = id,
    id_2 = id.1,
    proportion,
  ) %>% 
  # tidyr::complete(id_1, id_2, fill = list(proportion = 0))
  tidyr::pivot_wider(
    names_from = id_1,
    values_from = proportion,
    values_fill = 0
  )

Output:

# A tibble: 3 x 4
   id_2   `1`   `2`   `3`
  <dbl> <dbl> <dbl> <dbl>
1     1 1     0.316     0
2     2 0.270 1         0
3     3 0     0         1

Things to consider:

• Keep the areas as proportions, instead of percentages. It's usually better for calculations later.
• Stay long, and don't pivot. It's usually better for calculations later, because you can join on id_1 and id_2.
• If you pivot wide, you might want it as a Matrix, instead of a data.frame.

If the intersection calculation is really expensive compared to the rest of the code, this solution might be faster than my second one. Like @lovalery's solution, only one data.frame is passed to sf:st_intersection(). By my interpretation, the main difference is this solution stays long and uses a cross join to enumerate all the possible combinations, and pivots wide in the final step. It's easier (for me) to manipulate when things are long instead of wide.

intersections <-
  s.f.2 %>% 
  sf::st_intersection() %>% 
  dplyr::filter(2L <= n.overlaps) %>%
  dplyr::mutate(
    numerator = sf::st_area(.)
  ) %>%
  tibble::as_tibble() %>%
  tidyr::unnest_wider(origins) %>%
  dplyr::select(
    id_1 = `...1`, 
    id_2 = `...2`, 
    numerator
  )

denominators <- 
  s.f.2 %>% 
  tibble::as_tibble() %>% 
  dplyr::select(id, area) 
  
denominators %>% 
  dplyr::full_join(denominators, by = character()) %>% 
  dplyr::select(
    id_1 = id.x,
    id_2 = id.y,
    area = area.y
  ) %>% 
  dplyr::left_join(intersections, by = c("id_1" = "id_1",  "id_2" = "id_2")) %>%
  dplyr::left_join(intersections, by = c("id_1" = "id_2",  "id_2" = "id_1")) %>% 
  dplyr::mutate(
    numerator   = dplyr::coalesce(numerator.x, numerator.y, 0),
    proportion  = dplyr::if_else(id_1 == id_2, 1, numerator / area),
  ) %>% 
  dplyr::select(id_1, id_2, proportion) %>% 
  tidyr::pivot_wider(
    names_from  = id_1,
    values_from = proportion,
    values_fill = 0
  )

Output:

# A tibble: 3 × 4
   id_2   `1`   `2`   `3`
  <dbl> <dbl> <dbl> <dbl>
1     1 1     0.316     0
2     2 0.270 1         0
3     3 0     0         1