With tidyverse, you can arrive at this answers using usual group_by operations.

Sample data

I'm creating column names to make it easier to convert to tibble.

set.seed(123)
M <- matrix(sample(0:5, 100, TRUE),
            sample(0:5, 100, TRUE),
            ncol = 2,
            nrow = 100,
            dimnames = list(NULL, c("colA", "colB")))

Solution

library("tidyverse") 

as_tibble(M) %>% 
    arrange(colA, colB) %>% 
    group_by(colA, colB) %>% 
    summarise(num_pairs = n(), .groups = "drop") %>% 
    pivot_wider(names_from = colB, values_from = num_pairs) %>% 
    remove_rownames() 

Preview

# A tibble: 6 x 7
   colA   `0`   `1`   `2`   `4`   `5`   `3`
  <int> <int> <int> <int> <int> <int> <int>
1     0     4     4     4     2     4    NA
2     1     2     2     4     6     2    NA
3     2     6     4    NA     2     6    NA
4     3     2    NA    NA     4     6     2
5     4    NA     2     6    NA     2     4
6     5     6     2     4     4     2     2

Comments

You have asked:

I mean that I want to know how many pairs (0,0) , (0,1), ..., (0,5),... (5,5) are in both columns?

This answer gives you that, the question is how important is for you to have your results stored as a matrix? You can convert the results further into matrix by using as.matrix on what you get. Likely, I would stop after summarise(num_pairs = n(), .groups = "drop") as that gives very usable results, easy to subset join and so forth.

We can also use table

table(as.data.frame(M))

-output

# colB
#colA 0 1 2 3 4 5
#   0 4 4 4 0 2 4
#   1 2 2 4 0 6 2
#   2 6 4 0 0 2 6
#   3 2 0 0 2 4 6
#   4 0 2 6 4 0 2
#   5 6 2 4 2 4 2