On the Determinant and its Derivatives of the Rank-one Corrected Generator of a Markov Chain on a Graph
Jerzy A Filar Michael Haythorpe Walter Murray




J.A. Filar
Flinders University

jerzy.filar@flinders.edu.au
M. Haythorpe
Flinders University

michael.haythorpe@flinders.edu.au
W. Murray
Stanford University

walter@stanford.edu

On the Determinant and its Derivatives of the Rank-one Corrected Generator of a Markov Chain on a Graph.
    J.A. Filar, M. Haythorpe, W. Murray
    December 30, 2023
========================================================================================================

We present an algorithm to find the determinant and its
first and second derivatives of a
rank-one corrected generator matrix of a doubly stochastic Markov chain.
The motivation arises from the fact that the global minimiser of this
determinant solves the Hamiltonian cycle problem. It is essential for
algorithms that find global minimisers to evaluate both first and second
derivatives at every iteration. Potentially the computation of these
derivatives could require an overwhelming amount of work since  for the
Hessian N^2 cofactors are required. We show how the  doubly stochastic
structure and the properties of the objective may be exploited to calculate all cofactors from a single LU
decomposition.





§ INTRODUCTION




The Hamiltonian cycle problem (HCP) is an important graph theory problem that features prominently in
complexity theory because it is NP-complete <cit.>.
HCP has also gained recognition because two
special cases: the Knight's tour and the Icosian game, were
solved by Euler and Hamilton,
respectively. Finally, HCP is closely related to the well known Traveling Salesman problem.

The definition of HCP is the following: given a graph Γ containing N nodes,
determine whether any simple cycles of length N exist in the graph.
These simple cycles of length N are known as Hamiltonian cycles.
If Γ contains at least one Hamiltonian cycle (HC), we say that
Γ is a Hamiltonian graph. Otherwise, we say that Γ
is a non-Hamiltonian graph.

While there are many graph theory techniques that have been designed to solve HCP, another common approach is to
associate a variable x_ij with each arc (i,j) ∈Γ, and solve an associated optimisation problem. A convenient method of representing these constraints is to use a matrix P(x), where

    p_ij(x)     =    {[     x_ij,           (i,j) ∈Γ,;        0,                     ].

The discrete nature of HCP naturally lends itself to integer programming optimisation problems. However, arising from an embedding of HCP in a Markov decision process (a practice initiated by Filar et al <cit.>), continuous optimisation problems that are equivalent to HCP have been discovered in recent times. In particular, it was demonstrated in <cit.> that if we define

    A(P(x))     =     I - P(x) + 1/Nee^T,

where e is a column vector with unit entries, then HCP is equivalent to solving the following optimisation problem:

    min  -(A(P(x)))

subject to

    ∑_j ∈𝒜(i)x_ij    =     1,    i = 1, , N,
    ∑_i ∈𝒜(j)x_ij    =     1,    j = 1, , N,
    
    x_ij   ≥    0,   ∀ (i,j) ∈Γ,

where 𝒜(i) is the set of nodes reachable in one step from node i. Constraints (<ref>)–(<ref>) are called the doubly-stochastic constraints. For neatness, we refer to constraints (<ref>)–(<ref>) as the set 𝒟𝒮, and we call the objective function f(P(x)). Then, the above problem can be represented as follows:


    min{f(P(x))    |    x∈𝒟𝒮}.


Note that we need to find a global minimiser and, typically, there are many of them. However, the number of global minimisers is typically extremely small compared to the number of local minimisers. One consequence of multiple global minimisers is that there is a similarly large number of stationary points. To distinguish such stationary points from minimisers, algorithms (see for example <cit.>) require the use of second derivatives. It is not hard to appreciate that evaluating such derivatives will be expensive even for moderately-sized problems unless some special structure is identified. We exploit the structure of the Hessian and the fact the points at which it is evaluated are in 𝒟𝒮. We assume an algorithm to solve this problem starts at a feasible point and all iterates remain feasible. For problems with nonlinear objectives this is almost always the best approach. Finding a point in 𝒟𝒮 is simple and does not involve evaluating the objective or its derivatives.





§ PRELIMINARY RESULTS



Our work is based on some key properties of the LU factorisation of the matrix I - P. We adopt the following common notation. As noted we denote a vector of unit elements by e. In a given proof e will still be used even when different occurrences may be of different dimension since its dimension can be inferred. The vector e_j has all zero elements except the jth, which is unity. We will often denote the (i, j)th elements of a matrix, say W, by w_ij.

A matrix G is said to have property S if

    g_ii   ≥    0,   ∀ i = 1, , N,
    
        g_ij   ≤    0,   ∀ i ≠ j,
    
        Ge = G^Te     =     0.



    If G has property S the LU decomposition of G exists.


Following common convention the diagonal elements of L are 1. There is no loss of generality if we assume that g_11≠ 0. After one step of Gaussian elimination (GE) we obtain

    G = ([   1 0^T;   ℓ   I ])([   1 0^T;   0   Ĝ ])([ u_11  u^T;    0    I ]).

Note that, expanding this matrix multiplication gives

    G = ([     u_11      u^T;    u_11ℓ Ĝ +ℓ u^T ]),

and therefore u^T = [[ g_12 g_13    ⋯ g_1N ]] and ℓ = 1/u_11[[ g_21 g_31    ⋯ g_N1 ]]^T. Clearly, u_11 = g_11≥ 0, and therefore u^T ≤ 0 and ℓ≤ 0.
We have G e =  0, which implies

    ([   1 0^T;   0   Ĝ ])([ u_11  u^T;    0    I ]) e =  0.

It follows that Ĝ e = Ĝ^T e =  0. By definition we have

    ĝ_ij = g_i+1,j+1 - ℓ_iu_j,   ∀ i ≠ j.

Since l_i ≤ 0 and u_j ≤ 0 it follows that ĝ_ij≤ 0, ∀ i ≠ j. From this result and Ĝ e =  0 it follows that ĝ_ii≥ 0 and that Ĝ has property S. We can now proceed with the next step of GE. Note that if ĝ_11 = 0 we must have the first row and column of Ĝ be zero, which implies the corresponding row of U is zero and the off-diagonal elements of L are zero. Consequently, we can proceed until the leading diagonal element of Ĝ is nonzero.


 
        Regardless of the rank of G  we have u_NN = 0.


    After N-1 steps of GE Ĝ is a 1× 1 matrix. The only 1 × 1 matrix that satisfies property S is 0.




If G satisfies property S and G = LU, then Ue = 0.

The result follows immediately from the nonsingularity of L and the property that Ge = 0.

If G satisfies property S and G = LU, the off-diagonal elements of L and U are non-positive. Moreover, -1 ≤ l_ij≤ 0 for i  j. 




The first part follows immediately from the fact that u^T ≤ 0 and l≤ 0 in each iteration of GE. Note that when performing the GE the pivot used is always the largest in magnitude. This is a consequence of Ĝ^T e = 0 and the off-diagonal elements in the column of Ĝ being negative. Since the elements being eliminated are not bigger in magnitude than the pivot it follows that the magnitude of elements of L are not bigger than 1.



The above result has important consequences for the condition number of L.


If G satisfies property S, is  rank N - 1 and G = LU, then L^T e = e_N.


Since U has rank N-1 and u_NN = 0, the basis for the nullspace of U^T is κ e_N, where κ 0. From property S,


    G^T e = U^T (L^T e) = 0


and hence we obtain L^T e = κ e_N, for some κ. Since L has unit diagonal entries, we have κ = 1, which yields the result.
















    When G  satisfies property S and is rank N-1 then G̅≡ G +  e e_N^T is nonsingular. 




We have

    G̅= LU +  ee^T_N = L(U + v e^T_N)= LU̅,

where L v = e. Note that U̅ = U +v e^T_N is upper triangular. Moreover,
the (N,N)th element of U̅ is v_N. Since

    e^T Lv = e_N^T v = v_N = e^T e = N,

it implies that U̅, and hence G̅, is nonsingular.










Clearly G = I-P has property S. It should be noted that the existence of an LU decomposition for G has been demonstrated previously for irreducible P in Heyman <cit.>. In the context of this paper, we note that when P(x) corresponds to the strict interior of 𝒟𝒮, irreducibility applies.






§ COMPUTING (A(P)) USING AN LU DECOMPOSITION




When appropriate we suppress the (fixed) argument x in P(x), A(P(x)) and f(P(x)), and write simply P, A(P) (or just A) and f(P), respectively. The most efficient way to compute the determinant of a matrix is to compute its LU decomposition. Normally we obtain LU = Π A, where Π is a permutation matrix. However, we shall show that det (A) may be computed from either the determinant of  the leading principal minor of  I - P or of det(I - P + ee_N^T). Consequently, no permutation matrix is needed. This is of importance since unlike A we expect I - P to be sparse, which is the case of interest for finding HC.  Factorising sparse matrices is much faster than factorising dense matrices. The difference is even greater when it is unnecessary to do numerical pivoting. Knowing apriori where the fill-in is in the factors will enable more efficient data management and avoid indirect addressing.

We assume an algorithm to solve (<ref>) starts at a feasible point and all iterates remain feasible. For problems with nonlinear objectives this is almost always the best approach. Finding a point in 𝒟𝒮 is simple and does not involve evaluating the objective or its derivatives. It is known that the global minimiser of (<ref>) is -N. Indeed, from Theorem 4.1, Proposition 4.4 and Proposition 4.6 in Ejov et al <cit.>, we know that -(A(P(x))) ∈ [-N,0], where the upper bound is obtained when A(P(x)) is singular. However, we have no interest in finding the determinant of A when it is singular. For the rest of this paper we assume x ∈𝒟𝒮 and that A(P(x)) is nonsingular. Since A(P(x)) is merely a rank-one correction of I - P, and the latter is always singular for x ∈𝒟𝒮, the nonsingularity of A(P(x)) implies that I - P has rank N - 1.





 §.§ Product forms of A(P) and (A(P))

To calculate the objective function f(P), its gradient and Hessian, we begin by performing an LU decomposition
to obtain

    LU     =     G = I - P.

 While there may be a need to reorder the matrix to obtain a sparse factorisation. Since the pattern of nonzero elements in G is symmetric then symmetric pivoting may be used, which is equivalent to renumbering the nodes of the graph. It follows we may assume the effort to compute the LU factors is O(N^2). Note this reordering is done once since it depends only on the location of nonzeros and not their values. Whatever ordering is chosen the LU factors of G reordered exist without the need for pivoting based on numerical considerations.


The outline of the derivation of f(P) in terms of U is as follows.


  (1) We express A(P) as the product of three nonsingular factors.

  (2) We show that two of these factors have a determinant of 1.

  (3) We show that the third factor shares all but one eigenvalue with U, with the single different eigenvalue being N (rather than 0).

  (4) We express the determinant as a
product of the first N-1 diagonal elements of U, and N.
First, we express A as a product of L and another factor. Let v be an N × 1 vector, and U̅ be an N × N matrix, such that

    Lv = eU̅ := U + ve^T_N,

where e_N^T = [[ 0 ⋯ 0 1 ]].
Since L is
nonsingular, v 0 exists and is unique, and therefore U̅ is well-defined. The first N - 1 columns of U̅ are identical to those of U.  Consequently U̅ is also upper triangular. Since LU = G satisfies property S, it follows from Lemma <ref> that u_NN = 0, and therefore from (<ref>) that u̅_NN = v_N = N.
 Therefore,

    U̅ =  ∏_i=1^N u̅_ii =
    N ∏_i=1^N-1 u_ii,


where u_ii is the i-th diagonal element of U. Exploiting (<ref>) we may write

    A     =     (I - P) + 1/Nee^T   =   LU + 1/NLve^T   =   L(U + 1/Nve^T)
        =     L(U̅
    + v[1/Ne^T - e_N^T]).








































































Since U̅ is nonsingular we may define w to be the unique solution to the system

    U̅^T w    =    1/Ne - e_N.

Then, from (<ref>)




    A = L(U̅ + v[1/Ne^T - e_N^T]) = L(U̅ + vw^TU̅) = L(I + vw^T)U̅.

We take the determinant of (<ref>) to obtain

    (A)     =    (L)(I + vw^T)(U̅).

Note that, for any vectors c and d,

    (I +  cd^T)     =     1 + d^Tc.

This is because cd^T has one eigenvalue d^Tc of multiplicity 1 and an eigenvalue 0 of multiplicity N-1. Consequently,

    (I + vw^T)     =     1 + w^Tv,

which we substitute into expression (<ref>) above.


The inner-product w^Tv in (<ref>) satisfies

    w^T v    =     0.


From their respective definitions (<ref>) and (<ref>),

    w^T     =    (1/Ne^T - e_N^T)(U̅)^-1,
    v    =     L^-1e.

Then, from (<ref>)–(<ref>) we obtain

    w^T v    =    (1/Ne^T - e_N^T)(U̅)^-1L^-1e  =  (1/Ne^T - e_N^T)(LU̅)^-1e
        =    (1/Ne^T - e_N^T)(I - P + ee_N^T)^-1e.

Since I, P and ee_N^T are all stochastic matrices, we know that I - P + ee_N^T has row sums of 1 as well. Hence, its inverse also has row sums equal to 1, that is,

    (I - P + ee_N^T)^-1e    =    e.

Substituting (<ref>) into (<ref>), we obtain

    w^T v = (1/Ne^T - e_N^T)e = 0,

which concludes the proof.
We now derive the main theorem of this subsection.

Let LU denote the LU decomposition of I - P and u_ii be the diagonal elements of U then

    (A(P))     =     N ∏_i=1^N-1 u_ii.



From (<ref>), (<ref>) and Lemma <ref> we know that

    (A(P))     =    (L)(1 + 0)(U̅).

From the construction of the LU decomposition we know that (L) = 1 and using (<ref>)
gives

    (A(P))     =    (U̅)= N∏_i=1^N-1
    u_ii.




This concludes the proof.
Note that finding v and w is a simple process because L and U̅^T are lower-triangular matrices, so we can solve the systems of linear equations in (<ref>) and (<ref>) directly. 


 §.§ Finding the gradient g(P)

Next we use the LU decomposition found in Subsection <ref> to find the gradient of
f(P) = -(A(P)). Note that since variables of f(P) are entries x_ij of the probability transition
matrix P(x), we derive an expression for g_ij(P) := ∂ f(P)/∂ x_ij for each
x_ij such that (i,j) ∈Γ.

Consider vectors a_j and b_i satisfying the equations U̅^Ta_j = e_j and Lb_i = e_i, where e_j is a zero vector except for a unit entry in the j-th column. Then, we define Q := I - vw^T, where v and w^T are as in (<ref>)–(<ref>). We prove the following result
in this subsection:

    g_ij(P) = (A(P))(a_j^TQb_i),

where g_ij(P) is the gradient vector element corresponding to the arc (i,j) ∈Γ.

The outline of the derivation of (<ref>) is as follows.

  (1) We represent each element g_ij(P) of the gradient vector as a cofactor of A(P).

  (2) We construct an elementary matrix that transforms matrix A(P) into a matrix with determinant equal to the above cofactor of A(P).

  (3) We then express the element g_ij(P) of the gradient vector as the product of (A(P)) and the determinant of the elementary matrix, the latter of which is shown to be equal to a_j^TQb_i.
For any matrix V = (v_ij)_i,j=1^N,N it is well-known (e.g., see May <cit.>) that
∂(V)/∂ v_ij = (-1)^i+j(V^ij), where V^ij is the (i,j)-th
minor of V. That is, ∂(V)/∂ v_ij is the (i,j)-th cofactor of V. Since
the (i,j)-th entry of A(P) is simply a_ij = δ_ij - x_ij + 1/N (where δ_ij is the Kronecker delta that is 1 if i = j, and 0 otherwise), it now follows that

    g_ij(P) = ∂ f(P)/∂ x_ij = ∂[-A(P)]/∂ a_ijda_ij/dx_ij = (-1)^i+j(A^ij(P)).

However, rather than finding the cofactor we calculate (<ref>) by finding
the determinant of a modification of A in which the ith row  has been replaced with e_j^T. Since A is a full-rank matrix, it is
possible to perform row operations to achieve this. Suppose A is composed of rows r_1^T, r_2^T, , r_N^T. Then, we perform the following row operation:

    r_i^T    →   α_j(1)r_1^T + α_j(2)r_2^T +  +
    α_j(N)r_N^T,

where α_j(i) is the i-th element of vector α_j and A^Tα_j =
e_j. This row operation replaces the i-th row of A with α_j^T A = e_j^T, as desired.

In this case, from (<ref>), A^T = U̅^T(I + wv^T)L^T. Since A is nonsingular
α_j can be found directly:

    α_j     =    (A^T)^-1e_j   =  [U̅^T(I + wv^T)L^T]^-1e_j 
        =     (L^T)^-1(I + wv^T)^-1(U̅^T)^-1e_j.



    (I + wv^T)^-1    =     I - wv^T.

Consider

    (I + wv^T)(I - wv^T)     =     I - wv^T + wv^T - wv^Twv^T
        =     I - wv^Twv^T
        =     Iv^Tw = w^Tv =
    0,

Therefore, (I + wv^T)^-1 = (I - wv^T).
Taking the above result and substituting into (<ref>), we obtain

    α_j    =    (L^T)^-1(I - wv^T)(U̅^T)^-1e_j.

Next, we define an elementary matrix E_ij by

    E_ij    :=     I - e_ie_i^T + e_iα_j^T,

and note that it performs the desired row operation (<ref>) on A because

    E_ijA     =     A - e_ir_i^T + e_ie_j^T,

in effect replacing the i-th row of A with e_j^T. Therefore,

    g_ij(P)     =     (-1)^i+j(A^ij)   =  (E_ijA)   =  (E_ij)(A).

From (<ref>), we rewrite E_ij = I - e_i(e_i - α_j)^T. Then, from (<ref>)
we obtain

    (E_ij)     =     1 - (e_i -
    α_j)^Te_i   =   1 - e_i^Te_i + α_j^Te_i   =  α_j^Te_i.

Substituting (<ref>) into (<ref>) we obtain

    (E_ij)     =    e_j^T(U̅)^-1(I - vw^T)(L)^-1e_i.

For convenience we define Q := I - vw^T. Then

    (E_ij)     =     a_j^TQb_iU̅^Ta_j
    = e_jLb_i = e_i.

We now derive the main result of this subsection.
 The general gradient element of f(P) is given by

    g_ij(P) = ∂ f(P)/∂ x_ij    =    (A(P))(a_j^TQb_i).

Substituting (<ref>) into (<ref>) immediately yields the result.
Note that we can calculate all a_j and b_i in advance, by solving the systems of linear
equations in (<ref>), again in reduced row echelon form. Then, for the sake of efficiency we first calculate

    q̂_j^T     :=     a_j^TQ,   j = 1, , N,

and then calculate

    q̂_ij    :=    q̂_j^Tb_i,   i = 1, , N,   j = 1, , N.

This allows us to rewrite the formula for g_ij(P) as

    g_ij(P) = -f(P)q̂_ij.



 §.§ Finding the Hessian matrix H(P)

Here, we show that the LU decomposition found in Subsection <ref> can also be used to calculate the Hessian
of f(P) efficiently. Consider g_ij and q̂_ij as defined in (<ref>)
and (<ref>) respectively. We prove the following result in this subsection:

    H_[ij],[kℓ](P) := ∂^2 f(P)/∂ x_ij∂ x_kℓ = g_kjq̂_iℓ - g_ijq̂_kℓ,

where H_[ij],[kℓ] is the general element of the Hessian matrix corresponding to arcs (i,j) and
(k,ℓ)  ∈Γ.

The outline of the derivation is as follows.

  (1) We represent each element H_[ij],[kℓ](P) of the Hessian matrix as a cofactor of a minor of A(P).

  (2) We construct a second elementary matrix that in conjunction with E_ij (see (<ref>)) transforms matrix A(P) into one with a determinant equivalent to the (k,ℓ)-th cofactor of A^ij(P).

  (3) We then show that the general element of
the Hessian matrix is the product of (A(P)) and the determinants of the two elementary matrices.

  (4) Using results obtained from finding g(P) in Subsection <ref>, we obtain these values immediately.
We define A^[ij],[kℓ] to be the matrix A with rows i,k and columns j,ℓ removed. An argument similar to that for
g_ij(P) in the previous subsection can be made that finding

    H_[ij],[kℓ](P) = ∂^2
    f(P)/∂ x_ij∂ x_kℓ = (-1)^(i+j+k+ℓ+1)(A^[ij],[kℓ])i ≠ kj ≠ℓ,

is equivalent to finding the negative determinant of A with the ith and kth rows changed to e_j^T and e_ℓ^T respectively. That is,

    ∂^2 f(P)/∂ x_ij∂ x_kℓ    =    
    -(Ê_kℓE_ijA(P))
        =     -(Ê_kℓ)(E_ij)(A(P)),

where Ê_kℓ is an additional row operation constructed to change row k of E_ijA into e_ℓ^T. Note
that if i = k or j = ℓ, the matrix A^[ij],[kℓ] is no longer square and the determinant no longer exists.
If this occurs, we define H_[ij],[kℓ] := 0. If both i = k and j = ℓ, we also define H_[ij],[kℓ] := 0, as the determinant is linear in each element of A(P).

Consider E_ijA composed of rows r̂_1^T, r̂_2^T, , r̂_N^T. Then, we
perform the following row operation:

    r̂_k    →   γ_ℓ(1)r̂_1 + γ_ℓ(2)r̂_2 +  + γ_ℓ(N)r̂_N,

where (E_ijA)^Tγ_ℓ = e_ℓ. Then, similarly to (<ref>), we directly find γ_ℓ:

    γ_ℓ    =     (E_ij^T)^-1(L^T)^-1(I - wv^T)(U^T)^-1e_ℓ.

Next, in a similar fashion to (<ref>), we construct an elementary matrix Ê_kℓ

    Ê_kℓ    =     I - e_ke_k^T +e_kγ_ℓ^T
        =     I - e_k(e_k^T - γ_ℓ^T).

Then, we evaluate (Ê_kℓ):

    (Ê_kℓ)     =     1 - (e_k^T - γ_ℓ^T)e_k
        =     1 - 1 + γ_ℓ^Te_k 
        =    e_ℓ^T(U̅)^-1QL^-1(E_ij)^-1e_k.

Recall from (<ref>) that E_ij = I - e_i(e_i^T - α_j^T). We have

    (E_ij)^-1    =     I + 1/α_j^Te_ie_i(e_i^T - α_j^T).

Recall from (<ref>) that g_ij = (A)(E_ij), and from (<ref>)
that α_j^Te_i = (E_ij) ≠ 0, and therefore (<ref>) holds. Then,

    α_j^Te_i     =    g_ij/(A).

Substituting (<ref>) into (<ref>) we obtain

    (E_ij)^-1    =     I + (A)/g_ije_i(e_i^T - α_j^T),

and further substituting (<ref>) into (<ref>), we obtain

    (Ê_kℓ)     =    e_ℓ^T(U̅)^-1QL^-1(I + (A)/g_ije_i(e_i^T - α_j^T))e_k
        =    e_ℓ^T(U̅)^-1QL^-1(e_k + (A)/g_ije_ie_i^Te_k - (A)/g_ije_i α_j^Te_k).

Note that since i ≠ k, e_ie_i^Te_k = 0, and from (<ref>),
α_j^Te_k = g_kj/(A). Hence, from (<ref>) and
(<ref>) we obtain

    (Ê_kℓ)     =    e_ℓ^T(U̅)^-1QL^-1(e_k
    - g_kj/g_ije_i) 
        =    a_ℓ^TQ(b_k - g_kj/g_ij b_i).

We now derive the main result of this subsection.
The general element of the Hessian of f(P) is given by

    H_[ij],[kℓ]    =     g_kjq̂_iℓ - g_ijq̂_kℓ,

where q̂_iℓ and q̂_kℓ are defined in (<ref>).
From (<ref>) and (<ref>), we can see that H_[ij],[kℓ] = -(Ê_kℓ)g_ij.
Then, from (<ref>), (Ê_kℓ) = a_ℓ^TQ(b_k -
g_kj/g_ijb_i) and
so H_[ij],[kℓ] = - a_ℓ^TQ(b_kg_ij - b_ig_kj).

In order to improve computation time, we take advantage of the fact that we evaluate every q̂_ij
while calculating the gradient to rewrite the second order partial derivatives of f(P) as

    H_[ij],[kℓ]    =     g_kja_ℓ^TQb_i - g_ija_ℓ^TQb_k 
        =     g_kjq̂_iℓ - g_ijq̂_kℓ.

This concludes the proof.
Note that in practice, we do not calculate some g_kj's when calculating g(P) as an arc (k,j)
need not exist in the graph. In these cases we find g_jk using the gradient formula,
g_jk = -f(P)(q̂_jk), which remains valid despite arc (k,j) not appearing in the graph.


 §.§ Leading principal minor

It is, perhaps, interesting that instead of using the objective function f(P) = -(I - P + 1/Nee^T), it is also possible to use f^1(P) := -(G^NN(P)), the negative determinant of the
leading principal minor of I - P. The following, somewhat surprising, result justifies this claim.


  (1)f^1(P) = 1/Nf(P) = -1/N(I - P + 1/Nee^T).

  (2) If the graph is Hamiltonian, then

    min_P ∈𝒟𝒮 f^1(P) = -1.

First, we show part (1), that is, f^1(P) = 1/Nf(P). To find f^1(P), we construct LU = I - P as before, and define L̂, Û as:

    L̂ = [[   e_1^TL;        ⋮; e_N-1^TL;    e_N^T ]],  Û = [[   Ue_1      ⋯ Ue_N-1    e_N ]].

That is, L̂ is
the same as L with the last row replaced by e_N^T, and Û is the same as U with the last
column replaced with e_N. Then consider

    L̂Û    =    [[   e_1^TL;        ⋮; e_N-1^TL;    e_N^T ]][[   Ue_1      ⋯ Ue_N-1    e_N ]]
        =    [[     e_1^TLUe_1              ⋱   e_1^TLUe_N-1      e_1^TLe_N;              ⋮              ⋱              ⋮              ⋮;   e_N-1^TLUe_1              ⋯ e_N-1^TLUe_N-1    e_N-1^TLe_N;      e_N^TUe_1              ⋯    e_N^TUe_N-1       e_N^Te_N ]].

Since L is lower-triangular, e_i^TLe_N = 0 for all i ≠ N. Likewise, since U is
upper-triangular, e_N^TUe_j = 0 for all j ≠ N. Therefore the above matrix simplifies to

    L̂Û    =    [[     e_1^TLUe_1              ⋯   e_1^TLUe_N-1              0;              ⋮              ⋱              ⋮              ⋮;   e_N-1^TLUe_1              ⋯ e_N-1^TLUe_N-1              0;              0              ⋯              0              1 ]],

which is the same as LU with the bottom row and rightmost column removed, and a 1 placed in the bottom-right element.
Therefore, (L̂Û) = (G^NN(P)), and consequently

    f^1(P)     =    
    -(L̂)(Û).

Note that L̂ and Û are triangular matrices, so


    (L̂) = ∏_i = 1^N
    l̂_ii(Û) = ∏_i = 1^N û_ii.

However, only the last diagonal elements of L̂ and Û are different from L and U̅ (see (<ref>)) respectively, so

    (L̂) = l̂_NN∏_i = 1^N-1
    l_ii(Û) = û_NN∏_i = 1^N-1u̅_ii.

Now, since l̂_NN = l_NN = 1, we have

    (L̂) = (L) =
    1.

We also have û_NN = 1, but by Lemma <ref>, u̅_NN = N and hence

    (Û) = 1/N(U̅).

Therefore, substituting (<ref>) and (<ref>) into (<ref>) we obtain

    f^1(P)     =    
    -(L̂)(Û)
        =     -1/N(U̅)
        =    
    -1/N(I - P + 1/Nee^T) = 1/Nf(P).

Therefore, part (1) is proved.

The proof of part (2) of Theorem <ref> follows directly from the fact that
min f(P) = -N (proved in <cit.>),
and part (1).

Using the leading principal minor has the advantage that the rank-one modification
1/Nee^T is not required, which makes calculating the gradient and the Hessian even simpler than described in Subsection <ref> and Subsection
<ref> respectively. The derivation of the gradient and Hessian formulae for the negative determinant of the leading principal
minor follows the same process as that for the determinant function, except that the matrix Q = I - vw^T
is not required. 
The formulae for f^1(P), g^1(P) and H^1(P) then reduce to

    f^1     =     -∏_i = 1^N-1 u_ii,
    
    g^1_ij    =     -f^1(P)(a_j^1)^Tb_i^1,
     H^1_[i,j],[k,ℓ]    =     g^1_kj(a_ℓ^1)^Tb_i^1 -
    g^1_ij(a_ℓ^1)^Tb_k^1,

where

    L̂b_i^1     =    e_i,
    Û^Ta_j^1     =    e_j.


In practice, the determinant of the leading principal minor is used rather than
that of the whole matrix. It is simpler, more efficient and the optimal value is independent of the
graph. It eliminates the need to scale any parameters by the size of the graph. When f^1(P) is used in lieu of f(P) the corresponding gradient vector and Hessian matrix are denoted by g^1(P) and H^1(P), respectively.




§ LU DECOMPOSITION-BASED EVALUATION ALGORITHM




The algorithm for computing f^1(P), g^1(P), H^1(P) is given here, along with the complexity of each step of the algorithm. Let k denote the average degree of the graph, that is, there are kN edges.



Input: P       

Output: f^1(P), g^1(P), H^1(P)       

       

begin        Complexity

1) Perform LU decomposition to find LU = I - P.         O(kN^2)

       

2) Calculate L̂ and Û, using (<ref>).        O(N)

       

3) Calculate each (a_j^1)^T and b_i^1, using (<ref>) and (<ref>).        O(N^3)

       

4) Calculate each (a_j^1)^Tb_i^1.        O(N^3)

       

5) Calculate f^1(P) = -∏_i=1^N-1 u_ii.        O(N)

       

6) Calculate each g_ij^1(P) = -f^1(P)(a_j^1)^Tb_i^1.        O(kN)

       

7) Calculate each H^1_[ij],[kℓ](P) = {[ g_kj^1(a_ℓ^1)^Tb_i^1 - g_ij^1(a_ℓ^1)^Tb_k^1,                                                                    i ≠ k j ≠ℓ i,j,k,l ≠ N;                                           0,                                                                                           ].        O(k^2N^2)

end       


Function evaluations algorithm




If the graph is sparse, the complexity of the above algorithm is O(N^3). However, for sufficiently dense graphs (that is, k > √(N)) the complexity of the above algorithm is O(k^2N^2). Note that each element of the Hessian is calculated in O(1) time, because they simply involve scalar multiplication
where all of the scalars have already been calculated in earlier steps, that is, the gradient terms in step 6, and
each (a_i^1)^Tb_l^1 in step 4.

These bounds are considerably better than the O(k^3N^4) bound that applies if we simply perform an LU decomposition for each element in the Hessian and gradient.
Consider the following six-node cubic graph Γ_6.







The adjacency matrix of Γ_6 is

    [[ 0 1 0 1 1 0; 1 0 1 0 0 1; 0 1 0 1 0 1; 1 0 1 0 1 0; 1 0 0 1 0 1; 0 1 1 0 1 0 ]].




Consider a point x such that,


    P(x) = [[   0 2/3   0 1/6 1/6   0; 2/3   0 1/6   0   0 1/6;   0 1/6   0 2/3   0 1/6; 1/6   0 2/3   0 1/6   0; 1/6   0   0 1/6   0 2/3;   0 1/6 1/6   0 2/3   0 ]].




Performing the LU decomposition of I - P using MATLAB's  routine we obtain matrices L and U (given to four decimal places)

    L = [[       1       0       0       0       0       0; -0.6667       1       0       0       0       0;       0 -0.3000       1       0       0       0; -0.1667 -0.2000 -0.7368       1       0       0; -0.1667 -0.2000 -0.0351 -0.5556       1       0;       0 -0.3000 -0.2281 -0.4444 -1.0000       1 ]],
    U = [[       1 -0.6667       0 -0.1667 -0.1667       0;       0  0.5556 -0.1667 -0.1111 -0.1111 -0.1667;       0       0  0.9500 -0.7000 -0.0333 -0.2167;       0       0       0  0.4342 -0.2412 -0.1930;       0       0       0       0  0.8148 -0.8148;       0       0       0       0       0       0 ]].











Consequently, L̂ = [[   L^Te_1        ⋯ L^Te_N-1      e_N ]]^T and
Û = [[   Ue_1      ⋯ Ue_N-1    e_N ]] are simply



    L̂ = [[       1       0       0       0       0       0; -0.6667       1       0       0       0       0;       0 -0.3000       1       0       0       0; -0.1667 -0.2000 -0.7368       1       0       0; -0.1667 -0.2000 -0.0351 -0.5556       1       0;       0       0       0       0       0       1 ]],Û = [[       1 -0.6667       0 -0.1667 -0.1667       0;       0  0.5556 -0.1667 -0.1111 -0.1111       0;       0       0  0.9500 -0.7000 -0.0333       0;       0       0       0  0.4342 -0.2412       0;       0       0       0       0  0.8148       0;       0       0       0       0       0       1 ]].










For all i,j, we calculate the a_j^1 and b_i^1 vectors using (<ref>) and (<ref>). Namely,



    a_1^1 = [[      1;    1.2; 0.2105; 1.0303; 0.6818;      0 ]], a_2^1 = [[      0;    1.8; 0.3158; 0.9697; 0.5455;      0 ]], a_3^1 = [[      0;      0; 1.0526; 1.6970; 0.5455;      0 ]], a_4^1 = [[      0;      0;      0; 2.3030; 0.6818;      0 ]], a_5^1 = [[      0;      0;      0;      0; 1.2273;      0 ]], a_6^1 = [[ 0; 0; 0; 0; 0; 1 ]],


    b_1^1 = [[      1; 0.6667;    0.2; 0.4474; 0.5556;      0 ]], b_2^1 =
    [[      0;      1;    0.3; 0.4211; 0.4444;      0 ]], b_3^1 =
    [[      0;      0;      1; 0.7368; 0.4444;      0 ]], b_4^1 =
    [[      0;      0;      0;      1; 0.5556;      0 ]], b_5^1 =
    [[ 0; 0; 0; 0; 1; 0 ]], b_6^1 =
    [[ 0; 0; 0; 0; 0; 1 ]].

We can now represent each (a_j^1)^Tb_i^1 as the ij-th element of the matrix

    [(a_j^1)^T(b_i^1)]_i,j=1^N,N = [[ 2.6818      2 1.2727 1.4091 0.6818      0;      2 2.5455 1.2727 1.2727 0.5455      0; 1.2727 1.2727 2.5455      2 0.5455      0; 1.4091 1.2727      2 2.6818 0.6818      0; 0.6818 0.5455 0.5455 0.6818 1.2273      0;      0      0      0      0      0      1 ]].

Then, f^1(P) = ∏_i=1^N-1û_ii≈ -0.1867. Note that we can directly verify the preceding by confirming that (A(P)) ≈ 1.1204 = -6(f^1(P)).

The gradient vector is then found using (<ref>). Note that we are only interested in the gradient elements for the eighteen arcs in the graph; this yields, to three decimal places:

    g^1(P) ≈[[ 0.374 0.263 0.127 0.374 0.238     0 0.238 0.374     0 0.263 0.374 0.127 0.127 0.127     0     0     0     0 ]].

Finally, the Hessian is found using (<ref>), given here to two decimal places:

    H^1(P) ≈[[     0     0     0  0.53  0.13     0     0 -0.41     0  0.11 -0.44 -0.09  0.02 -0.11     0     0     0     0;     0     0     0  0.11 -0.03     0  0.41     0     0  0.97  0.11  0.16  0.16     0     0     0     0     0;     0     0     0  0.02 -0.03     0  0.04 -0.11     0  0.16 -0.09     0  0.53  0.24     0     0     0     0;  0.53  0.11  0.02     0     0     0  0.13 -0.44     0     0 -0.41 -0.11     0 -0.09     0     0     0     0;  0.13 -0.03 -0.03     0     0     0  0.91  0.13     0  0.41     0  0.04  0.04 -0.03     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0;     0  0.41  0.04  0.13  0.91     0     0     0     0 -0.03  0.13 -0.03 -0.03  0.04     0     0     0     0; -0.41     0 -0.11 -0.44  0.13     0     0     0     0  0.11  0.53  0.02 -0.09     0     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0;  0.11  0.97  0.16     0  0.41     0 -0.03  0.11     0     0     0     0     0  0.16     0     0     0     0; -0.44  0.11 -0.09 -0.41     0     0  0.13  0.53     0     0     0     0 -0.11  0.02     0     0     0     0; -0.09  0.16     0 -0.11  0.04     0 -0.03  0.02     0     0     0     0  0.24  0.53     0     0     0     0;  0.02  0.16  0.53     0  0.04     0 -0.03 -0.09     0     0 -0.11  0.24     0     0     0     0     0     0; -0.11     0  0.24 -0.09 -0.03     0  0.04     0     0  0.16  0.02  0.53     0     0     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0;     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0 ]].





















Support for this work was provided by Australian Research Council (DP0666632 and DP0984470), the Office of Naval Research (Grant: N00014-02-1-0076) and the Army (Grant:  W911NF-07-2-0027-1). We would also like to thank V. Ejov for useful discussions and the referees whose comments corrected some errors and prompted us to improve the presentation.


9
detpaperV. Ejov, J.A. Filar, W. Murray and G.T. Nguyen. Determinants and longest cycles of graphs. SIAM Journal on Discrete Mathematics, 22(3):1215–1225, 2009.
krassJ.A. Filar and D. Krass. Hamiltonian cycles and Markov chains. Mathematocs of Operations Research, 19:223–237, 1994.
garey1976M.R. Garey, D.S. Johnson and R.E. Tarjan. The planar Hamiltonian circuit problem is NP-complete. SIAM Journal on Computing, 5(4):704–714, 1976.
mikeM. Haythorpe. Markov Chain Based Algorithms for the Hamiltonian Cycle Problem. PhD thesis, University of South Australia, 2010. Available at: http://www.stanford.edu/group/SOL/dissertations/michael-haythorpe-thesis.pdf.
hey1D. Heyman. A decomposition theorem for infinite stochastic matrices. Journal of Applied Probability, 32:893–901, 1995.
May1K. May. Derivatives of Determinants and Other Multilinear Functions. Mathematics Magazine, 38(5):207–208, 1965.

