On the Determinant and its Derivatives of the Rank-one Corrected Generator of a Markov Chain on a Graph

TL;DR

This paper introduces an efficient algorithm to compute the determinant and its derivatives of a rank-one corrected generator matrix in a Markov chain, facilitating solutions to the Hamiltonian cycle problem.

Contribution

It presents a novel method leveraging the matrix structure to compute all cofactors from a single LU decomposition, reducing computational complexity.

Findings

Efficient computation of derivatives for the generator matrix.

Reduction of cofactor calculations to a single LU decomposition.

Application to solving the Hamiltonian cycle problem.

Abstract

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.

Equations

$$p_{ij}(x)$$

$$A(P(x))$$

$$\min -\det (A(P(x)))$$

$$\sum _{j\in \mathcal{A}(i)}{x}_{ij}$$

$$\sum _{i\in \mathcal{A}(j)}{x}_{ij}$$

$$x_{ij}$$

$$\min \{f(P(x))|x\in \mathcal{D}\mathcal{S}\}.$$

$$g_{ii}$$

$$g_{ij}$$

$$Ge=G^{T}e$$

$$G=\left(\begin{array}[]{cc}1&{0}^{T}\\ {\ell}&I\end{array}\right)\left(\begin{array}[]{cc}1&{0}^{T}\\ {0}&\hat{G}\end{array}\right)\left(\begin{array}[]{cc}u_{11}&{u}^{T}\\ {0}&I\end{array}\right).$$

$$G=\left(\begin{array}[]{cc}u_{11}&{u}^{T}\\ u_{11}{\ell}&\;\;\hat{G}+{\ell u^{T}}\end{array}\right),$$

$$\left(\begin{array}[]{cc}1&{0}^{T}\\ {0}&\hat{G}\end{array}\right)\left(\begin{array}[]{cc}u_{11}&{u}^{T}\\ {0}&I\end{array}\right){e}={0}.$$

$${\hat{g}}_{ij}=g_{i+1,j+1}-{\ell }_i{u}_j,\forall i\ne j.$$

$$G^{T}e=U^T(L^{T}e)=0$$

$$\bar{G}=LU+e{e}_N^T=L(U+v{e}_N^T)=L\bar{U},$$

$$e^{T}Lv=e_N^{T}v=v_N=e^{T}e=N,$$

$$LU$$

$$Lv=e\text{mbox},and\bar{U}:=U+v{e}_N^T,$$

$$\text{hfill}\det \bar{U}=\prod _{i=1}^N{\bar{u}}_{ii}=N\prod _{i=1}^{N-1}{u}_{ii},$$

$$A$$

$${\bar{U}}^{T}w$$

$$A=L\left(\bar{U}+v\left[\frac{1}{N}{e}^T-e_N^T\right]\right)=L\left(\bar{U}+v{w}^T\bar{U}\right)=L\left(I+v{w}^T\right)\bar{U}.$$

$$\det (A)$$

$$\det (I+c{d}^T)$$

$$\det (I+v{w}^T)$$

$$w^{T}v$$

$$w^T$$

$$v$$

$$w^{T}v$$

$$(I-P+e{e}_N^T{)}^{-1}e$$

$$w^{T}v=\left(\frac{1}{N}{e}^T-e_N^T\right)e=0,$$

$$\det (A(P))$$

Full text

11institutetext: J.A. Filar 22institutetext: Flinders University

22email: jerzy.filar@flinders.edu.au 33institutetext: M. Haythorpe 44institutetext: Flinders University

44email: michael.haythorpe@flinders.edu.au 55institutetext: W. Murray 66institutetext: Stanford University

66email: 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

1 Introduction

The Hamiltonian cycle problem (HCP) is an important graph theory problem that features prominently in complexity theory because it is NP-complete garey1976 . 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 $\Gamma$ 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 $\Gamma$ contains at least one Hamiltonian cycle (HC), we say that $\Gamma$ is a Hamiltonian graph. Otherwise, we say that $\Gamma$ 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)\in\Gamma$, and solve an associated optimisation problem. A convenient method of representing these constraints is to use a matrix $P({x})$, where

[TABLE]

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 krass ), continuous optimisation problems that are equivalent to HCP have been discovered in recent times. In particular, it was demonstrated in detpaper that if we define

[TABLE]

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

[TABLE]

subject to

[TABLE]

where $\mathcal{A}(i)$ is the set of nodes reachable in one step from node $i$. Constraints (5)–(7) are called the doubly-stochastic constraints. For neatness, we refer to constraints (5)–(7) as the set $\mathcal{DS}$, and we call the objective function $f(P({x}))$. Then, the above problem can be represented as follows:

[TABLE]

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 mike ) 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 $\mathcal{DS}$. 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 $\mathcal{DS}$ is simple and does not involve evaluating the objective or its derivatives.

2 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 $j$th, which is unity. We will often denote the $(i,j)$th elements of a matrix, say $W$, by $w_{ij}$.

Definition 1

A matrix $G$ is said to have property $S$ if

[TABLE]

Theorem 2.1

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

Proof

Following common convention the diagonal elements of $L$ are 1. There is no loss of generality if we assume that $g_{11}\neq 0$. After one step of Gaussian elimination (GE) we obtain

[TABLE]

Note that, expanding this matrix multiplication gives

[TABLE]

and therefore {u^{T}}=\left[\begin{array}[]{cccc}g_{12}&\;g_{13}&\cdots&\;g_{1N}\end{array}\right] and {\ell}=\displaystyle\frac{1}{u_{11}}\left[\begin{array}[]{cccc}g_{21}&\;g_{31}&\cdots&\;g_{N1}\end{array}\right]^{T}. Clearly, $u_{11}=g_{11}\geq 0$, and therefore ${u}^{T}\leq 0$ and ${\ell}\leq 0$. We have $G{e}={0}$, which implies

[TABLE]

It follows that $\hat{G}{e}=\hat{G}^{T}{e}={0}$. By definition we have

[TABLE]

Since $l_{i}\leq 0$ and $u_{j}\leq 0$ it follows that $\hat{g}_{ij}\leq 0,\;\forall i\neq j$. From this result and $\hat{G}{e}={0}$ it follows that $\hat{g}_{ii}\geq 0$ and that $\hat{G}$ has property $S$. We can now proceed with the next step of GE. Note that if $\hat{g}_{11}=0$ we must have the first row and column of $\hat{G}$ 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 $\hat{G}$ is nonzero.

Lemma 1

Regardless of the rank of $G$ we have $u_{NN}=0$.

Proof

After $N-1$ steps of GE $\hat{G}$ is a $1\times 1$ matrix. The only $1\times 1$ matrix that satisfies property $S$ is [math].

Corollary 1

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

Proof

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

Lemma 2

If $G$ satisfies property $S$ and $G=LU$, the off-diagonal elements of $L$ and $U$ are non-positive. Moreover, $-1\leq l_{ij}\leq 0$ for $i\neq j$.

Proof

The first part follows immediately from the fact that ${u}^{T}\leq 0$ and ${l}\leq 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 $\hat{G}^{T}e=0$ and the off-diagonal elements in the column of $\hat{G}$ 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$.

Lemma 3

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

Proof

Since $U$ has rank $N-1$ and $u_{NN}=0$, the basis for the nullspace of $U^{T}$ is $\kappa e_{N}$, where $\kappa\neq 0$. From property $S$,

[TABLE]

and hence we obtain $L^{T}e=\kappa e_{N}$, for some $\kappa$. Since $L$ has unit diagonal entries, we have $\kappa=1$, which yields the result.

Lemma 4

When $G$ satisfies property $S$ and is rank $N-1$ then $\bar{G}\equiv G+{e}{e}_{N}^{T}$ is nonsingular.

Proof

We have

[TABLE]

where $Lv=e$. Note that $\bar{U}=U+ve^{T}_{N}$ is upper triangular. Moreover, the $(N,N)$th element of $\bar{U}$ is $v_{N}$. Since

[TABLE]

it implies that $\bar{U}$, and hence $\bar{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 hey1 . In the context of this paper, we note that when $P(x)$ corresponds to the strict interior of $\mathcal{DS}$, irreducibility applies.

3 Computing $\det(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=\Pi A$, where $\Pi$ 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 (8) 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 $\mathcal{DS}$ is simple and does not involve evaluating the objective or its derivatives. It is known that the global minimiser of (8) is $-N$. Indeed, from Theorem 4.1, Proposition 4.4 and Proposition 4.6 in Ejov et al detpaper , we know that $-\det(A(P(x)))\in[-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\in\mathcal{DS}$ 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\in\mathcal{DS}$, the nonsingularity of $A(P(x))$ implies that $I-P$ has rank $N-1$.

3.1 Product forms of $A(P)$ and $\det(A(P))$

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

[TABLE]

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

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

We show that two of these factors have a determinant of 1. 3. (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. (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\times 1$ vector, and $\bar{U}$ be an $N\times N$ matrix, such that

[TABLE]

where {e}_{N}^{T}=\left[\begin{array}[]{cccc}0&\cdots&0&1\end{array}\right]. Since $L$ is nonsingular, ${v}\neq 0$ exists and is unique, and therefore $\bar{U}$ is well-defined. The first $N-1$ columns of $\bar{U}$ are identical to those of $U$. Consequently $\bar{U}$ is also upper triangular. Since $LU=G$ satisfies property $S$, it follows from Lemma 1 that $u_{NN}=0$, and therefore from (9) that $\bar{u}_{NN}=v_{N}=N$. Therefore,

[TABLE]

where $u_{ii}$ is the $i$-th diagonal element of $U$. Exploiting (11) we may write

[TABLE]

Since $\bar{U}$ is nonsingular we may define ${w}$ to be the unique solution to the system

[TABLE]

Then, from (13)

[TABLE]

We take the determinant of (15) to obtain

[TABLE]

Note that, for any vectors ${c}$ and ${d}$,

[TABLE]

This is because ${c}{d}^{T}$ has one eigenvalue ${d}^{T}{c}$ of multiplicity 1 and an eigenvalue 0 of multiplicity $N-1$. Consequently,

[TABLE]

which we substitute into expression (16) above.

Lemma 5

The inner-product ${w}^{T}{v}$ in (18) satisfies

[TABLE]

Proof

From their respective definitions (14) and (11),

[TABLE]

Then, from (19)–(20) we obtain

[TABLE]

Since $I$, $P$ and ${e}{e_{N}}^{T}$ are all stochastic matrices, we know that $I-P+{e}{e}_{N}^{T}$ has row sums of 1 as well. Hence, its inverse also has row sums equal to 1, that is,

[TABLE]

Substituting (22) into (21), we obtain

[TABLE]

which concludes the proof.

We now derive the main theorem of this subsection.

Theorem 3.1

Let $LU$ denote the LU decomposition of $I-P$ and $u_{ii}$ be the diagonal elements of $U$ then

[TABLE]

Proof

From (16), (18) and Lemma 5 we know that

[TABLE]

From the construction of the LU decomposition we know that $\det(L)=1$ and using (12) gives

[TABLE]

This concludes the proof.

Remark 1

Note that finding ${v}$ and ${w}$ is a simple process because $L$ and $\bar{U}^{T}$ are lower-triangular matrices, so we can solve the systems of linear equations in (14) and (11) directly.

3.2 Finding the gradient g(P)

Next we use the LU decomposition found in Subsection 3.1 to find the gradient of $f(P)=-\det(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):=\frac{\partial f(P)}{\partial x_{ij}}$ for each $x_{ij}$ such that $(i,j)\in\Gamma$.

Consider vectors $a_{j}$ and $b_{i}$ satisfying the equations $\bar{U}^{T}a_{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-{v}{w}^{T}$, where ${v}$ and ${w}^{T}$ are as in (19)–(20). We prove the following result in this subsection:

[TABLE]

where $g_{ij}(P)$ is the gradient vector element corresponding to the arc $(i,j)\in\Gamma$.

The outline of the derivation of (23) is as follows.

1. (1)

We represent each element $g_{ij}(P)$ of the gradient vector as a cofactor of $A(P)$. 2. (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. (3)

We then express the element $g_{ij}(P)$ of the gradient vector as the product of $\det(A(P))$ and the determinant of the elementary matrix, the latter of which is shown to be equal to $a_{j}^{T}Qb_{i}$.

For any matrix $V=\left(v_{ij}\right)_{i,j=1}^{N,N}$ it is well-known (e.g., see May May1 ) that $\displaystyle\frac{\partial\det(V)}{\partial v_{ij}}=(-1)^{i+j}\det(V^{ij})$, where $V^{ij}$ is the $(i,j)$-th minor of $V$. That is, $\displaystyle\frac{\partial\det(V)}{\partial v_{ij}}$ is the $(i,j)$-th cofactor of $V$. Since the $(i,j)$-th entry of $A(P)$ is simply $a_{ij}=\delta_{ij}-x_{ij}+\frac{1}{N}$ (where $\delta_{ij}$ is the Kronecker delta that is 1 if $i=j$, and 0 otherwise), it now follows that

[TABLE]

However, rather than finding the cofactor we calculate (24) by finding the determinant of a modification of $A$ in which the $i$th 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}$, $\ldots$, $r_{N}^{T}$. Then, we perform the following row operation:

[TABLE]

where $\alpha_{j}(i)$ is the $i$-th element of vector $\alpha_{j}$ and $A^{T}\alpha_{j}={e}_{j}$. This row operation replaces the $i$-th row of $A$ with $\alpha_{j}^{T}A={e}_{j}^{T}$, as desired.

In this case, from (15), $A^{T}=\bar{U}^{T}(I+{w}{v}^{T})L^{T}$. Since $A$ is nonsingular $\alpha_{j}$ can be found directly:

[TABLE]

Lemma 6

[TABLE]

Proof

Consider

[TABLE]

Therefore, $\left(I+{w}{v}^{T}\right)^{-1}=\left(I-{w}{v}^{T}\right).$

Taking the above result and substituting into (26), we obtain

[TABLE]

Next, we define an elementary matrix $E_{ij}$ by

[TABLE]

and note that it performs the desired row operation (25) on A because

[TABLE]

in effect replacing the $i$-th row of $A$ with ${e}_{j}^{T}$. Therefore,

[TABLE]

From (28), we rewrite $E_{ij}=I-{e}_{i}({e}_{i}-\alpha_{j})^{T}$. Then, from (17) we obtain

[TABLE]

Substituting (27) into (30) we obtain

[TABLE]

For convenience we define $Q:=I-{v}{w}^{T}$. Then

[TABLE]

We now derive the main result of this subsection.

Proposition 1

The general gradient element of $f(P)$ is given by

[TABLE]

Proof

Substituting (32) into (29) immediately yields the result.

Remark 2

Note that we can calculate all $a_{j}$ and $b_{i}$ in advance, by solving the systems of linear equations in (32), again in reduced row echelon form. Then, for the sake of efficiency we first calculate

[TABLE]

and then calculate

[TABLE]

This allows us to rewrite the formula for $g_{ij}(P)$ as

[TABLE]

3.3 Finding the Hessian matrix H(P)

Here, we show that the LU decomposition found in Subsection 3.1 can also be used to calculate the Hessian of $f(P)$ efficiently. Consider $g_{ij}$ and $\hat{q}_{ij}$ as defined in (33) and (35) respectively. We prove the following result in this subsection:

[TABLE]

where $H_{[ij],[k\ell]}$ is the general element of the Hessian matrix corresponding to arcs $(i,j)$ and $(k,\ell)\in\Gamma$.

The outline of the derivation is as follows.

1. (1)

We represent each element $H_{[ij],[k\ell]}(P)$ of the Hessian matrix as a cofactor of a minor of $A(P)$. 2. (2)

We construct a second elementary matrix that in conjunction with $E_{ij}$ (see (28)) transforms matrix $A(P)$ into one with a determinant equivalent to the $(k,\ell)$-th cofactor of $A^{ij}(P)$. 3. (3)

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

Using results obtained from finding $g(P)$ in Subsection 3.2, we obtain these values immediately.

We define $A^{[ij],[k\ell]}$ to be the matrix $A$ with rows $i,k$ and columns $j,\ell$ removed. An argument similar to that for $g_{ij}(P)$ in the previous subsection can be made that finding

[TABLE]

is equivalent to finding the negative determinant of $A$ with the $i$th and $k$th rows changed to ${e}_{j}^{T}$ and ${e}_{\ell}^{T}$ respectively. That is,

[TABLE]

where $\hat{E}_{k\ell}$ is an additional row operation constructed to change row $k$ of $E_{ij}A$ into ${e}_{\ell}^{T}$. Note that if $i=k$ or $j=\ell$, the matrix $A^{[ij],[k\ell]}$ is no longer square and the determinant no longer exists. If this occurs, we define $H_{[ij],[k\ell]}:=0$. If both $i=k$ and $j=\ell$, we also define $H_{[ij],[k\ell]}:=0$, as the determinant is linear in each element of $A(P)$.

Consider $E_{ij}A$ composed of rows $\hat{r}_{1}^{T}$, $\hat{r}_{2}^{T}$, $\ldots$, $\hat{r}_{N}^{T}$. Then, we perform the following row operation:

[TABLE]

where $(E_{ij}A)^{T}\gamma_{\ell}={e}_{\ell}.$ Then, similarly to (26), we directly find $\gamma_{\ell}$:

[TABLE]

Next, in a similar fashion to (28), we construct an elementary matrix $\hat{E}_{k\ell}$

[TABLE]

Then, we evaluate $\det(\hat{E}_{k\ell})$:

[TABLE]

Recall from (28) that $E_{ij}=I-{e}_{i}({e}_{i}^{T}-\alpha_{j}^{T})$. We have

[TABLE]

Recall from (29) that $g_{ij}=\det(A)\det(E_{ij})$, and from (30) that $\alpha_{j}^{T}{e}_{i}=\det(E_{ij})\neq 0$, and therefore (42) holds. Then,

[TABLE]

Substituting (43) into (42) we obtain

[TABLE]

and further substituting (44) into (41), we obtain

[TABLE]

Note that since $i\neq k$, ${e}_{i}{e}_{i}^{T}{e}_{k}={0}$, and from (43), $\alpha_{j}^{T}{e}_{k}=\displaystyle\frac{g_{kj}}{\det(A)}$. Hence, from (45) and (32) we obtain

[TABLE]

We now derive the main result of this subsection.

Proposition 2

The general element of the Hessian of $f(P)$ is given by

[TABLE]

where $\hat{q}_{i\ell}$ and $\hat{q}_{k\ell}$ are defined in (35).

Proof

From (37) and (29), we can see that $H_{[ij],[k\ell]}=-\det(\hat{E}_{k\ell})g_{ij}$. Then, from (46), $\det(\hat{E}_{k\ell})=a_{\ell}^{T}Q(b_{k}-\displaystyle\frac{g_{kj}}{g_{ij}}b_{i})$ and so $H_{[ij],[k\ell]}=-a_{\ell}^{T}Q(b_{k}g_{ij}-b_{i}g_{kj})$.

In order to improve computation time, we take advantage of the fact that we evaluate every $\hat{q}_{ij}$ while calculating the gradient to rewrite the second order partial derivatives of $f(P)$ as

[TABLE]

This concludes the proof.

Remark 3

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)(\hat{q}_{jk})$, which remains valid despite arc $(k,j)$ not appearing in the graph.

3.4 Leading principal minor

It is, perhaps, interesting that instead of using the objective function $f(P)=-\det\left(I-P+\frac{1}{N}{e}{e}^{T}\right)$, it is also possible to use $f^{1}(P):=-\det(G^{NN}(P))$, the negative determinant of the leading principal minor of $I-P$. The following, somewhat surprising, result justifies this claim.

Theorem 3.2

1. (1)

$f^{1}(P)=\displaystyle\frac{1}{N}f(P)=-\displaystyle\frac{1}{N}\det\left(I-P+\displaystyle\frac{1}{N}{e}{e}^{T}\right)$. 2. (2)

If the graph is Hamiltonian, then

[TABLE]

Proof

First, we show part (1), that is, $f^{1}(P)=\displaystyle\frac{1}{N}f(P)$. To find $f^{1}(P)$, we construct $LU=I-P$ as before, and define $\hat{L}$, $\hat{U}$ as:

[TABLE]

That is, $\hat{L}$ is the same as $L$ with the last row replaced by ${e}_{N}^{T}$, and $\hat{U}$ is the same as $U$ with the last column replaced with ${e}_{N}$. Then consider

[TABLE]

Since $L$ is lower-triangular, ${e}_{i}^{T}L{e}_{N}=0$ for all $i\neq N$. Likewise, since $U$ is upper-triangular, ${e}_{N}^{T}U{e}_{j}=0$ for all $j\neq N$. Therefore the above matrix simplifies to

[TABLE]

which is the same as $LU$ with the bottom row and rightmost column removed, and a 1 placed in the bottom-right element. Therefore, $\det(\hat{L}\hat{U})=\det(G^{NN}(P))$, and consequently

[TABLE]

Note that $\hat{L}$ and $\hat{U}$ are triangular matrices, so

[TABLE]

However, only the last diagonal elements of $\hat{L}$ and $\hat{U}$ are different from $L$ and $\bar{U}$ (see (11)) respectively, so

[TABLE]

Now, since $\hat{l}_{NN}=l_{NN}=1$, we have

[TABLE]

We also have $\hat{u}_{NN}=1$, but by Lemma 4, $\bar{u}_{NN}=N$ and hence

[TABLE]

Therefore, substituting (73) and (74) into (71) we obtain

[TABLE]

Therefore, part (1) is proved.

The proof of part $(2)$ of Theorem 3.2 follows directly from the fact that $\min f(P)=-N$ (proved in detpaper ), and part $(1)$.

Remark 4

Using the leading principal minor has the advantage that the rank-one modification $\frac{1}{N}{e}{e}^{T}$ is not required, which makes calculating the gradient and the Hessian even simpler than described in Subsection 3.2 and Subsection 3.3 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-{v}{w}^{T}$ is not required.

The formulae for $f^{1}(P)$, $g^{1}(P)$ and $H^{1}(P)$ then reduce to

[TABLE]

where

[TABLE]

Remark 5

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.

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

[TABLE]

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>\sqrt{N}$) the complexity of the above algorithm is $O(k^{2}N^{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 $\left(a_{i}^{1}\right)^{T}b_{l}^{1}$ in step 4.

These bounds are considerably better than the $O(k^{3}N^{4})$ bound that applies if we simply perform an LU decomposition for each element in the Hessian and gradient.

Example 1

Consider the following six-node cubic graph $\Gamma_{6}$.

The adjacency matrix of $\Gamma_{6}$ is

[TABLE]

Consider a point ${x}$ such that,

[TABLE]

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

[TABLE]

Consequently, \hat{L}=\left[\begin{array}[]{cccc}L^{T}{e}_{1}&\cdots&L^{T}{e}_{N-1}&{e}_{N}\end{array}\right]^{T} and \hat{U}=\left[\begin{array}[]{cccc}U{e}_{1}&\cdots&U{e}_{N-1}&{e}_{N}\end{array}\right] are simply

[TABLE]

For all $i,j$, we calculate the $a_{j}^{1}$ and $b_{i}^{1}$ vectors using (78) and (79). Namely,

[TABLE]

[TABLE]

We can now represent each $(a_{j}^{1})^{T}b_{i}^{1}$ as the $ij$-th element of the matrix

[TABLE]

Then, $f^{1}(P)=\prod\limits_{i=1}^{N-1}\hat{u}_{ii}\approx-0.1867$. Note that we can directly verify the preceding by confirming that $\det(A(P))\approx 1.1204=-6(f^{1}(P))$.

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

[TABLE]

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

[TABLE]

Acknowledgements.

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.