Voting power (4): Speeding up the computation

Introduction and recapitulation

Recall from previous posts that we have considered two power indices for computing the power of a voter in a weighted system; that is, the ability of a voter to influence the outcome of a vote. Such systems occur when the voting body is made up of a number of "blocs": these may be political parties, countries, states, or any other groupings of people, and it is assumed that within every bloc all members will vote the same way. Indeed, in some legislatures, voting according to "the party line" is a requirement of membership.

Examples include the American Electoral College, in which the "voters" are the states; the Australian Senate, the European Union, the International Monetary Fund, and many others.

Given a set \(N=\{1,2,\ldots,n\}\) of voters and their weights \(w_i\), and a quota \(q\) required to pass any motion, we have represented this as

\[ [q;w_1,w_2,\ldots,w_n] \]

and we define a winning coalition as any subset \(S\subset N\) of voters for which

\[ \sum_{i\in S}w_i\ge q. \]

It is convenient to define a characteristic function \(v\) on all subsets of \(N\) so that \(v(S)=1\) if \(S\) is a winning coalition, and 0 otherwise. Given a winning coalition \(S\), a voter \(i\in S\) is necessary if \(v(S)-v(S-\{i\})=1\).

For any voter \(i\), the number of winning coalitions for which that voter is necessary is

\[ \beta_i = \sum_S\left(v(S)-v(S-\{i\})\right) \]

where the sum is taken over all winning coalitions. Then the Banzhaf power index is this value normalized so that the sum of all indices is unity:

\[ B_i=\frac{\beta_i}{\sum_{i=1}^n \beta_i}. \]

The Shapley-Shubik power index is defined by considering all permutations \(p\) of \(N\). Taken cumulative sums from the left, a voter \(p_k\) is pivotal if this is the first voter for which te cumulative sum is at least \(q\). For each voter \(i\), let \(\sigma_i\) be the number of permutations for which \(i\) is pivotal. Then

\[ S_i=\frac{1}{n!}\sigma_i \]

which ensures that the sum of all indices is unity.

Although these two indices seem very different, there is in fact a deep connection. Consider any permutation \(p\) and suppose that \(i=p_k\) is the pivotal voter, This voter will also be pivotal in all permutations for which \(i=p_k\) and the values to the right and left of \(i\) stay there: there will be \((k-1)!(n-k)!\) such permutations. However, we can consider the values up to and including \(i=p_k\) as a winning coalition for which \(i\) is necessary, which means we can write

\[ S_i=\sum_S\frac{(n-k)!(k-1)!}{n!}\left(v(S)-v(S-\{i\})\right) \]

which can be compared to the Banzhaf index above as being similar and with a different weighting function. Note that the above expression can be written as

\[ S_i=\sum_S\left(k{\dbinom n k}\right)^{-1}\left(v(S)-v(S-\{i\})\right) \]

which uses smaller numbers. For example, if \(n=50\) then \(n!\approx 3.0149\times 10^{64}\) but the largest binomial value is only approximately \(1.264\times 10^{14}\).

Computing with polynomials

We have also seen that if we define

\[ f_i(x) = \prod_{m\ne i}(1+x^{w_m}) \]

then

\[ \beta_i = \sum_{j=q-w_i}^{q-1}a_j \]

where \(a_j\) is the coefficient of \(x_j\) in \(f_i(x)\).

Similarly, if

\[ f_i(x,y) = \prod_{m\ne i}(1+yx^{w_m}) \]

then

\[ S_i=\sum_{k=0}^{n-1}\frac{k!(n-1-k)!}{n!}\sum_{j=q-w_i}^{q-1}c_{jk} \]

where \(c_{jk}\) is the coefficient of \(x^jy^k\) in the expansion of \(f_i(x,y)\).

In a previous post we have shown how to implement these in Python, using the Sympy library. However, Python can be slow, and using Cython is not trivial. We thus here show how to use Julia and its Polynomials package.

Using Julia for speed

The Banzhaf power indices can be computed almost directly from the definition:

 1using Polynomials
 2
 3function px(n)
 4    return(Polynomial([0,1])^n)
 5end
 6
 7function banzhaf(q,w)
 8    n = length(w)
 9    inds = vec(zeros(Int64,1,n))
10    for i in 1:n
11        p = Polynomial([1])
12        for j in 1:i-1
13            p = p * px(v[j])
14        end
15        for j in i+1:n
16            p = p*px(v[j])
17        end
18        inds[i] = sum(coeffs(p)[k] for k in q-w[i]+1:q)
19    end
20    return(inds)
21end

This will actually return the vector of \(\beta_i\) values which can then be easily normalized.

The function px is a "helper function" that simply returns the polynomial \(x^n\).

For the Shapley-Shubik indices, the situation is a bit trickier. There are indeed some Julia libraries for multivariate polynomials, but they seem (at the time of writing) to be not fully functional. However, consider the polynomials

\[ f_i(x,y)=\prod_{m\ne i}(1+yx^{w_m}) \]

from above. We can consider this as a polynomial of degree \(n-1\) in \(y\), whose coefficients are polynomials in \(x\). So if

\[ f_i(x,y) = 1 + p_1(x)y + p_2(x)y^2 +\cdots + p_{n-1}(x)y^{n-1} \]

then \(f_i(x,y)\) can be represented as a vector of polynomials

\[ [1,p_1(x),p_2(x),\ldots,p_{n-1}(x)]. \]

With this representation, we need to perform a multiplication by \(1+yx^p\) and also determine coefficients.

Multiplication is easy, noting at once that \(1+yx^p\) is linear in \(y\), and so we use the expansion of the product

\[ (1+ay)(1 + b_1y + b_2y^2 + \cdots + b_{n-1}y^{n-1}) \]

to

\[ 1 + (a+b_1)y + (ab_1+b_2)y^2 + \cdots + (ab_{n-2}+b_{n-1})y^{n-1} + ab_{n-1}y^n. \]

This can be readily programmed as:

1function mulp1(n,p)
2    p0 = Polynonial(0)
3    px = Polynomial([0,1])
4    c1 = cat(p,p0,dims=1)
5    c2 = cat([p0,p,dims=1)
6    return(c1 + px^n .* c2)
7end

The first two lines aren't really necessary, but they do make the procedure easier to read.

And we need a little program for extracting coefficients, with a result of zero if the power is greater than the degree of the polynomial (Julia's coeff function simply produces a list of all the coefficients in the polynomial.)

1function one_coeff(p,n)
2    d = degree(p)
3    pc = coeffs(p)
4    if n <= d
5        return(pc[n+1])
6    else
7        return(0)
8    end
9end

Now we can put all of this together in a function very similar to the Python function for computing the Shapley-Shubik indices with polynomials:

 1function shapley(q,w)
 2    n = length(w)
 3    inds = vec(zeros(Float64,1,n))
 4    for i in 1:n
 5        p = Polynomial(1)
 6        for j in 1:i-1
 7            p = mulp1(w[j],p)
 8        end
 9        for j in i+1:n
10            p = mulp1(w[j],p)
11        end
12        B = vec(zeros(Float64,1,n))
13        for j in 1:n
14            B[j] = sum(one_coeff(p[j],k) for k in q-w[i]:q-1)
15        end
16        inds[i] = sum(B[j+1]/binomial(n,j)/(n-j) for j in 0:n-1)
17    end
18    return(inds)
19end

And a quick test (with timing) of the powers of the states in the Electoral College; here ecv is the number of electors of all the states, in alphabetical order (the first states are Alabama, Alaska, Arizona, and the last states are West Virginia, Wisconsin, Washington):

1ecv = [9, 3, 11, 6, 55, 9, 7, 3, 3, 29, 16, 4, 4, 20, 11, 6, 6, 8, 8, 4, 10, 11, 16, 10, 6,
210, 3, 5, 6, 4, 14, 5, 29, 15, 3, 18, 7, 7, 20, 4, 9, 3, 11, 38, 6, 3, 13, 12, 5, 10, 3]
3@time(s = shapley(270,ecv));
40.722626 seconds (605.50 k allocations: 713.619 MiB, 7.95% gc time)

This is running on a Lenovo X1 Carbon, 3rd generation, using Julia 1.5.3. The operating system is a very recently upgraded version of Arch Linux, and currently using kernel 5.10.3.