# 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:

```
using Polynomials
function px(n)
return(Polynomial([0,1])^n)
end
function banzhaf(q,w)
n = length(w)
inds = vec(zeros(Int64,1,n))
for i in 1:n
p = Polynomial([1])
for j in 1:i-1
p = p * px(v[j])
end
for j in i+1:n
p = p*px(v[j])
end
inds[i] = sum(coeffs(p)[k] for k in q-w[i]+1:q)
end
return(inds)
end
```

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:

```
function mulp1(n,p)
p0 = Polynonial(0)
px = Polynomial([0,1])
c1 = cat(p,p0,dims=1)
c2 = cat([p0,p,dims=1)
return(c1 + px^n .* c2)
end
```

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.)

```
function one_coeff(p,n)
d = degree(p)
pc = coeffs(p)
if n <= d
return(pc[n+1])
else
return(0)
end
end
```

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:

```
function shapley(q,w)
n = length(w)
inds = vec(zeros(Float64,1,n))
for i in 1:n
p = Polynomial(1)
for j in 1:i-1
p = mulp1(w[j],p)
end
for j in i+1:n
p = mulp1(w[j],p)
end
B = vec(zeros(Float64,1,n))
for j in 1:n
B[j] = sum(one_coeff(p[j],k) for k in q-w[i]:q-1)
end
inds[i] = sum(B[j+1]/binomial(n,j)/(n-j) for j in 0:n-1)
end
return(inds)
end
```

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):

```
ecv = [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,
10, 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]
@time(s = shapley(270,ecv));
0.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.