Estimating parameters of a directed weighted graph model with beta-distributed edge-weights
Marianna Bolla Ahmed Elbanna Jozsef Mala
math.ST stat.TH 62F10, 62B05



 Deep Learning Models for Wireless Signal Classification with Distributed Low-Cost Spectrum Sensors

    Sreeraj Rajendran, Student Member, IEEE,
        Wannes Meert, Member, IEEE
        Domenico Giustiniano, Senior Member, IEEE,
        Vincent Lenders, Member, IEEE
        and Sofie Pollin, Senior Member, IEEE.
        
        






    December 30, 2023
===============================================================================================================================================================================================================================================





§ ABSTRACT


We introduce a directed, weighted random graph model, where the edge-weights
are independent and beta-distributed with parameters depending on their
endpoints. We will show that the row- and column-sums of the transformed
edge-weight matrix are sufficient statistics for the parameters, and
use the theory of exponential families to prove that the ML estimate of
the parameters exists and is unique. Then an algorithm to find this
estimate is introduced together with convergence proof that uses properties
of the digamma function. Simulation results and applications are also presented.

Keywords: exponential family,
sufficient statistics, ML estimation, digamma function, successive
approximation

MSC2010: 62F10, 62B05. 




§ INTRODUCTION


The theory of ML estimation in the following types of
exponential family random graph models
has frequently been investigated in the last decade,
see, e.g., <cit.>.
The graph has n vertices, and the adjacency relations between them are
given by the n× n random edge-weight matrix =(w_ij )
of zero diagonal.
If  is symmetric, then we have an  undirected
graph; otherwise, our graph is directed, where w_ij
is the nonnegative weight assigned to the i→ j edge according to the
model.


We assume that the edge-weights (above or out of the main diagonal) are
completely independent (but their distribution usually depends on
different parameters),
and have an exponential family distribution _.
So the likelihood function has the general form

    L_ ( )
    
    = e^< ,  ( ) >   -   Z ( ) · h( ) ,

with the canonical parameter ,
log-partition (cumulant) function Z ( ), and
canonical sufficient statistic .
In these random graph models, components of = ( ) are the row-sums
and/or column-sums of  or some -related matrix, i.e., they are
vertex-degrees or in- and out-degrees of the observed undirected or directed,
weighted or unweighted graph (in the weighted case, the edge-weights
may undergo a suitable transformation).
Also, h ( ) is usually 1 over the support of the
likelihood function, indicating that given the canonical sufficient statistics,

the joint distribution of the entries
is uniform (microcanonical) in these models.

To make inferences on the parameters, typically we have only one observation
for the graph.
It may seem that it is a one-element sample, but there are the adjacencies
that form the sample; the number of them is n 2
in the undirected, and n(n-1) in the directed case. The number of parameters,
contained in , is n in the undirected and 2n in the directed case.
The parameters
can be considered as affinities or potentials of the vertices to make
ties in the undirected, and to emanate or adsorb edges in the
directed case. It is important that we divide the components of the
canonical parameter  of the underlying distribution
of the ij or i→ j edge between the connected vertices,
like α_i +α_j in the undirected and
α_i +β_j in the directed case (i j),
see <cit.>.











In regular exponential families (Θ is open), the ML equation
∇_ln L_ ( ) =
is equivalent to

    ∇_ Z( ) =   .

Since ∇_ Z( ) = _,
the ML equation (<ref>) means that the canonical sufficient statistic
is made equal
to its expectation. But when is it possible? Now we briefly summarize
existing theoretical  results on this issue.














Let M = {_  :  ∈Θ} denote the
so-called mean parameter space in the model; it is necessarily convex.
Let  M^0 denote its interior.
When the canonical statistic is also complete,
and hence, minimal sufficient, the representation (<ref>) is
minimal (i.e., the model is not overparametrized).


In exponential family, the gradient mapping ∇ Z:  Θ→ M
is one-to-one if and only if the exponential family representation is
minimal.


In a minimal exponential family, the gradient mapping ∇ Z is onto
M^0.

By Propositions <ref> and <ref>, any parameter in M^0
is uniquely realized by the _θ distribution
for some θ∈Θ. Also, in a regular and minimal
exponential family, M is an open set and  is identical
to M^0.

As the ML estimate of   is the solution of (<ref>),
we have the following.


Assume, the (canonical) parameter space Θ is open. Then there exists a
solution ∈Θ to the ML equation
∇_ Z( ) = if and only if ∈ M^0;
further, if such a solution exists, it is also unique.


Note that in regular and minimal exponential families, M^0 is also the
interior of T, which is the convex hull of all possible values of
, see, e.g., <cit.>.

In the case of discrete distributions, it frequently happens that the boundary
of T has positive measure. For instance, the so-called threshold graphs
are located on the boundary of the polyhedron, determined by the
Erdős–Gallai conditions, in the model
of <cit.> which uses Bernoulli distributed entries.
However, in the case of an
absolutely continuous _ distribution,







the boundary of T
has zero Lebesgue measure, and so, probability zero with respect to the
 _ measure.
Therefore, in view of Proposition <ref>,
the ML equation has a unique solution with probability 1.


The organization of the paper is as follows. In Section <ref>, we
introduce a model for directed edge-weighted graphs and prove that
a unique ML estimate of the parameters exists.
In Section <ref>, we define an iterative algorithm to find
this solution, and prove its convergence with a convenient starting.
In Section <ref>,
the algorithm is applied to randomly generated and real-word data.
In Appendix A, properties of the digamma function, whereas in Appendix B,
the boundary of our
M is discussed.  The long proof of the main
convergence theorem of the iteration algorithm, introduced in
Section <ref>, is presented in Appendix C.

We remark that edge-weighted graphs of uniformly bounded edge-weights are
prototypes of real-world networks, see e.g., <cit.>.
Without loss of generality, if the
edge-weights are transformed into the [0,1] interval, the
beta-distribution
for them, with varying parameters, is capable to model a wide range of
possible probability densities on them. This indicates the soundness of
the model to be introduced in Section <ref>.




§ A  RANDOM GRAPH MODEL WITH BETA-DISTRIBUTED
EDGE-WEIGHTS


Let =(w_ij) be the n× n (usually not symmetric)
edge-weight matrix of a random directed graph
on n vertices: w_ii =0 (i=1,… ,n) and w_ij∈ [0,1] is the
weight of the i→ j edge (i j). Our model is the following:
the i j weight obeys a beta-distribution  with parameters a_i >0 and
b_j >0. The parameters are collected in =(a_1 ,… a_n ) and
=̱ (b_1 ,… b_n ), or briefly, in = ( ,)̱.
Here a_i can be thought of as the potential of the vertex i to send
messages out, and b_i is its resistance to receive messages in.

The likelihood function is factorized as

    L_ , ( )    =
     ∏_i jΓ (a_i  +b_j )/Γ (a_i ) Γ (b_j )
     w_ij^a_i -1 (1-w_ij)^b_j -1
        = C( ,)̱∏_i jexp
     [ (a_i -1 ) ln w_ij + (b_j -1) ln (1-w_ij) ] 
        = exp
     [ ∑_i=1^n (a_i -1 ) ∑_j iln w_ij +
       ∑_j=1^n (b_j -1 ) ∑_i jln (1-w_ij ) -Z ( , )̱],

where C( ,)̱ is the normalizing constant, and
Z ( , )̱ = -ln C( ,)̱ is the log-partition (cumulant) function.
Since the likelihood function depends on  only through the row-sums of
the n× n matrix = ( ) of general entry ln w_ij and the
column-sums  of the n× n matrix =  ()
of general entry ln (1-w_ij ),
by the Neyman–Fisher factorization theorem, the row-sums
R_1 ,… ,R_n of  and column-sums  C_1 ,… ,C_n of  are
sufficient statistics for the parameters. Moreover,
=( , ) = (R_1 ,… ,R_n , C_1 ,… ,C_n) is the canonical
sufficient statistic, which is also minimal.
Note that  contains the log-weights of the original graph,
while  contains the
the log-weights of the complement graph of edge-weight matrix 
with entries 1-w_ij (i j ).

The first factor in the Neyman–Fisher factorization
(including gamma-functions)
depends only on the parameters and on the sample through these
sufficient statistics,
whereas the seemingly not present other
factor – which would merely depend on  – is constantly
1, indicating that the conditional
joint distribution of the entries, given the row- and column-sums of the
log-weight and log-complement matrix
is uniform (microcanonical) in this model.
So under the conditions on the margins of  and , the
directed graphs coming from the above model are uniformly distributed.

The system of likelihood equations is obtained by making the derivatives of
L_ , ( ) with respect to the parameters equal to 0:

    ∂  L_ , ( )/∂ a_i    =
    
    
     ∑_j iψ (a_i  +b_j ) - (n-1) ψ (a_i) +R_i =0,
        i=1,… ,n; 
    ∂ L_ , ( )/∂ b_j    =
    
    
     ∑_i jψ (a_i  +b_j ) - (n-1) ψ (b_j) +C_j =0 ,
       j=1,… ,n .

Here ψ (x)=∂lnΓ (x)/∂ x =
Γ' (x)/Γ (x) for x>0
is the digamma function. For its properties, see Appendix A.



To apply the theory of Section <ref>,
we utilize that the parameter space
Θ⊂_+^2n is open, akin to the canonical
parameter space, (-1 ,∞ )^2n.
Note that the canonical parameter  is, in fact,
(' ,'̱ )=' = -,
where ∈^2n is the vector of all 1 coordinates.
With it, the log-partition function  is

    Z (', ' )=
     -∑_j iΓ ( a_i'  + b_j' +2 ) +
      ∑_j iΓ (a_i' +1)
      +∑_i jΓ (b_j' +1) .

In view of (<ref>), the ML equation is equivalent to

    ∂ Z (', ' )/∂a_i'    =
     -∑_j iψ (a_i'  +b_j' +2 ) +
      (n-1) ψ (a_i' +1)  =R_i ,    i=1,… ,n; 
    ∂ Z (', ' )/∂b_j'    =
     -∑_i jψ (a_i'  +b_j ' +2 ) +
      (n-1) ψ (b_j' +1)  =C_j  ,    i=1,… ,n  .

But this system of equations is the same as (<ref>), in terms of
the parameter ' instead of .

In view of  Section <ref>,
the mean parameter space M consists of parameters
(A_1 ,… A_n ,B_1 ,… ,B_n) obtained by the gradient mapping, that is,

    A_i = A_i ( ,)̱ =- ∑_j i[ ψ (a_i  +b_j ) -
      ψ (a_i  ) ] ,    i=1,… ,n; 
    
     B_j = B_j ( ,)̱ =- ∑_i j[ ψ (a_i  +b_j  ) -
       ψ (b_j ) ]  ,    j=1,… ,n .

M is an open set, whose boundary
is determined by the limit properties between the digamma and the log functions,
see Appendix B for details. There we also find a correspondence between
the points on the boundary of M and those on the boundary
of the  convex hull T
of the possible sufficient statistics =( , ) within _-^2n.
It is interesting that while the boundary points of M do not belong
to the open set M, the boundary points of T do belong to T,
and can be realized as row- and column-sums of the ( ) and ( )
matrices with a  of off-diagonal entries in  (0,1).
However, this boundary has 0 probability, and so, any canonical
sufficient statistic  of the observed graph  is in M,
with probability 1.
Therefore, by Proposition <ref>,  we can state the following.

The system of the ML equations (<ref>) has a unique
solution =(, ̱̂), with probability 1.



Later we will use the following
trivial upper bound for the sum of row- and column-sums (of the  and 
 matrices):

    ∑_i=1^n R_i + ∑_j=1^n C_j =
     ∑_i=1^n ∑_j iln w_ij
     + ∑_j=1^n ∑_i jln (1-w_ij) =
      ∑_i jln [w_ij (1-w_ij )] ≤ -2ln 2   n (n-1)

due to rearranging the terms and the relation w_ij (1-w_ij )≤ 1/4 for
w_ij∈ [0,1] with equality if and only if w_ij=1/2 (i j).
For finer estimates see Appendix B.

Also note that the Hessian of the system of ML equations (consisting of
the second order partial derivatives of L_) at
) does not contain the  sufficient statistics any more, therefore
the negative of it is the Fisher-information matrix at .
Because of
the regularity conditions, the information matrix is positive, and so,
the Hessian is negative definite.


This is also an indication of the existence of a unique ML estimate.




§ ITERATION ALGORITHM TO FIND THE PARAMETERS


To use a fixed point iteration,
now we rewrite the system of likelihood equations in the form
=f( ), where = ( ,)̱, as follows:

    a_i    = ψ^-1[ 1/n-1 R_i +1/n-1∑_j iψ  (a_i  +b_j ) ]
      =:g_i ( ,)̱,    i=1,… ,n 
    
     b_j    = ψ^-1[ 1/n-1 C_j +1/n-1∑_i jψ  (a_i  +b_j ) ]
      =:h_j ( ,)̱,    j=1,… ,n .

Here g_i's and h_j's are the coordinate functions of
f =(g,h) :^2n→^2n.
Then, starting at ^(0), we use the successive approximation
^(it ) :=f (^(it -1) for it =1,2,…, until convergence.
Now the the statement of convergence of the above iteration to the theoretically
guaranteed unique  (see Theorem <ref>) follows.



Let =(,̱̂ ) be the unique solution of the
ML equation (<ref>). Then the above
mapping f=(g,h) is a contraction in some closed
neighborhood K of , and so, starting at any
^(0)∈ K, the fixed point of the iteration
^(it ) =f (^(it -1)) exists and is .

The prof of this theorem is to be found in Appendix C.

Since K is only theoretically guaranteed, we need some practical
considerations about the choice of ^(0 ), which should be
adapted to the sufficient statistics.
In the sequel, for two vectors =(x_1,…,x_n),
=(y_1,…,y_n)  we use the notation
> if x_i>y_i for each i=1,…,n. Likewise,
≥ is the shorthand for x_i≥ y_i for each i=1,…,n.

Recall that f=(g,h) is the mapping (<ref>) of the fixed point
iteration, and
=(,̱̂)>0 is the (only) solution of the
equation
f()=, where ∈^2n is the vector of all 0 coordinates.



Let

    M:=max{max_i∈{ 1,…,n }(-R_i/n-1),
    max_i∈{ 1,…,n }(-C_i/n-1)}

and ε>0 be the (only) solution of the equation
ψ(2x)-ψ(x)=M. Then
(,̱̂)≥ε.


Proof.
In view of  (<ref>) we have that M≥ln 2. Since
equality in (<ref>) is attained with probability 0, we have that
M>ln 2 with probability 1. Therefore, by
Lemma <ref> of Appendix A,
there exists an ε, with probability 1, such that
ψ(2ε )-ψ(ε )=M.

Without loss of generality we can assume that

    â_i_0=min{min_i∈{ 1,…,n }â_i,
     min_i∈{ 1,…,n }b̂_i}.

Then by the ML equation, the monotonicity of ψ, and
Lemma <ref> of Appendix A, we get

    (n-1)M≥ -R_i_0=∑_j≠ i_0ψ(â_i_0+b̂_j)-(
    n-1)ψ(â_i_0)≥
    (n-1)[ψ(2â_i_0)-ψ(â_i_0)] .

Therefore, â_i_0≥ε,
whence â_i,b̂_i≥ε holds
for every i=1,…,n. 
□

 With the solution ε
of ψ(2x)-ψ(x)=M of (<ref>),
we have f(ε )≥ε.


Proof.

    g_i(ε )
    =ψ^-1(ψ(2ε)+R_i/n-1)≥ψ^-1(ψ(2ε)-M)=ε .

 Likewise,

    h_i(ε )
    =ψ^-1(ψ(2ε)+C_i/n-1)≥ψ^-1(ψ(2ε)-M)=ε .    □

It is also clear that we have the following.

If (,)̱≥ (, )>, then f(,)̱≥ f(, ).



With ε satisfying ψ(2ε )-ψ(ε )=M
 of (<ref>), and  starting at
^(0) =ε, the sequence ^(it) of the iteration
^(it)=f(^(it-1)) for it→∞ converges at
a geometric rate 
to the unique solution ( ,̱̂ ) of the ML equation.


Proof.
From Propositions <ref> and <ref>
we obtain that the sequence ^(it) is coordinate-wise increasing.
Moreover,
it is clear that (^(it) ) is bounded from above by
(,̱̂),
due to Proposition <ref>.
Therefore, the convergence of  ^(it) follows, and  by the
continuity of f, the limit is clearly a fixed point of f.
However, in view of Section <ref>, the solution of the ML equation
is a fixed point of f, and it cannot be else but the unique solution
(,̱̂), guaranteed by Theorem <ref>.
Further, from Theorem <ref> we get that the rate of convergence is (at least)
geometric. 
□

Therefore, a good starting can be chosen by these considerations. Also note
that at the above ^(0) and possibly at its first (finitely many)
iterates,
f is usually not a contraction. It becomes a contraction only when
some iterate ^(it_0 ) gets into the neighborhood K of

of Theorem <ref>, which is inevitable in view of the convergence of
the sequence ^(it). So, Theorem <ref> is literally applicable
only if we start the iteration at  ^(it_0 ). In practice,
however, we do not know the theoretically guaranteed neighborhood K.
The practical merit of Theorem <ref> is just that it offers a
realizable starting.



§ APPLICATIONS


First we generated a random directed edge-weighted graph on n=100 vertices.
The edge-weight matrix  had zero diagonal, and the off-diagonal entries
w_ij's were independent. Further, for i j,
the weight w_ij was generated according to
beta-distribution  with parameters a_i >0 and b_j >0, where a_i's and
b_j's were chosen randomly in the interval [1,5].

Then we estimated the parameters based on , and plotted the
a_i , â_i  (i=1,… ,n)  and
b_j , b̂_j  (j=1,… ,n) pairs.

Figure <ref>  shows a good fit between them.











We also applied the
algorithm to migration data between 34 countries. Here w_ij is proportional
to the number of people
in thousands who moved from country i to country j (to find jobs)
during the year 2011,
and it is normalized so that be in the interval (0,1).

The estimated parameters are in Table <ref>.

In this context, a_i's are related to the  emigration and
and b_i's to the counter-immigration potentials.
When a_i is large, country i has a relatively large potential
for emigration. On the contrary, when
b_i is large, country i tends to have a relatively large
resistance against immigration.



It should be noted again that edge-weighted graphs of this type very
frequently model real-world directed networks.





§ APPENDIX




 §.§ A.  Properties of the digamma function


Though, we do not use it explicitly, the following approximation of the
digamma function  ψ (x)=∂lnΓ (x)/∂ x =
Γ' (x)/Γ (x) (x>0)
is interesting for its own right.


ψ (x) = ln (x-1/2 ) +  O( 1/x^2 ) for x>1.

The statement of the lemma easily follows by Taylor expansion.


1/ψ'(x+y)>1/ψ'(x)+1/ψ'(y) for x,y>0.


Proof.
First we prove that the function
u(x)=1/ψ'(x), x∈ (0,∞) is strictly
convex. Indeed,  one can easily see that

    u”(x)=-ψ”'(x)[ψ'(x)]^2+2[ψ”(x)]^2ψ'(x)/[ψ'(x)]^4 ,

and this is positive due to ψ' (x)>0 and the fact that

    [ψ”(x)]^2/ψ”'(x)ψ'(x)>1/2 .

Latter one is a particular case
of Corollary 2.3 in <cit.>.

Now, in view of lim_x→0ψ'(x)=∞, we can extend u continuously
to 0 by setting
u(0)=0. Then u is still strictly convex,
and therefore, for every x,y>0 we have
u(x)=u(y/x+y· 0+x/x+y·(x+y))<
y/x+yu(0)+x/x+yu(x+y).
Consequently,

    u(x)< x/x+y u(x+y),

and likewise,

    u(y)< y/x+yu(x+y).

Adding (<ref>) and (<ref>) together,
we get the statement of the lemma. □


The function ψ(2x)-ψ(x), x∈ (0,∞)
is decreasing and its range is (ln 2,∞).


Proof.
It is easily seen by the identity
ψ(2x)=1/2ψ(x)+1/2ψ(x+1/2)+ln 2
which can be found in <cit.>. □

In the last lemma we collect some limiting properties of the digamma function
and its derivative, see, e.g., <cit.> for details.


The digamma function ψ is a strictly concave, smooth function on
(0,∞) that satisfies the following limit relations:

    lim_x→ 0+ψ (x) =-∞,   lim_x→∞ψ (x) =∞,   lim_x→∞ψ' (x) =0, 
       lim_x→∞ (ψ (x) -ln x ) =0,     lim_x→ 0+ (ψ (2x) -ψ (x)) =∞ .






 §.§ B. Considerations on the boundary of the mean parameter space


In Section <ref>, we saw that the mean parameter space M
consists of 2n-tuples
(A_1 ,… A_n ,B_1 ,… ,B_n) obtained from the parameters
( ,)̱ = (a_1 ,… a_n ,b_1 ,… ,b_n) of the underlying
beta-distributions by Equations (<ref>).

Denoting by L ( ,)̱ =(A_1 (, )̱ , … ,A_n ( ,)̱,
B_1 ( ,)̱,… , B_n ( ,)̱) this dependence, i.e., the
Θ→ M (one-to-one) mapping, a boundary
point L̅ =(A̅_1 ,… ,A̅_n ,B̅_1 ,… ,B̅_n)
of M can be obtained as
L̅ =lim_k→∞ L (^k ,^̱k ), where
^k = (a^k_1 ,… ,a^k_n ), ^̱k = (b^k_1 ,… ,b^k_n ), and

    lim_k→∞ a^k_i   = a̅_i ∈ [0,∞],    i=1,… n; 
    lim_k→∞ b^k_j   = b̅_j ∈ [0,∞],    j=1,… n .

In view of Lemma <ref>, only the
a̅_i,b̅_j =∞ cases have relevance. 
The sequence (^k ,^̱k ) can be chosen such that

    lim_k→∞a^k_i/b^k_j =x_ij  with  
     0<x_ij < ∞ ,   for    i j.

Then, using <ref>),

    A̅_i    = -lim_k→∞∑_j i[ ψ (a^k_i+
      b^k_j)- ψ (a^k_i ) ]
     =lim_k→∞∑_j i[ ln (a^k_i ) -
      ln (a^k_i  +b^k_j) ] =
     ∑_j ilnx_ij/1+x_ij,     i=1,… ,n; 
    B̅_j    = -lim_k→∞∑_i j[ ψ (a^k_i+
       b^k_j) - ψ (b^k_j ) ]
     =lim_k→∞∑_i j[ ln (b^k_j ) -
      ln (a^k_i  +b^k_j) ] = ∑_i jln1/1+x_ij ,
        j=1,… ,n .

These equations show that the boundary point L̅ of M
contains – in its coordinates –
the row- and column-sums of the matrices ( ) ad ( ) respectively
(see Section <ref>), where the general off-diagonal entry of the
 n× n edge-weight matrix  is x_ij/1+x_ij.

Observe that 2n-1 x_ij's can be chosen free, and all the others
are obtainable from them. To see this, consider the complete bipartite graph
on vertex classes (a_1 ,… ,a_n) and (b_1 ,… ,b_n),
where to the edge connecting a_i and b_j we assign a_i /b_j.
Choose a minimal spanning tree of this graph (it contains 2n-1 edges), and
consider the sequence of (^k ,^̱k )'s satisfying condition (<ref>).
Then, as k→∞,
the  x_ij's of the edges not included in the spanning tree
can be obtained from the x_ij's  of the 2n-1
edges included in  the spanning tree. Therefore the row- and
column-sums of the  edge-weight matrix  of entries
x_ij/1+x_ij (i j)
are on a (2n-1)-dimensional manifold
in _-^2n, so they are on the boundary of the convex hull T
of the possible sufficient statistics ( , ).
However,
this boundary has zero Lebesgue measure, and so, zero probability
with respect to the underlying absolutely continuous distribution.





 §.§ C.  Proof of Theorem <ref>


It suffices to prove that some induced matrix norm of the matrix of the
first  derivatives  of f at  is strictly less than 1.
We prove this for the L_1-norm. From (<ref>) we obtain that

    ∂ g_i/∂ a_i (, b̂ ) =
    1/n-1∑_j≠
    iψ'(â_i+b̂_j)/ψ'[ ψ^-1(1/n-1∑_j≠ iψ(â_i+b̂_j)+R_i/n-1)] .

From (<ref>) we have 1/n-1∑_j≠ iψ(â_i+b̂_j)+R_i/n-1=ψ(â_i).
Substituting it into (<ref>), we get

    ∂ g_i/∂ a_j(,b̂)= {[ ∑_s≠ i1/n-1ψ'(â_i+b̂_s)/ψ'(â_i)                             j=i;                               0                          j≠ i . ].

Likewise,

    ∂ g_i/∂ b_j(,b̂)= {[                         0                       j=i; 1/n-1ψ'(â_i+b̂_j)/ψ'(â_i)                    j≠ i . ].

Further,

    ∂ h_i/∂ a_j(,b̂)= {[                          0                        j=i; 1/n-1ψ'(â_j+b̂_i)/ψ'(b̂_i)                       j≠ i ].

and

    ∂ h_i/∂ b_j ( ,̱̂ )= {[ ∑_s≠ i1/n-1ψ'(â_s+b̂_i)/ψ'(b̂_i)                              j=i;                                0                           j≠ i . ].

Observe that J ( ,b̂ ) has nonnegative entries. Therefore,
its L_1-norm is the maximum of its column-sums.
The jth column-sum of J ( ,b̂ ) is equal to

    ∑_s≠ j1/n-1ψ'(â_j+b̂_s)/ψ'(â_j)
    + ∑_s≠ j1/n-1ψ'(â_j+b̂_s)/ψ'(b̂_s)
    =1/n-1∑_s≠ jψ'(â_j+b̂_s)(
    1/ψ'(â_j)+1/ψ'(b̂_s))

for j=1,… ,n; and likewise, the (n+j)th column-sum of
J ( ,b̂ ) is

    1/n-1∑_s≠ jψ'(â_s+b̂_j)(
    1/ψ'(â_s)+1/ψ'(b̂_j))

for j=1,… ,n.

As (<ref>) and (<ref>) are of similar appearance, it suffices
to prove that the right hand side of  (<ref>) is less than 1.
But ψ'(â_j+b̂_s)(
1/ψ'(â_j)+1/ψ'(b̂_s))<1 holds by Lemma <ref>,
and we have n-1 terms in the summation.

Since f: ^2n→^2n is continuously differentiable in a
neighborhood of =( ,̱̂), Theorem 3
of <cit.> implies that there is a closed neighborhood K of
 such that f is a contraction on K.
In particular, the fixed point iteration f (^(it -1) )=
^(it ) (it→∞ )
converges for every ^(0 )∈ K to ,
which is the unique solution of (<ref>). □


1

handbook
Abramowitz, M., Stegun, I. A., eds.,
Handbook of Mathematical Functions with
Formulas, Graphs, and Mathematical Tables (10th ed.).
New York: Dover (1972), pp. 258-259.

alzer
Alzer, H., Wells, J., Inequalities for the polygamma functions,
SIAM J. Math. Anal. 29 (6) (1998), 1459-1466.

Bernardo
Bernardo, J. M., Psi (digamma) function. Algorithm AS 103,
Applied Statistics 25 (1976), 315–317.

Bolla13
Bolla, M., Spectral clustering and biclustering. Wiley (2013).

Bolla
Bolla, M., Elbanna, A., Estimating parameters of a probabilistic heterogeneous
block model via the EM algorithm,
Journal of Probability and Statistics (2015),
Article ID 657965.

Chatterjee
Chatterjee, S., Diaconis, P. and Sly, A.,  Random graphs with a given degree
sequence, Ann. Stat. 21 (2010), 1400–1435.

Grasmair
Grasmair, M., Fixed point iterations, https://wiki.math.ntnu.no/_media/ma2501/2014v/fixedpoint.pdf





Hillar13
C. J. Hillar, A. Wibisono, Maximum entropy distributions on
graphs, arXiv:1301.3321v2 (2013).





Lauritzen
Lauritzen, S. L., Graphical Models. Oxfor Univ. Press (1995).




Yan Yan, T., Leng, C., Zhu, J.,
Asymptotics in directed exponential random graph models with an
increasing bi-degree sequence, Ann. Stat. (2016)
44, 31-57.

Wainwright
M. Wainwright,. M. I. Jordan, Graphical models, exponential families, and
variational inference, Foundations and Trends in Machine Learning
1 (1-2), 1-305 (2008). 




