---
title: "Power measures"
author: "Massimo Franceschet"
output:
  ioslides_presentation:
    css: ../style.css
    incremental: yes
editor_options: 
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(cache = TRUE, message = FALSE, warning = FALSE)
```

## Overview

* in some circumstances centrality - the quality of being connected to central ones - has limited utility in 
predicting the locus of power in networks
* consider exchange networks, where the relationship in the network involves the transfer of valued items (i.e., information, time, money, energy) 
* a set of exchange relations is positive if exchange in one relation promotes exchange in others and negative if exchange in one relation inhibits exchange in others
* in **negative exchange networks**, power comes from being connected to those who have few options. Being connected to those who have many possibilities reduces one's power

## Overview

* think, for instance, to a social network in which time is the exchanged value
* imagine that every actor has a limited time to listen to others and that each actor divides its time between its neighbors.
* clearly, exchange of time in one relation precludes the exchange of the same time in other relations 
* what are the actors that receive most attention?
* these are the nodes that are connected to many neighbors with few options, since they receive almost full attention from all their neighbors
* on the other hand, actors connected to few neighbors with a lot of possibilities receive little consideration, since their neighbors are mostly busy with others


## Power

The power thesis is as follows:

>A node is powerful if it is connected with powerless nodes.  


## Math

Power $x_{i}$ of node $i$ is given by:

$$x_i =  \sum_k \frac{a_{k,i}}{x_k}$$

that is:

$$x = x^{\div} A$$ 

The above states two properties of power:

>1. the larger the degree of a node, the larger its power. Hence, the more ties an actor has, the more powerful the actor is
>2. the smaller the power of neighbors of a node, the larger its power. For equal number of ties, actors that are linked to powerless others are powerful; on the other hand, actors that are tied to powerful others are powerless. 

## Math

Power can be computed with the following **iterative process**. 

Let $x_0 = e$, where $e$ is the vector of all 1's. For $k \geq 1$:

$$
x_k  = x_{k-1}^{\div} A
$$

For $k \geq 0$, let $r_k = x_{2k}$ and $c_k = x_{2k+1}$. 

>* If $r_k$ converges, then it does to a vector $r$ such that $\alpha r = r^{\div} A$, with $\alpha > 0$. It follows that $\sqrt{\alpha} r = (\sqrt{\alpha} r)^{\div} A$, hence $\sqrt{\alpha} r$ is a solution of the power equation $x  = x^{\div} A$. 
>* Similarly, if $c_k$ converges, it does to a vector $c$ such that $\beta c = c^{\div} A$ and $\sqrt{\beta} c$ is a again solution of the power equation $x  = x^{\div} A$. 

It holds that solutions $r$ and $c$ are proportional: $r = \gamma c$, for some  $\gamma > 0$.

## Math

* the 1st iteration $$x_1 = e A$$ is, for a given node $i$, the degree of $i$
* the 2nd iteration $$x_2  = (e A)^{\div} A$$ is, for a given node $i$, the sum of reciprocals of degrees of neighbors of $i$
* in general, power of node $i$ is positively influenced by the degrees of nodes that reach $i$ with an even number of edges  
* it is negatively influenced by the degrees of nodes that reach $i$ with an odd number of edges 
* nodes that reach $i$ with paths of even length as well as with paths of odd length play a dual role with respect to $i$
* notice, however, that such an influence devalues with the length of the path: it is stronger when the path is shorter

## Math

By [Sinkhorn-Knopp theorem](https://msp.org/pjm/1967/21-2/pjm-v21-n2-p14-s.pdf) (on symmetric matrices): 

> Power exists, and the corresponding iterative method converges, if and only if the graph $G_A$ associated with matrix $A$ is regularizable. 

The regularizability problem is equivalent to the following linear programming feasibility problem:

$$
\begin{array}{l}
B w = re \\
w \geq 1
\end{array}
$$

where: 

>* $B$ is the incidence matrix of $G_A$, that is, an $n \times m$ matrix such that $b_{i,l} = 1$ if $i$ belongs to edge $l$ and $b_{i,l} = 0$ otherwise, with $n$ and $m$ the number of nodes and edges of the graph
>* $w$ is a vector of length $m$ with edge weight variables  
>* $r$ is a variable for the regularization degree. 

Notice that this notion of incidence matrix for an arbitrary graph is different from the incidence matrix for a *bipartite* graph. 

Let's code a function that creates the incidence matrix of a graph.

```{r}
library(igraph)

# create incidence matrix for a graph
make_incidence = function(g) {
  n = vcount(g)
  m = ecount(g)
  # get edges as a matrix
  E = get.edges(g, E(g))
  B = matrix(0, nrow = n, ncol = m)
  # build incidence matrix
  for (i in 1:m) {
    B[E[i,1], i] = 1
    B[E[i,2], i] = 1
  }  
  return(B)
}

g = make_star(4, mode = "undirected")
plot(g)
make_incidence(g)


```

## Play

Consider a star graph with center 1 and periphery 2 and 3: 2--1--3. 

Such a graph is not regularizable and therefore the power equation $x = x^{\div} A$ has no solution. Such an equation corresponds to a nonlinear system with 3 equations, one for each node of the graph. 

Try solving the system manually and understand why there is no finite solution.

## Solution

The nonlinear system is:

$$
\begin{array}{lcl}
x_1 & = & \frac{1}{x_2} + \frac{1}{x_3} = 2 x_1 \\
x_2 & = & \frac{1}{x_1} \\
x_3 & = & \frac{1}{x_1}
\end{array}
$$

which has no finite solution. Indeed, the only solutions are $x = (0, \infty, \infty)$ and $x = (\infty, 0, 0)$.


## Perturbation

If matrix $A$ is not regularizable, it can be fully perturbed as follows: 

$$A^{F}_{\epsilon} = A + \epsilon E.$$ 

This corresponds to making the graph complete. A less intrusive diagonal perturbation is sufficient: 

$$A^{D}_{\epsilon} = A + \epsilon I$$ 

where $I$ is the diagonal matrix. 

This corresponds to addiing a loop with weight $\epsilon$ to each node in the graph. 


## Play

Prove that the graph associated to the perturbation matrices $A^{F}_{\epsilon} = A + \epsilon E$ and 
$A^{D}_{\epsilon} = A + \epsilon I$ are regularizable.

## Solution

The graph of $A^{F}_{\epsilon}$ is complete, hence regular, hence regularizable.

The graph of $A^{D}_{\epsilon}$ has one self-loop for every node. This can be used to adjust the weighting of the edges to obtain a solution as follows.

Let $k_i$ be the unweighted degree of node $i$ and let $K$ be the maximum unweighted degree of a node. A regularization solution is the following. Set a weight of 1 to all edges different from a self-loop and let $l_i$ be the weight of the self-loop attached to node $i$. Set $l_i = K - k_i + 1 > 0$. Then it is easy to see that all nodes have weighted degree equal to $K$. Indeed, node $i$ has degree: $$k_i - 1 + l_i = k_i - 1 + K - k_i + 1 = K$$


## Code

The following user-defined function `regularify` checks if a graph is regularizable and in case gives the regularization solution. 

It uses the <a href="https://cran.r-project.org/web/packages/lpSolve/index.html">lpSolve</a> and <a href="https://cran.r-project.org/web/packages/lpSolveAPI/index.html">lpSolveAPI</a> packages for solving linear, integer and mixed integer programs. 

In order to find a positive solution $w > 0$, the program solves the linear problem $$\hat{B} (\hat{w}, r) = -d$$ with $\hat{w} \geq 0$,  where:

>* $\hat{B}$ is the incidence matrix $B$ with one additional column $-e$, 
>* $d$ is a vector with node degrees, and 
>* $(\hat{w}, r)$ is a vector with weight variables $\hat{w}$ and an additional variable $r$ for the regularization degree. 

Setting $w = \hat{w} + 1$, the problem corresponds to $B w = r e$, with the constrain that $w \geq 1$. 


```{r}
library(lpSolve)
library(lpSolveAPI)

regularify = function (g) {
  n = vcount(g)
  m = ecount(g)
  B = make_incidence(g)
  # objective function
  obj = rep(0, m + 1)
  # constraint matrix
  con = cbind(B, rep(-1, n))
  # direction of constraints
  dir = rep("=", n)
  # right hand side terms
  rhs = -degree(g)
  # solve the LP problem
  sol = lp("max", obj, con, dir, rhs)
  # get solution
  if (sol$status == 0) {
    s = sol$solution
    # weights
    w = s[1:m] + 1
    # weighted degree
    d = s[m+1]
  }
  # return the solution
  if (sol$status == 0) {
    return(list(weights = w, degree = d)) 
  }
  else {
    return(NULL)   
  }
}
```



The following user-defined function `power` computes power using a direct computation on a regularizable adjacency matrix A:

```{r}
# Compute power x = (1/x) A 
#INPUT
# A = graph adjacency matrix
# t = precision
# OUTPUT
# A list with:
# vector = power vector
# iter = number of iterations

power = function(A, t) {
  n = dim(A)[1];
  # x_2k
  x0 = rep(0, n);
  # x_2k+1
  x1 = rep(1, n);
  # x_2k+2
  x2 = rep(1, n);
  diff = 1
  eps = 1/10^t;
  iter = 0;
  while (diff > eps) {
    x0 = x1;
    x1 = x2;
    x2 = (1/x2) %*% A;
    diff = sum(abs(x2 - x0));
    iter = iter + 1;
  } 
  # it holds now: alpha x2 = (1/x2) A
  alpha = ((1/x2) %*% A[,1]) / x2[1];
  # hence sqrt(alpha) * x2 = (1/(sqrt(alpha) * x2)) A
  x2 = sqrt(alpha) %*% x2;
  return(list(vector = as.vector(x2), iter = iter))
}  

```

## Play

1. make a star graph and check that it is not regularizable
2. perturb the graph matrix and compute power of the perturbed matrix

## Solution

```{r}
g = make_star(4, mode = "undirected")
plot(g)
```


```{r}
# The graph is not regularizable
regularify(g)

# Use diagonal perturbation
A = as_adjacency_matrix(g)
I = diag(0.15, vcount(g))
(AI = A + I)

# compute power
power(AI, 6)

```


## Dig deeper

> * [A Theory on Power in Networks](../biblio/power.pdf)
> * Bargaining and power in networks. Chapter 12 in book [Networks, crowds and markets](../biblio/networks-book.pdf)



















