arXiv:1701.02201v4 [stat.CO] 23 Oct 2017

A fast algorithm for maximal propensity score
matching
Pavel S. Ruzankin∗1
1

Sobolev Institute of Mathematics, Novosibirsk, Russia
1
Novosibirsk State University, Novosibirsk, Russia

Abstract
We present a new algorithm which detects the maximal number of
matched disjoint pairs satisfying a given caliper when the matching
is done with respect to a scalar index (e.g., propensity score), and
constructs a corresponding matching. We assume the caliper to be a
Lipschitz function of the observations. If the observations are ordered
with respect to the index then the matching needs O(N ) operations,
where N is the total number of objects to be matched. The case of
1-to-n matching is also considered.
Keywords: propensity score matching, matching with caliper.

1

One-to-one matching

We consider matching disjoint pairs of objects from two groups, which we
will call, using common terminology, treated and control objects. In other
words, a control object can be matched to no more than one treated object
and vice versa. We will consider only one-dimensional distance, such as in
propensity score matching, when the distance between objects is the distance
between points on the real line corresponding to these objects, assuming each
object is somehow projected to a unique point on the real line. We will call
these points propensity scores of the objects for the sake of clarity. However
no assumptions are made on how these points are related to the objects.
Let Xi , i = 1, ..., K, and Yj , j = 1, ..., L, be the propensity scores of
treated and control objects, K and L being the total numbers of treated and
∗

email: ruzankin@math.nsc.ru

1

control objects, respectively. Xj and Yj may take any values on the real line,
not necessarily on the interval (0, 1). Let N = K + L. Let c = c(x, y) ≥ 0
be the caliper for our matching, i.e., we match only pairs (i, j) such that
|Xi − Yj | ≤ c(Xi , Yj ). We will assume that the caliper is Lipschitz in both
arguments with constants 1, i.e., for all x, y, t,
|c(x, y) − c(x + t, y)| ≤ |t|,
|c(x, y) − c(x, y + t)| ≤ |t|.

(1)
(2)

We will consider some less restrictive conditions on the caliper in Sec. 4.
A natural problem is to find the maximal number of pairs that can be
matched. Though this problem can be solved employing network flow optimization algorithms (e.g., see Hansen and Klopfer (2006)), the known algorithms have complexity not less than O(N 2 ) if no assumptions on sparsity are
made. This approach to matching problems was used, e.g., by Rosenbaum
(2012) and Pimentel et al. (2015).
Our main goal is to introduce a fast algorithm detecting the maximal
number of matched pairs and constructing a corresponding matching. Our
algorithm has complexity O(N) when both the treated and control objects
are sorted with respect to the propensity score:
X1 ≤ X2 ≤ · · · ≤ XK

and Y1 ≤ Y2 ≤ · · · ≤ YL .

(3)

Thus once we have sorted the observations (which takes O(N log N) or less
operations), we can reasonably fast solve the inverse problem of finding the
minimal constant caliper suitable for using q percent of data for a given q.
For instance, if the propensity score belongs to the interval (0, 1) then l runs
of the algorithm (O(lN) operations) yield the accuracy of 2−l for the minimal
caliper.
From now on we assume that relation (3) holds.
Let us now introduce the algorithm. The variable M will contain the
current number of matched pairs. After the algorithm finishes, M contains
the maximal number of matched pairs. Am and Bm store the index numbers
of treated and control object, respectively, in the m-th matched pair.
We present the algorithm as the following pseudocode:

2

M := 0
i := 1
j := 1
while (i ≤ K and j ≤ L)
if (|Xi − Yj | ≤ c(Xi , Yj ))
M := M + 1
AM := i
BM := j
i := i + 1
j := j + 1
else
if (Xi < Yj )
i := i + 1
else
j := j + 1
end if
end if
end while
As we see, the algorithm just walks through all the observations and successively collects all feasible pairs.
The algorithm requires O(N) operations since in each iteration of the
while-loop the variable i or j or both are increased. Certainly, to apply the
algorithm, first we must sort the observations with respect to the propensity
score, which requires O(N log N) operations or even less when modern radix
sort algorithms are used.
In Sec. 3 we prove that the algorithm produces the maximal possible
number of matched pairs.

2

1-to-n matching

The algorithm can be modified for 1-to-n matching. We assume that a treated
object is to be matched to no more than n control objects, and a control
object must not be matched to more than one treated object. Our algorithm
maximizes the number of matched control objects or, in other words, the
number of matched pairs.
The following pseudocode uses the same variables as above. Di is the
number of controls matched to the i-th treated object. The variable k corresponds to the current number of controls matched to the i-th treated object.

3

M := 0
i := 1
j := 1
k := 0
Di := 0 for all i = 1, ..., K
while (i ≤ K and j ≤ L)
if (|Xi − Yj | ≤ c(Xi , Yj ))
k := k + 1
M := M + 1
AM := i
BM := j
Di := k
if (k = n)
k := 0
i := i + 1
end if
j := j + 1
else
if (Xi < Yj )
k := 0
i := i + 1
else
j := j + 1
end if
end if
end while
The complexity is still O(N) and does not depend on n since, as above,
in each iteration of the while-loop the variable i or j or both are increased.

3

Optimality of the algorithm

We offer the following two proofs for the maximality of the number of pairs
matched by the algorithms above.

3.1

The first proof

First consider one-to-one matching. We will prove by induction that, under
conditions (1) and (2), the first algorithm yields the maximal number of
matched pairs.

4

There exists a matching M satisfying the caliper (i.e., |Xi −Yj | ≤ c(Xi , Yj )
for all (i, j) ∈ M) and containing the maximal number of matched pairs.
Consider the first step. If X1 < Y1 − c(X1 , Y1 ) then

X1 < Y1 − c(X1 , Y1 ) + Yj − Y1 + c(X1 , Y1 ) − c(X1 , Yj ) ≡ Yj − c(X1 , Yj ).
for all j, since Yj − Y1 + c(X1 , Y1) − c(X1 , Yj ) ≥ 0 by (2), and, hence, the
first treated object can not be used for matching. Analogously if Y1 < X1 −
c(X1 , Y1 ) then the first control object is not suitable for matching by (1).
Thus first steps of the algorithm skip the observations that can not be used
for matching.
After the above operation we can assume, for the sake of convenience, that
|X1 − Y1 | ≤ c(X1 , Y1 ). Let us show that matching now the first treated with
the first control object, as the algorithm does, does not reduce the maximal
number of matched pairs, if we match the maximal number of pairs for the
remaining 2, . . . , K-th treated and 2, . . . , L-th control objects.
If the first treated or the first control object are not matched in M then
removing from M a possible pair with the first treated or the first control
object and then adding (1, 1) to M does not change the number of pairs
in M. Thus, in this case, matching the pair (1, 1) and then matching the
maximal number of pairs for the 2, . . . , K-th treated and 2, . . . , L-th control
objects yields the total maximal number of matched pairs.
The case when M contains the pair (1, 1) is clear.
It remains to consider the case when M contains some pairs (1, j1 ) and
(i1 , 1), where i1 6= 1 and j1 6= 1. In this case we have Xi1 ≤ Y1 + c(Xi1 , Y1 )
and Yj1 ≤ X1 + c(X1 , Yj1 ). Therefore
Xi1 − Yj1 ≤ Y1 + c(Xi1 , Y1 ) − Yj1

= c(Xi1 , Yj1 ) − Yj1 − Y1 + c(Xi1 , Yj1 ) − c(Xi1 , Y1 )
≤ c(Xi1 , Yj1 )
by (2) and analogously Yj1 − Xi1 ≤ c(Xi1 , Yj1 ) by (1). Hence,
|Xi1 − Yj1 | ≤ c(Xi1 , Yj1 ).
Thus, removing from M the pairs (1, j1 ) and (i1 , 1) and adding the pairs (1, 1)
and (i1 , j1 ) does not change the number of pairs in M. Again, matching the
pair (1, 1) and then matching the maximal number of pairs for the 2, . . . , Kth treated and 2, . . . , L-th control objects yields the total maximal number
of matched pairs.
Applying the above argument to the remaining 2, . . . , K-th treated and
2, . . . , L-th control observations proves the validity of the first algorithm by
induction.
5

For the case of 1-to-n matching it suffices to consider the second algorithm as the first one applied to observations where we take n identical
treated objects instead of each corresponding treated object from the original
observations, i.e., we “repeat” each treated object n times.

3.2

The second proof

For the second proof we use a theorem on the Monge-Kantorovich mass
transfer problem. Though this proof is suitable only for the case of constant
caliper c = const, it illustrates the connection between matching and the
mass transfer problem.
Let P and Q be finite (nonnegative) continuous measures on the real line
with Borel σ-algebra with P (R) = Q(R). Let
ρ(P, Q) = inf {U({(x, y) : |x − y| > c})},
U

where the supremum is taken over all (nonnegative) continuous measures on
R2 with marginals P and Q: U(A, R) = P (A), U(R, A) = Q(A) for all Borel
A. Put
F (t) = P((−∞, t)), G(t) = P((−∞, t)).
Ruzankin (2001) proved the following theorem.
Theorem A. The equalities hold:
ρ(P, Q) =
=

lim S(y) − P (R)

y→∞

lim T (y) − Q(R),

y→∞

where the functions S and T are specified by the relations
lim S(y) = lim T (y) = 0,

y→−∞

y→−∞

dS(y) = max{dF (y), T (y + dy − c) − S(y)},
dT (y) = max{dG(y), S(y + dy − c) − T (y)}

(4)

(5)
(6)

for all y and the assumption that the functions S and T are left continuous.
The functions S and T exist and are uniquely defined.
The relations (5), (6) and the left-continuity condition mean that
S(y + w) = S(y) + P ([y, y + w))
+ sup (T (y + v − z) − S(y) − P ([y, y + v)))+ ,
0<v≤w

6

(7)

T (y + w) = T (y) + Q([y, y + w))
+ sup (S(y + v − z) − T (y) − Q([y, y + v)))+

(8)

0<v≤w

for all y and w > 0, where t+ = max{t, 0}.
Put
µ(P, Q) = sup{U({(x, y) : |x − y| ≤ c})},

(9)

U

where the supremum is taken over all continuous measures on R2 with marginals
P and Q. We have
µ(P, Q) = P (R) − ρ(P, Q) = 2P (R) − lim S(y) = 2Q(R) − lim T (y).
y→∞

y→∞

The measure U0 that achieves the supremum in (9) can be built as follows.
Put
V (y) = F (y) − T (y + c) + G(y + c),
W (y) = G(y) − S(y + c) + F (y + c).

(10)
(11)

The functions V, W, F − V, G − W are left continuous and nondecreasing.
Let the measure Z on R2 be defined by
Z((−∞, x) × (−∞, y)) = min{V (x), W (y)}.

(12)

Let the measure R be an arbitrary measure with marginals F − V and G − V
(here we use a function of y instead of the corresponding measure of (−∞, y)).
Put U0 = Z + R. Then the measure U0 achieves the supremum in (9).
Here the measure Z is responsible for an optimal “mass transfer”:
Z({(x, y) : |x − y| ≤ c}) = Z(R2 ) = µ(P, Q)
since V (y − c) ≤ W (y) ≤ V (y + c) for all y.
To prove the optimality of the above algorithms we are to consider the
case of discrete P and Q with supports consisting of finite numbers of points.
First consider one to one matching. Let the measure P̃ be concentrated
at the points X1 ≤ X2 ≤ · · · ≤ XK and the measure Q̃ be concentrated on
the points Y1 ≤ Y2 ≤ · · · ≤ YL with
P̃ ({y}) = #{j : Xj = y},

Q̃({y}) = #{j : Yj = y},

where the # sign denotes the number of elements of a set.

7

If K 6= L then to use Theorem A we have to extend one of the measures
P̃ or Q̃. Put
D = max{XK , YL } + 2c,
P (A) = P̃ (A) + (Q(R) − P (R))+ ID (A),
Q(A) = Q̃(A) + (P (R) − Q(R))+ ID (A),
where ID (A) = 1 if D ∈ A and is zero otherwise. Now P (R) = Q(R) and we
can apply Theorem A.
The function S can increase only at the points Xj , Yj + c, D + 2c, D + 3c
while the function T can increase only at the points Yj , Xj +c, D +2c, D +3c.
Relations (7), (8) can be rewritten as
S(y + 0) = max{T (y − c + 0), S(y) + P ({y})},
T (y + 0) = max{S(y − c + 0), T (y) + Q({y})}.

(13)
(14)

We can consider the problem of maximal one-to-one matching of the
points Xi to the points Yj within the caliper c as the problem of achieving
the supremum in (9), where U({(x, y)}) = k ≥ 0 for |x − y| ≤ c means that
exactly k points Xi = x are matched to k points Yj = y. Relations (10)–(11)
mean that S(y) − F (y) counts the “spare” part of Q((−∞, y − c)), which can
not be matched. Analogously T (y) − G(y) counts the part of P ((−∞, y − c))
that is not matched. On the other hand, relations (13)–(14) mean that the
matching, which corresponds to the measure Z, is done according to the
algorithms above. Since we can not match more than µ(P, Q) pairs, the
former algorithm yields the maximal number of matched pairs.
For the case of 1-to-n matching we can put
P̃ ({y}) = n #{j : Xj = y},

Q̃({y}) = #{j : Yj = y}.

and repeat the above argument.

4

The case of piecewise Lipschitz caliper

In this section we will describe an algorithm which yields a maximal number
of pairs under somewhat weaker conditions on the caliper. We will consider
one-to-one matching though it is easy to modify the algorithm below for the
case of 1-to-n matching just like it was done for the first algorithm.
We will assume that, first, the caliper is “Lipschitz-nondecreasing”:
c(x, y + t) ≥ c(x, y) − t for all t > 0, x, y,
c(x + t, y) ≥ c(x, y) − t for all t > 0, x, y
8

(15)
(16)

and, second, the caliper is piecewise Lipschitz in both arguments: there exist
disjoint intervals [a1 , a2 ), ..., [aU −1 , aU ) covering the domain of Xi and disjoint
intervals [b1 , b2 ), ..., [bV −1 , bV ) covering the domain of Yj such that, for each
u = 1, ..., U − 1,
c(x + t, y) ≤ c(x, y) + t for all x, x + t, y ∈ [au , au+1 ), t > 0

(17)

and, for each v = 1, ..., V − 1,
c(x, y + t) ≤ c(x, y) + t for all x, y, y + t ∈ [bv , bv+1 ), t > 0.

(18)

For example, if c(x, y) = f (x)+g(y), where f (x) and g(y) are nondecreasing step functions, or if c(x, y) = e−|x| (1 − (y − ⌊y⌋)), where ⌊y⌋ denotes the
greatest integer not greater than y, then conditions (15)–(18) are satisfied.
Let us now introduce an algorithm for a caliper satisfying (15)–(18). As
above, M is the current number of matched pairs. After the algorithm finishes, M is the maximal number of matched pairs. Am and Bm store the index numbers of treated and control object, respectively, in the m-th matched
pair.
I1 , ..., IU are increasing numbers such that Xi ∈ [au , au+1 ) whenever Iu ≤
i < Iu+1 ; and increasing numbers J1 , ..., JV are such that Yj ∈ [bv , bv+1 )
whenever Jv ≤ j < Jv+1 . Computing I1 , ..., IU and J1 , ..., JV given a1 , ..., aU ,
b1 , ..., bV , X1 , ..., XK , and Y1 , ..., YL requires O(N) operations. If some of the
intervals [au , au+1 ) contain no observations Xi then we are to take the number
of intervals [Iu , Iu+1 ) lesser than the number of intervals [au , au+1 ), but, to
simplify notations, we use the same U to enumerate Iu , u = 1, ..., U. The
same is done for the intervals [Jv , Jv+1 ).
For each u = 1, ..., U − 1, we will have Su = i if and only if i ∈ [Iu , Iu+1],
the observations XIu , ..., Xi−1 are already matched or discarded, and either
i = Iu+1 or Xi is currently neither matched nor discarded. Symmetrically,
for each v = 1, ..., V − 1, we will have Tv = j if and only if j ∈ [Jv , Jv+1 ],
the observations YJv , ..., Yj−1 are already matched or discarded, and either
j = Iu+1 or Yj is currently neither matched nor discarded.
We will assume that “and” in the if-statement means that the second
condition is checked only if the first one is true.

9

M := 0
Su := Iu for all u = 1, ..., U
Tv := Jv for all v = 1, ..., V
i := 1
j := 1
u0 := 1
v0 := 1
function increment i()
Su0 := Su0 + 1
i := i + 1
if (i = Iu0+1 )
u1 := u0 + 1
while (u1 < U and Su1 = Iu1+1 ) u1 := u1 + 1
u0 := u1
i := Su0
end if
end function
function increment j ()
Tv0 := Tv0 + 1
j := j + 1
if (j = Jv0+1 )
v1 := v0 + 1
while (v1 < V and Tv1 = Jv1+1 ) v1 := v1 + 1
v0 := v1
j := Tv0
end if
end function

10

while (i ≤ K and j ≤ L)
if (Xi < Yj )
for (v = v0, ..., V − 1)
if (Tv < Jv+1 and |Xi − YTv | ≤ c(Xi , YTv ))
M := M + 1
AM := i
BM := Tv
increment i()
if (v0 = v)
increment j ()
else
Tv := Tv + 1
end if
next while
end if
end for
increment i()
else
for (u = u0, ..., U − 1)
if (Su < Iu+1 and |XSu − Yj | ≤ c(XSu , Yj ))
M := M + 1
AM := Su
BM := j
increment j ()
if (u0 = u)
increment i()
else
Su := Su + 1
end if
next while
end if
end for
increment j ()
end if
end while
The complexity of the last algorithm is O((U + V )N) since each iteration
of the while-loop requires O(U + V ) operations and in each iteration the
variable i or j or both are increased.
The proof for the maximality of the number of matched pairs almost
repeats the proof for the first algorithm. The main difference that if, say,
at some step Xi < Yj then we have to check sequentially whether Xi can
11

be matched to each group {YJv , ..., YJv+1−1 }, v = 1, ..., V − 1, of the control
observations. As above, by (18) it is sufficient for each group to check whether
Xi can be matched to the first unmatched element of the group. Relations
(15) and (16) ensure that matching Xi to the first unmatched element of the
first suitable group does not diminish the number of matched pairs below its
maximal value.

Acknowledgments
The author is grateful to Ben B. Hansen for useful remarks helping to improve
the paper. The algorithms are planned to be added to the Optmatch R
package developed by Ben B. Hansen, Mark Fredrickson et al.

References
1. Hansen, B. B. and Klopfer, S. O. (2006), “Optimal full matching and
related designs via network flows”, Journal of Computational and Graphical Statistics, 15, No.3, 609–627.
2. Pimentel, S. D., Kelz, R. R., Silber, J. H, and Rosenbaum, P. R. (2015),
“Large, Sparse Optimal Matching With Refined Covariate Balance in
an Observational Study of the Health Outcomes Produced by New Surgeons,” Journal of the American Statistical Association, 110, No. 510,
517–527.
3. Rosenbaum, P. R. (2012), “Optimal Matching of an Optimally Chosen
Subset in Observational Studies,” Journal of Computational and Graphical Statistics, 21, No. 1, 57–71.
4. Ruzankin, P. S. (2001), “Construction of the optimal joint distribution
of two random variables,” Theory Probab. Appl., 46, No.2, 316–334.

12

