Matching Soulmates in R


In Matching Soulmates, my co-authors and I study a a recursive process in matching that forms coalitions that are mutually most-preferred by members of that coalition. We call this process the iterated matching of soulmates: IMS. We show that mechanisms that implement IMS have strong properties among those who are matched by IMS.

How many people can “usually” be matched by IMS? This depends a lot on the environment and structure of preferences. But what about in totally unstructured environments?

Here is some R code that takes advantage of R's array-centric functions and some custom operators to count the number of people that can be matched as soulmates in 10000 random 10-person stable roommates problems.

The core of this code is contained in the whos_a_soulmate function. This function is somewhat non-standard R code. It is inspired by the type of programming normally done in one of R's predecessors APL. Here is the function:

whos_a_soulmate <- function(p){which((p%r%rank%>%hadamard%r%min)==1)}

To see how this works, let's first setup an example preference matrix.

p <- matrix(c(3,1,2,1,3,2,2,1,3),3,3,byrow=TRUE)
##      [,1] [,2] [,3]
## [1,]    3    1    2
## [2,]    1    3    2
## [3,]    2    1    3

Plyer 1 likes 2 best, player 2 likes 1 best, player 3 likes 2 best. Note that 1 and 2 are soulmates. We should find that 1 and 2 are soulmates in the first round of IMS.

The first step of the function ensures the matrix “p” is a ranking matrix. We apply the “rank” function to the “p” matrix by row using the by-row operator created here. Since “p” is already a ranking matrix, this does not do anything. It is more important once some players have been removed and we need to re-rank.

p %r% rank
##      [,1] [,2] [,3]
## [1,]    3    1    2
## [2,]    1    3    2
## [3,]    2    1    3

Now we take take Hadamard product of the preference matrix with its own transpose by piping it to our “hadamard” function with the built-in R pipe.

p %r% rank %>% hadamard
##      [,1] [,2] [,3]
## [1,]    9    1    4
## [2,]    1    9    2
## [3,]    4    2    9

Note how off diagonal 1's in this matrix represent positions where both players prefer eachother most. We want to pull these players who have a 1 in their row. To do this, we pipe what we have to the “min” function using our by-row operator. If a 1 is present, that player has a soulmate.

p %r% rank %>% hadamard %r% min
##      [,1] [,2] [,3]
## [1,]    1    1    2

Note that 1 and 2 have a “1” here, indicating that players 1 and 2 have a soulmate from this group. Select just these players, we compare this vector to 1. The indicies that map to “TRUE” are players who have a soulmate.

p %r% rank %>% hadamard %r% min == 1
##      [,1] [,2]  [,3]

We see now that players 1 and 2 have a soulmate. We can now remove these players by using our remove operator “%rm%”.

p <- p %rm% whos_a_soulmate(p) %>% as.matrix
##      [,1]
## [1,]    3

What we get is the remainder of the preference matrix “p” after removing soulmates. We only have player 3 left. Note that player 3 originally ranked being “alone” as third-best. That's why a 3 shows up here. Let's pipe this to rank again to have the player(s) rerank their potential partners.

p %r% rank
##      [,1]
## [1,]    1

Now we can run the process again.

p <- p %rm% whos_a_soulmate(p)
## <0 x 0 matrix>

Since being alone is the only outcome left, player 3 is a solemate-group-of-one and is removed. This demonstrates how the “IMS” function proceeds. It calls “remove_soulmates”. If anyone is removed, “IMS” calls itself. If there is no one to remove or no one further can be removed, the function returns the remaining preference matrix.


#Operator "m %r% f" applies function "f" to matrix "m" by rows. 
`%r%` <- function(m,f){t(apply(m,1,}

#Operator "m %rm% i" removes rows and column indicies "i" from matrix "m". 
`%rm%` <- function(m,i){if(length(i)>0){m[-i,-i]}else{m}}  

#Operator applies function "f" across list "l" with sapply.
`%s%` <- function(l,f){sapply(l,}

#Returns Hadamard product of matrix m with its transpose.
hadamard <- function(m){m * t(m)}

#Given preference matrix, returns vector of player indicies who are first-order soulmates.
whos_a_soulmate <- function(p){which((p%r%rank%>%hadamard%r%min)==1)}

#Removes first order soulmates from preference matrix p.
remove_soulmates <- function(p){p %rm% whos_a_soulmate(p) %>% as.matrix}

#Recursively applies remove_soulmates until there's none left to remove.
ims <- function(p){if(dim(p)[1]==0 || identical(remove_soulmates(p),p)){p}else{ims(remove_soulmates(p))}}

#Create a random roommates preference matrix with "dim" players.
create_preference <- function(dim){
p <- matrix(runif(dim*dim),dim,dim,byrow=TRUE)
diag(p)<- 1
p %r% rank

#Set up parameters. "dim" is the number of players. "n" is the number of random trials. 
dim <- 10
n <- 10000

#Create List of "n" Preference Matricies with "dim" players.
preference_list <- lapply(1:n,function(x){create_preference(dim)})

#Apply IMS to the list of preferences then count the number of players remaining.
remaining <- preference_list %s% ims %s% dim 

#Get the number of players removed by IMS.
removed <- dim - remaining[1,]

#Make a table of frequencies of the number of removed players.
removed %>% table/n
## .
##      0      2      4      6     10 
## 0.5365 0.2794 0.1079 0.0323 0.0439