Embedding Conjugate Gradient in Learning Random Walks for Landscape Connectivity Modeling in Conservation

1 Introduction

With many animal species experiencing a population decline due to habitat loss and fragmentation, efforts to conserve wildlife and to protect biodiversity are forefront challenges for computational sustainability [Gomes, 2009]. The difficulty of the task intensifies with the necessity to accurately predict the effects of broad-scale ecological processes from different actions implemented by conservation efforts [McRae and Beier, 2007].

Towards these set of challenges, artificial intelligence has been employed to understand the relationship between environmental variables and animal movement behavior. Early approaches in this domain were based on maximizing graph-theoretic notions such as the least-cost paths, multi-commodity flows, steiner trees, etc. [Sheldon et al.,

2010; Bras et al., 2013; Dilkina et al., 2013] but fall short in capturing the stochasticity of animal decision-making.

A newer model called the Circuitscape [McRae et al., 2008] incorporates circuit theory and random walks to model animals as individual agents who traverse to neighboring land parcels with probabilities proportional to resistance values, i.e., the cost of moving. However, existing methods indirectly learn random walk parameters; these parameters are derived from separate statistical models that estimate the landscape cost [Inman et al., 2013; McClure et al., 2016], which can lead to inaccurate random walk models that may not resemble true movement dynamics. Moreover, ecologists have stated a dire need to develop automated procedures to weigh different environmental factors in constructing Circuitscape models [Dickson et al., 2019].

We develop a principled computational framework that learns random walk models directly from data. Our core idea is a moment matching strategy, which correlates the model's statistical properties, such as the hitting and commuting times, with those observed empirically. By exploiting the connection of random walks with circuit theory, we show that learning the probabilistic process can be reduced to finding the optimal graph Laplacian matching empirical observations. The graph Laplacian is the central notion that connects statistical properties of random walks with concepts of voltages, currents, and resistances. Following a gradient descent strategy, our forward pass to compute the hitting and commuting times requires inverting the graph Laplacian. The backward pass, therefore, requires an efficient computation scheme to flow the gradients through the matrix inverse.

Computationally, we propose a neural network capable of back-propagating gradients through the matrix inverse in an end-to-end fashion. Our first idea, DirectInv, was developed through the implicit function theorem, motivated by the work of [Amos and Kolter, 2017; Ling et al., 2019]. Nevertheless, the direct application of the implicit function theorem leads to a Kronecker computation, which is quartic in the size of the graph Laplacian. We then developed CGInv, which back-propagates the gradients through a neural network encoding each layer as one conjugate gradient iteration. While the time complexity of full CGInv is lower than DirectInv, CGInv can be further expedited by initializing the current execution with the root from the previous one, which is relatively close to the correct solution. This decreases the number of iterations to convergence and further speeds up the forward and backward computations.

To showcase our computational framework, we designed two sets of experiments to learn random walk models from animal movement data. In the first experiment, we simulated data from Dilkina et al. (2013) and learnt optimal models from empirical observations. Our second experiment learned landscape-resistance values by correlating them with genetic similarities among animal individuals in different locations via the Mantel test. In both applications, our approach succeeds in recovering the ground-truth animal movement patterns. Our contributions can be summarized in the following way:

1. We formulated an end-to-end computational framework using a moment matching strategy to learn random walks from data.
2. We developed an efficient and powerful technique called CGInv to back-propagate gradients through the matrix inverse, via embedding each conjugate gradient step as a neural network layer.
3. We successfully demonstrated the effectiveness of our framework in recovering the ground-truth animal movement patterns from realistic datasets.

2 RandomWalks and Circuit Theory

We describe the behavior of a single animal individual as a random walker traversing over a predefined landscape. The landscape is modeled as an undirected graph G = (V,E) comprised of a node set V representing various locations and an edge set E representing direct, physical connections between two locations. The animal individual is fitted with unique transition probabilities which are species-specific. The transition probability pu,v from node u to node v is encoded in edge (u, v) ∈ E through an edge weight.

Circuit Theory. The Circuitscape model builds an electrical circuit described by the graph G. The edge weights feature resistance values, ru,v, that reflect the cost to travel along a given edge (u,v) from u to v. The reciprocal of resistance, 1/ru,v, yields a conductance value denoted as cu,v. We connect the Circuitscape model with the random walk model by letting the transition probability of a random walker at node u traveling to node v be cu,v/Cu where Cu = �v∈N(u) cu,v represents the sum of edge conductance values from all the neighboring nodes v ∈ N(u) of node u.

Ohm's Law. The Ohm's law relates voltages, currents, and resistances in the Circuitscape model through a set of linear equations. Let CG be the conductance matrix, where its u, v-th entry cu,v represents the conductance value of edge (u, v). Let DG be a diagonal matrix, where the u-th diagonal entry is Cu, the sum of the conductance values of the edges incident to node u. The graph Laplacian of G is defined as LG = DG-CG. Denote i = (i1, . . . , in)T as the vector of currents, in which ik is the amount of current that has been injected to node k 1 . Denote v = (v1, . . . , vn)T as the vector of voltages, in which vk is the voltage value of node k 2 . Then, Ohm's Law can be written in the following matrix form:

1 Current that flows out of a node is denoted by negative values

An important extension of this principle is used in defining effective resistance, r e f f u,v , between two nodes u and v. The r e f f u,v measures the resistance between u and v by injecting 1 ampere of current into node u, extracting 1 ampere from node v, and measuring V(u)-V(v). Mathematically, let v st be the voltage vector corresponding to the current vector i st = 1s -1t with only 1 at the s and t-th entries and 0 elsewhere. Applying Ohm's law , the effective resistance r e f f s,t is exactly the voltage difference between the s-th and the t-th entry of v st , i.e.,

Here L-1 G represents the pseudo-inverse. In this paper, we assume that graph G is always connected. Therefore LG has exactly one zero eigenvalue, corresponding to the eigenvector (1, . . . ,1)T implying LG to be highly ill-conditioned. When we perform matrix inversions in this paper, we remove the singularity of LG by grounding a single node in the voltage vector to 0. Then, we remove the row and column in LG that corresponds to the grounded node. This reduction yields in a well-conditioned, invertible matrix that can be handled.

Hitting and Commuting Times. An important notion in random walks is the hitting time. The hitting time Huv from u to v represents the expected number of steps for a random walker on G to reach v starting from u. From the hitting time, the commuting time Juv is defined as Juv =Huv +Hvu, which represents the expected number of steps to start from u, travel to v, and come back at u. The following two theorems connect effective resistance and Ohm's Law to the commuting and hitting times respectively.

Theorem 1. [Chandra et al., 1996] Let F = �x,y 1/rx,y = �x,y cx,y. The commuting time for a random walk between node s and t is equal to F r e f f s,t .

Theorem 2. [Chandra et al., 1996] Suppose we inject �y∈N(x) 1 rxy =�y∈N(x) cxy amount of current to each node x ∈ V and extract �x∈V �y∈N(x) 1 rxy amount of current from node v, then the voltage difference between node u and node v, Fuv is exactly the hitting time between u, v, Huv. In mathematical notation, Huv = (1u -1v)T L-1 G iuv, where the v-th entry of iuv is �x̸=v �y∈N(x) cxy, and the x-th entry of iuv (x ̸= v) is �y∈N(x) cxy.

3 Learning RandomWalks

Learning the migration of wild-life animals is crucial towards designing effective conservation strategies. Maintaining landscape connectivity allows animals to roam freely, promoting the dispersal of their gene flows and increasing the robustness of the species against environmental changes.

2 All voltages should be referenced with respect to a given node, for example, node 0. In other words, we assume that node 0 is grounded.

To understand landscape connectivity, the Circuitscape model was developed to use circuit theory and the concept of effective resistance to accurately represent landscape-to-resistance costs [McRae, 2006; McRae et al., 2008]. While circuit-theory based models are powerful, many approaches in inferring landscape connectivity parameters rely on expert designed values or subjective methods in evaluating the contribution of different environmental factors.

To combat these challenges, we developed a principled framework whose foundation lies in the method of moments technique to precisely estimate landscape-to-resistance values from data. We use method of moments approach to build the appropriate random walk model by finding resistance values that best mimic animal movement patterns.

We first learn a Circuitscape model representing animals' dispersal patterns based on partially observed trajectories. As a second application, we learn landscape features by modeling gene flow as done in previous work [Cushman et al., 2006; McRae and Beier, 2007]. We correlate the genetic distances of animal individuals observed across a landscape with their migratory patterns using the Circuitscape model. We show that both applications can be solved with a single computational framework.

3.1 Learning RandomWalks fromPartial Trajectories

As animals disperse across their habitats, they arrive at different observable locations at different times. Since the path taken by the animal between two observed locations is unknown, the time series data represents partial trajectory data as a sequence of time-location pairs:

Here, a single animal individual visits location l1 at time t1, then visits location l2 at time t2, and repeats this process until time tk. The behavior learning problem is: how can we identify the probabilistic model of the random walk given partial trajectory data listed as above?

We adopt the method of moments technique to learn the true probability distribution by matching the empirically observed hitting times with those generated from the Circuitscape model. From data, we can observe instances of the hitting time between any two nodes u and v. An instance of the hitting time represents the number of steps an animal individual was observed to travel starting from u and arriving at v for the first time. The following is a more rigorous definition:

The time tu represents the animal's timestamp at location lu and the time tv represents the animal's timestamp at location lv. As the number of observations increase, the true expected hitting time should be well approximated by the empirically computed hitting times.

Therefore, the core idea is to match the expected hitting times computed by the Circuitscape model to the N instances of the hitting times observed empirically:

This optimization equation translates into finding the optimal L∗G that minimizes the sum of the squared differences for all N instances of the hitting time observed in the data.

3.2 Mantel Coefficient Optimization

A fundamental challenge within ecology and conservation biology is to use genetic approaches to understand the dispersal of animal species over long periods of time [Cushman et al., 2006; McRae and Beier, 2007]. One approach employed by Cushman et al. (2006) is to manually find landscape-resistance costs through Mantel tests among the genetic distances between pairs of animal individuals and their least-cost paths. Nevertheless, they used a brute-force approach to consider all combinations of environmental features. We extend the study by Cushman et al. (2006) by providing a principled way to find the optimal combination of environmental features via correlation of the genetic distances with commuting times inferred from the resistance surface.

Here, we assume as input a genetic distance matrix that describes the genetic similarity between all pairs of animal individuals across a landscape. Similar to Cushman et al. (2006), we use the Bray-Curtis percentage dissimilarity measure of genetic distance between individuals where 0 indicates identical composition and 1 indicates completely different individuals. Then, we compute the commuting times between all pairs of individuals using the Circuitscape model and find the LG that maximizes the Mantel correlation coefficient with the genetic distance matrix [Legendre and Legendre, 2012]. This translates to finding the graph Laplacian LG which maximizes the following objective:

Here, Gkl refers to specific indexed locations within the genetic distance matrix G. Similarly, Jkl refers to specific matrix elements within the commuting time matrix J. The mean and standard deviation of the genetic matrix are ¯G and sG respectively. The mean and standard deviation of the commuting time are ¯ J and sJ respectively. Since G and J are symmetric, Equation (4) only operates on the upper-triangular section of both matrices where d = n(n1)/2 is the number of items in the sum.

3.3 Training via Stochastic Gradient Descent

Analytically solving the optimization problems in Equations (3) and (4) is infeasible due to over-constrained linear systems. Gradient-based learning offers a lucrative alternative that yields a straightforward methodology. In both scenarios, the gradients from the optimization equations need to propagate back to the physical resistance values to perform stochastic gradient descent. We can see that the hitting times Hli,li+1 from Equation (3) and commuting time terms Jkl , ¯ J, and sJ from Equation (4) both rely on the L-1 G from Theorems (2) and (1) respectively. We must propagate one step further to the LG to reach the conductance values, which are reciprocals of resistance. By finding the gradients through LG, we can then implement a stochastic gradient descent framework to perform updates on the resistance values.

Figure 1: (Left, orange and blue box) The conjugate gradient iteration to find the root for Ax = b. Every operation in one iteration is fully differentiable. One conjugate gradient iteration takes the input of rk, pk and xk, and computes rk+1, pk+1 and xk+1 as the output for the next iteration. After n iterations, xn+1 is guaranteed to be the root of Ax = b. (Right, red dashed box) Back-propagate gradient through one conjugate gradient iteration. Here, we assume having access to the partial derivative of the loss function w.r.t. the output rk+1, pk+1 and xk+1, ¶ l ¶rk+1 , ¶l ¶ pk+1 and ¶l ¶xk+1 . We compute ¶l ¶rk , ¶l ¶ pk and ¶l ¶xk w.r.t. the input. The computation is dominated by matrix vector multiplication, which scales O(n 2 ). The last line shows ¶l ¶A can be computed via ¶l ¶rk and ¶l ¶ak .

4 Propagating Gradients through the Matrix Inverse Efficiently

Since hitting and commuting times rely on inverting LG in Ohm's Law in Equation (1), a key ingredient in the stochastic gradient descent framework is to compute the following partial derivative:

In our machine learning framework, successfully computing this quantity dictates the efficiency of the gradient descent methodology. We attacked this problem from two different angles. We initially started our approach, DirectInv, by using the implicit function theorem to derive a gradient rule. We iterated upon this idea to develop our second approach, CGInv, which inverts the process of conjugate gradient. The following two sections detail both approaches.

4.1 DirectInv using Implicit Function Theorem

To compute the ¶A-1 ¶A , we reformulate the matrix A and A-1 into vectors:

Here, ki, j is the i, j-th entry of the matrix A and li, j is the i, j-th entry of A-1 . We define ¶A-1/¶A as a n 2 -by-n 2 matrix, in which the (r,s)-th entry is ¶vec(A-1 )r ¶vec(A)s . In other words, the ((i-1)n+ j,(u-1)n+v)-th entry of the matrix ¶A-1 /¶A is ¶A-1 u,v/¶Ai, j. We can derive ¶A-1 /¶A in a closed form, reflected by the following theorem:

Theorem 3. Suppose matrix A is full rank, then

Here, ⊗ means the Kronecker product. Our derivation is based on the following implicit function theorem:

Theorem 4. Let f : R n+m → R m be a continuously differentiable function, and let R n+m have coordinates (x,y). Fix a point (a,b) = (a1, . . . ,an,b1, . . . ,bm) with f (a,b) = 0, where 0 ∈ R m is the zero vector. If the Jacobian matrix Jf ,y(a,b) is invertible, then there exists an open set U of R n containing a such that there exists a unique continuously differentiable function g : U →R m such that

and

Moreover, the partial derivatives of g in U are given by

Through this derivation, we can see that the ¶L-1 G /¶LG =

-L-1 G ⊗ (L-1 G )T can be computed through the Kronecker product. DirectInv produces a gradient computation that requires O(n 4 ) multiplication operations to compute the Kronecker product. Due to its large time complexity, DirectInv is limited to application on small-scale instances.

4.2 CGInv: Inverting Conjugate Gradient

Our second approach, CGInv, relies on back-propagating through a neural network representing conjugate gradient, motivated from the properties of LG. The structure of LG is symmetric, diagonally-dominant, and positive semi-definite. Additionally, LG is quite sparse due to the limited number of edges incident on each node. Therefore, a matrix inversion technique needs to only work with non-zero elements of the matrix. By removing the row and column corresponding to the node with grounded voltage, LG is a well-conditioned, positive-definite matrix suitable for conjugate gradient. An illustration of conjugate gradient is shown in Figure 1. The method of conjugate gradients is initialized with a vector x0 that serves as a guess. Following the iterative procedure as shown in the blue box of Figure 1, the k-th conjugate gradient iteration uses the current approximation of the root xk and refines the solution to xk+1. Essentially, conjugate gradient modifies the solution vector by taking correct length steps in the computed search directions, pk+1, from the residual vectors rk+1 [Shewchuk, 1994]. Hence, the method converges in at most n iterations to find the true solution. The time complexity of each conjugate gradient iteration scales O(n 2 ) and yields a worst-case running time of O(n 3 ) if all n iterations are used.

Differentiating Conjugate Gradient. The depth of conjugate gradient is dependent on how far away the initial guess is from the true root, but this variability never alters the routine in each iteration. Therefore, we can treat conjugate gradient as a chain of identical processes that can be unrolled one iteration at a time. This description can naturally be extended as treating each conjugate gradient iteration as a single layer in a neural network and connecting them sequentially. Since matrix-vector operations dominate each iteration, we can easily back-propagate through the neural network to compute the gradient ¶L-1 G ¶LG . The red dashed box of Figure 1 represents the partial derivatives of each quantity in a single conjugate gradient iteration with respect to the loss function. We assume access to the partial derivatives ¶ l ¶rk+1 , ¶l ¶ pk+1 , and ¶l ¶xk+1 to compute derivatives for ¶l ¶rk , ¶l ¶ pk , and ¶l ¶xk in the previous iteration. These derivatives are linked together to compute the partial derivative of the loss function with respect to matrix A, ¶l ¶A , where A represents LG. The computation of derivatives between two consecutive layers scales O(n 2 ). Since the gradient through the matrix inverse is neatly represented as a feed-forward neural network, we simply embed the gradient computation into the stochastic gradient descent framework.

5 Experimental Results

In ecological research, landscape-resistance costs are usually combinations of different factors such as roads, slope, elevation, land cover, and range limitations [Cushman et al., 2006; McRae and Beier, 2007]. We generated multiple structured 10×10 resistance grids fromDilkina et al. (2013) to mimic these costs. Then, we computed a single resistance surface by adopting a linear parametrization assigning a specific weight to each factor. These resistance surfaces served as our ground-truth configurations.

5.1 Results for Learning RandomWalks from Partial Trajectories

To comprehensively test our model performance in learning random walks from partial trajectories, we simulated an animal individual for a 1000 time-steps based on a unique linear parametrization of three structured instances. To emulate partial trajectories, we experimented with 4 different scenarios that consisted of 80%, 60%, 40%, and 20% data availability. For each scenario, we executed 5 different variations of CGInv to test the effectiveness of restricting the maximum number of conjugate gradient iterations. We also ran DirectInv, but the runtime was much slower than CGInv, so we omitted graphing the results.

Figure 2 represents the results of different experiments. To observe that training converges for all 5 variations of CGInv, the graph on the left plots the 40% data availability experiment as a representative example of the single-step KL divergence plateauing near 0. We measured the KL-divergence between the learnt and true probabilities of moving from one-step in every location and summed all the individual values into an aggregate measure called the single-step KL divergence. As expected, all measures between the probability transition matrices, hitting time matrices, and weight vectors improve as data availability scales upward. Through the graphs in Figure 2, we verified that our computational framework fitted with CGInv variations can effectively learn the ground-truth values in three different aspects: transition probabilities, hitting times, and the weight vector.

From a computational efficiency standpoint, we noticed that running CGInv with only 2 conjugate gradient iterations produces massive benefits in total run-time while maintaining similar accuracy to other training variations. The benefits of running smaller number of conjugate gradient iterations lie in the fact that one can compute many more gradient updates and the added stochasticity of the gradient updates can help in optimization. Additionally, optimization algorithms are observed to converge faster in terms of total computation if the computed gradients are estimates of the true gradient [Goodfellow et al., 2016]. From this perspective, variations of CGInv with less number of conjugate gradient iterations have the potential to be more efficient and effective in training.

5.2 Mantel Coefficient Optimization Results

To execute Mantel optimization experiments, we first simulated tagging 100 animal individuals randomly initialized across 100 locations expressed as a 10×10 grid. The ground-truth resistance surface was a linear combination of 10 different structured instances. To compute the genetic distance matrix, we programmed a simple genetics simulator that generates Bray-Curtis distances between two animal individuals based on their true commuting time added with Gaussian noise. The magnitude of the Bray-Curtis distances were linearly proportional to the length of the commuting time. Since the Mantel coefficient detects a linear relationship between the two ecological distance matrices, we created the genetic distance matrix based on linear methods so our framework has an opportunity to recover the true resistance grid. To test the efficiency of CGInv, we initialized 5 variations of the training algorithm from CG 2 to CG 100 to recover the original set of weights.

The results of the Mantel optimization experiments are shown in Figure 3. We omitted the results of CG 2 from the graph because it failed to converge to the ground-truth value. However, our computational framework is successful in maximizing the Mantel coefficient for other variations of CGInv as training progresses. An interesting result in the table of Figure 3 is CG 25 obtains a remarkably accurate solution vector remarkably even though the training curve is masked by CG 75 and nearly identical to other CGInv variations. Compared to CG 100, 75, and 50, the output of CG 25 contains noisier solutions due to fewer allowed conjugate gradient iterations. Back-propagating through such solutions might yield more stochastic gradients that allow the training algorithm to converge on more optimal extremum.

Surprisingly, we observed that reducing the number of conjugate gradient iterations seems to have no substantial effect on the efficiency of CGInv. We performed subsequent experiments to measure the average number of conjugate gradient iterations on the initial epochs for CG 100, which is going to be made available in the full version of the paper. The average number of conjugate gradient iterations started at 50 and started dropping by 1 or 2 iterations every couple epochs or so. Due to this observation, we realized that the training progress across different variations of CGInv produced minimal improvements in training time efficiency.

6 Conclusion

Ecological models such as the Circuitscape provide novel and intuitive approaches in evaluating the role of landscape features in genetic differentiation and animal movement. However, the true capacity of these models can be unlocked only when we master the ability to fully utilize ecological data, computational power, and statistical learning. Our computational framework is one such example integrating core concepts from multiple disciplines.

Through the method of moment matching, we showed that our framework successfully learns random walks in two different ecological setups. Wrapping conjugate gradient as layers in a neural network allowed us to scale the efficiency of our approach as evidenced in experiments. Additionally, our computational framework can be readily deployed in real-world ecological studies. As future work, we look forward to embedding the learned random walk model in landscape optimization studies to fuel conservation efforts.

Acknowledgements

Devulapalli acknowledges the support from the Purdue Doctoral Fellowship. This research was supported by NSF grants IIS-1850243, CCF-1918327, CMMI-1763108, CMMI-1935451, as well as by Microsoft AI for Earth grants.