Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Algebraic connectivity and Fiedler vectors of undirected graphs.
3"""
5import networkx as nx
6from networkx.utils import (
7 not_implemented_for,
8 np_random_state,
9 reverse_cuthill_mckee_ordering,
10)
12__all__ = [
13 "algebraic_connectivity",
14 "fiedler_vector",
15 "spectral_ordering",
16 "spectral_bisection",
17]
20class _PCGSolver:
21 """Preconditioned conjugate gradient method.
23 To solve Ax = b:
24 M = A.diagonal() # or some other preconditioner
25 solver = _PCGSolver(lambda x: A * x, lambda x: M * x)
26 x = solver.solve(b)
28 The inputs A and M are functions which compute
29 matrix multiplication on the argument.
30 A - multiply by the matrix A in Ax=b
31 M - multiply by M, the preconditioner surrogate for A
33 Warning: There is no limit on number of iterations.
34 """
36 def __init__(self, A, M):
37 self._A = A
38 self._M = M
40 def solve(self, B, tol):
41 import numpy as np
43 # Densifying step - can this be kept sparse?
44 B = np.asarray(B)
45 X = np.ndarray(B.shape, order="F")
46 for j in range(B.shape[1]):
47 X[:, j] = self._solve(B[:, j], tol)
48 return X
50 def _solve(self, b, tol):
51 import numpy as np
52 import scipy as sp
54 A = self._A
55 M = self._M
56 tol *= sp.linalg.blas.dasum(b)
57 # Initialize.
58 x = np.zeros(b.shape)
59 r = b.copy()
60 z = M(r)
61 rz = sp.linalg.blas.ddot(r, z)
62 p = z.copy()
63 # Iterate.
64 while True:
65 Ap = A(p)
66 alpha = rz / sp.linalg.blas.ddot(p, Ap)
67 x = sp.linalg.blas.daxpy(p, x, a=alpha)
68 r = sp.linalg.blas.daxpy(Ap, r, a=-alpha)
69 if sp.linalg.blas.dasum(r) < tol:
70 return x
71 z = M(r)
72 beta = sp.linalg.blas.ddot(r, z)
73 beta, rz = beta / rz, beta
74 p = sp.linalg.blas.daxpy(p, z, a=beta)
77class _LUSolver:
78 """LU factorization.
80 To solve Ax = b:
81 solver = _LUSolver(A)
82 x = solver.solve(b)
84 optional argument `tol` on solve method is ignored but included
85 to match _PCGsolver API.
86 """
88 def __init__(self, A):
89 import scipy as sp
91 self._LU = sp.sparse.linalg.splu(
92 A,
93 permc_spec="MMD_AT_PLUS_A",
94 diag_pivot_thresh=0.0,
95 options={"Equil": True, "SymmetricMode": True},
96 )
98 def solve(self, B, tol=None):
99 import numpy as np
101 B = np.asarray(B)
102 X = np.ndarray(B.shape, order="F")
103 for j in range(B.shape[1]):
104 X[:, j] = self._LU.solve(B[:, j])
105 return X
108def _preprocess_graph(G, weight):
109 """Compute edge weights and eliminate zero-weight edges."""
110 if G.is_directed():
111 H = nx.MultiGraph()
112 H.add_nodes_from(G)
113 H.add_weighted_edges_from(
114 ((u, v, e.get(weight, 1.0)) for u, v, e in G.edges(data=True) if u != v),
115 weight=weight,
116 )
117 G = H
118 if not G.is_multigraph():
119 edges = (
120 (u, v, abs(e.get(weight, 1.0))) for u, v, e in G.edges(data=True) if u != v
121 )
122 else:
123 edges = (
124 (u, v, sum(abs(e.get(weight, 1.0)) for e in G[u][v].values()))
125 for u, v in G.edges()
126 if u != v
127 )
128 H = nx.Graph()
129 H.add_nodes_from(G)
130 H.add_weighted_edges_from((u, v, e) for u, v, e in edges if e != 0)
131 return H
134def _rcm_estimate(G, nodelist):
135 """Estimate the Fiedler vector using the reverse Cuthill-McKee ordering."""
136 import numpy as np
138 G = G.subgraph(nodelist)
139 order = reverse_cuthill_mckee_ordering(G)
140 n = len(nodelist)
141 index = dict(zip(nodelist, range(n)))
142 x = np.ndarray(n, dtype=float)
143 for i, u in enumerate(order):
144 x[index[u]] = i
145 x -= (n - 1) / 2.0
146 return x
149def _tracemin_fiedler(L, X, normalized, tol, method):
150 """Compute the Fiedler vector of L using the TraceMIN-Fiedler algorithm.
152 The Fiedler vector of a connected undirected graph is the eigenvector
153 corresponding to the second smallest eigenvalue of the Laplacian matrix
154 of the graph. This function starts with the Laplacian L, not the Graph.
156 Parameters
157 ----------
158 L : Laplacian of a possibly weighted or normalized, but undirected graph
160 X : Initial guess for a solution. Usually a matrix of random numbers.
161 This function allows more than one column in X to identify more than
162 one eigenvector if desired.
164 normalized : bool
165 Whether the normalized Laplacian matrix is used.
167 tol : float
168 Tolerance of relative residual in eigenvalue computation.
169 Warning: There is no limit on number of iterations.
171 method : string
172 Should be 'tracemin_pcg' or 'tracemin_lu'.
173 Otherwise exception is raised.
175 Returns
176 -------
177 sigma, X : Two NumPy arrays of floats.
178 The lowest eigenvalues and corresponding eigenvectors of L.
179 The size of input X determines the size of these outputs.
180 As this is for Fiedler vectors, the zero eigenvalue (and
181 constant eigenvector) are avoided.
182 """
183 import numpy as np
184 import scipy as sp
186 n = X.shape[0]
188 if normalized:
189 # Form the normalized Laplacian matrix and determine the eigenvector of
190 # its nullspace.
191 e = np.sqrt(L.diagonal())
192 D = sp.sparse.dia_array((1 / e, 0), shape=(n, n)).tocsr()
193 L = D @ L @ D
194 e *= 1.0 / np.linalg.norm(e, 2)
196 if normalized:
198 def project(X):
199 """Make X orthogonal to the nullspace of L."""
200 X = np.asarray(X)
201 for j in range(X.shape[1]):
202 X[:, j] -= (X[:, j] @ e) * e
204 else:
206 def project(X):
207 """Make X orthogonal to the nullspace of L."""
208 X = np.asarray(X)
209 for j in range(X.shape[1]):
210 X[:, j] -= X[:, j].sum() / n
212 if method == "tracemin_pcg":
213 D = L.diagonal().astype(float)
214 solver = _PCGSolver(lambda x: L @ x, lambda x: D * x)
215 elif method == "tracemin_lu":
216 # Convert A to CSC to suppress SparseEfficiencyWarning.
217 A = sp.sparse.csc_array(L, dtype=float, copy=True)
218 # Force A to be nonsingular. Since A is the Laplacian matrix of a
219 # connected graph, its rank deficiency is one, and thus one diagonal
220 # element needs to modified. Changing to infinity forces a zero in the
221 # corresponding element in the solution.
222 i = (A.indptr[1:] - A.indptr[:-1]).argmax()
223 A[i, i] = np.inf
224 solver = _LUSolver(A)
225 else:
226 raise nx.NetworkXError(f"Unknown linear system solver: {method}")
228 # Initialize.
229 Lnorm = abs(L).sum(axis=1).flatten().max()
230 project(X)
231 W = np.ndarray(X.shape, order="F")
233 while True:
234 # Orthonormalize X.
235 X = np.linalg.qr(X)[0]
236 # Compute iteration matrix H.
237 W[:, :] = L @ X
238 H = X.T @ W
239 sigma, Y = sp.linalg.eigh(H, overwrite_a=True)
240 # Compute the Ritz vectors.
241 X = X @ Y
242 # Test for convergence exploiting the fact that L * X == W * Y.
243 res = sp.linalg.blas.dasum(W @ Y[:, 0] - sigma[0] * X[:, 0]) / Lnorm
244 if res < tol:
245 break
246 # Compute X = L \ X / (X' * (L \ X)).
247 # L \ X can have an arbitrary projection on the nullspace of L,
248 # which will be eliminated.
249 W[:, :] = solver.solve(X, tol)
250 X = (sp.linalg.inv(W.T @ X) @ W.T).T # Preserves Fortran storage order.
251 project(X)
253 return sigma, np.asarray(X)
256def _get_fiedler_func(method):
257 """Returns a function that solves the Fiedler eigenvalue problem."""
258 import numpy as np
260 if method == "tracemin": # old style keyword <v2.1
261 method = "tracemin_pcg"
262 if method in ("tracemin_pcg", "tracemin_lu"):
264 def find_fiedler(L, x, normalized, tol, seed):
265 q = 1 if method == "tracemin_pcg" else min(4, L.shape[0] - 1)
266 X = np.asarray(seed.normal(size=(q, L.shape[0]))).T
267 sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
268 return sigma[0], X[:, 0]
270 elif method == "lanczos" or method == "lobpcg":
272 def find_fiedler(L, x, normalized, tol, seed):
273 import scipy as sp
275 L = sp.sparse.csc_array(L, dtype=float)
276 n = L.shape[0]
277 if normalized:
278 D = sp.sparse.dia_array(
279 (1.0 / np.sqrt(L.diagonal()), 0), shape=(n, n)
280 ).tocsc()
281 L = D @ L @ D
282 if method == "lanczos" or n < 10:
283 # Avoid LOBPCG when n < 10 due to
284 # https://github.com/scipy/scipy/issues/3592
285 # https://github.com/scipy/scipy/pull/3594
286 sigma, X = sp.sparse.linalg.eigsh(
287 L, 2, which="SM", tol=tol, return_eigenvectors=True
288 )
289 return sigma[1], X[:, 1]
290 else:
291 X = np.asarray(np.atleast_2d(x).T)
292 M = sp.sparse.dia_array((1.0 / L.diagonal(), 0), shape=(n, n)).tocsr()
293 Y = np.ones(n)
294 if normalized:
295 Y /= D.diagonal()
296 sigma, X = sp.sparse.linalg.lobpcg(
297 L, X, M=M, Y=np.atleast_2d(Y).T, tol=tol, maxiter=n, largest=False
298 )
299 return sigma[0], X[:, 0]
301 else:
302 raise nx.NetworkXError(f"unknown method {method!r}.")
304 return find_fiedler
307@not_implemented_for("directed")
308@np_random_state(5)
309@nx._dispatchable(edge_attrs="weight")
310def algebraic_connectivity(
311 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None
312):
313 r"""Returns the algebraic connectivity of an undirected graph.
315 The algebraic connectivity of a connected undirected graph is the second
316 smallest eigenvalue of its Laplacian matrix.
318 Parameters
319 ----------
320 G : NetworkX graph
321 An undirected graph.
323 weight : object, optional (default: None)
324 The data key used to determine the weight of each edge. If None, then
325 each edge has unit weight.
327 normalized : bool, optional (default: False)
328 Whether the normalized Laplacian matrix is used.
330 tol : float, optional (default: 1e-8)
331 Tolerance of relative residual in eigenvalue computation.
333 method : string, optional (default: 'tracemin_pcg')
334 Method of eigenvalue computation. It must be one of the tracemin
335 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
336 or 'lobpcg' (LOBPCG).
338 The TraceMIN algorithm uses a linear system solver. The following
339 values allow specifying the solver to be used.
341 =============== ========================================
342 Value Solver
343 =============== ========================================
344 'tracemin_pcg' Preconditioned conjugate gradient method
345 'tracemin_lu' LU factorization
346 =============== ========================================
348 seed : integer, random_state, or None (default)
349 Indicator of random number generation state.
350 See :ref:`Randomness<randomness>`.
352 Returns
353 -------
354 algebraic_connectivity : float
355 Algebraic connectivity.
357 Raises
358 ------
359 NetworkXNotImplemented
360 If G is directed.
362 NetworkXError
363 If G has less than two nodes.
365 Notes
366 -----
367 Edge weights are interpreted by their absolute values. For MultiGraph's,
368 weights of parallel edges are summed. Zero-weighted edges are ignored.
370 See Also
371 --------
372 laplacian_matrix
374 Examples
375 --------
376 For undirected graphs algebraic connectivity can tell us if a graph is connected or not
377 `G` is connected iff ``algebraic_connectivity(G) > 0``:
379 >>> G = nx.complete_graph(5)
380 >>> nx.algebraic_connectivity(G) > 0
381 True
382 >>> G.add_node(10) # G is no longer connected
383 >>> nx.algebraic_connectivity(G) > 0
384 False
386 """
387 if len(G) < 2:
388 raise nx.NetworkXError("graph has less than two nodes.")
389 G = _preprocess_graph(G, weight)
390 if not nx.is_connected(G):
391 return 0.0
393 L = nx.laplacian_matrix(G)
394 if L.shape[0] == 2:
395 return 2.0 * float(L[0, 0]) if not normalized else 2.0
397 find_fiedler = _get_fiedler_func(method)
398 x = None if method != "lobpcg" else _rcm_estimate(G, G)
399 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
400 return float(sigma)
403@not_implemented_for("directed")
404@np_random_state(5)
405@nx._dispatchable(edge_attrs="weight")
406def fiedler_vector(
407 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None
408):
409 """Returns the Fiedler vector of a connected undirected graph.
411 The Fiedler vector of a connected undirected graph is the eigenvector
412 corresponding to the second smallest eigenvalue of the Laplacian matrix
413 of the graph.
415 Parameters
416 ----------
417 G : NetworkX graph
418 An undirected graph.
420 weight : object, optional (default: None)
421 The data key used to determine the weight of each edge. If None, then
422 each edge has unit weight.
424 normalized : bool, optional (default: False)
425 Whether the normalized Laplacian matrix is used.
427 tol : float, optional (default: 1e-8)
428 Tolerance of relative residual in eigenvalue computation.
430 method : string, optional (default: 'tracemin_pcg')
431 Method of eigenvalue computation. It must be one of the tracemin
432 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
433 or 'lobpcg' (LOBPCG).
435 The TraceMIN algorithm uses a linear system solver. The following
436 values allow specifying the solver to be used.
438 =============== ========================================
439 Value Solver
440 =============== ========================================
441 'tracemin_pcg' Preconditioned conjugate gradient method
442 'tracemin_lu' LU factorization
443 =============== ========================================
445 seed : integer, random_state, or None (default)
446 Indicator of random number generation state.
447 See :ref:`Randomness<randomness>`.
449 Returns
450 -------
451 fiedler_vector : NumPy array of floats.
452 Fiedler vector.
454 Raises
455 ------
456 NetworkXNotImplemented
457 If G is directed.
459 NetworkXError
460 If G has less than two nodes or is not connected.
462 Notes
463 -----
464 Edge weights are interpreted by their absolute values. For MultiGraph's,
465 weights of parallel edges are summed. Zero-weighted edges are ignored.
467 See Also
468 --------
469 laplacian_matrix
471 Examples
472 --------
473 Given a connected graph the signs of the values in the Fiedler vector can be
474 used to partition the graph into two components.
476 >>> G = nx.barbell_graph(5, 0)
477 >>> nx.fiedler_vector(G, normalized=True, seed=1)
478 array([-0.32864129, -0.32864129, -0.32864129, -0.32864129, -0.26072899,
479 0.26072899, 0.32864129, 0.32864129, 0.32864129, 0.32864129])
481 The connected components are the two 5-node cliques of the barbell graph.
482 """
483 import numpy as np
485 if len(G) < 2:
486 raise nx.NetworkXError("graph has less than two nodes.")
487 G = _preprocess_graph(G, weight)
488 if not nx.is_connected(G):
489 raise nx.NetworkXError("graph is not connected.")
491 if len(G) == 2:
492 return np.array([1.0, -1.0])
494 find_fiedler = _get_fiedler_func(method)
495 L = nx.laplacian_matrix(G)
496 x = None if method != "lobpcg" else _rcm_estimate(G, G)
497 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
498 return fiedler
501@np_random_state(5)
502@nx._dispatchable(edge_attrs="weight")
503def spectral_ordering(
504 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None
505):
506 """Compute the spectral_ordering of a graph.
508 The spectral ordering of a graph is an ordering of its nodes where nodes
509 in the same weakly connected components appear contiguous and ordered by
510 their corresponding elements in the Fiedler vector of the component.
512 Parameters
513 ----------
514 G : NetworkX graph
515 A graph.
517 weight : object, optional (default: None)
518 The data key used to determine the weight of each edge. If None, then
519 each edge has unit weight.
521 normalized : bool, optional (default: False)
522 Whether the normalized Laplacian matrix is used.
524 tol : float, optional (default: 1e-8)
525 Tolerance of relative residual in eigenvalue computation.
527 method : string, optional (default: 'tracemin_pcg')
528 Method of eigenvalue computation. It must be one of the tracemin
529 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
530 or 'lobpcg' (LOBPCG).
532 The TraceMIN algorithm uses a linear system solver. The following
533 values allow specifying the solver to be used.
535 =============== ========================================
536 Value Solver
537 =============== ========================================
538 'tracemin_pcg' Preconditioned conjugate gradient method
539 'tracemin_lu' LU factorization
540 =============== ========================================
542 seed : integer, random_state, or None (default)
543 Indicator of random number generation state.
544 See :ref:`Randomness<randomness>`.
546 Returns
547 -------
548 spectral_ordering : NumPy array of floats.
549 Spectral ordering of nodes.
551 Raises
552 ------
553 NetworkXError
554 If G is empty.
556 Notes
557 -----
558 Edge weights are interpreted by their absolute values. For MultiGraph's,
559 weights of parallel edges are summed. Zero-weighted edges are ignored.
561 See Also
562 --------
563 laplacian_matrix
564 """
565 if len(G) == 0:
566 raise nx.NetworkXError("graph is empty.")
567 G = _preprocess_graph(G, weight)
569 find_fiedler = _get_fiedler_func(method)
570 order = []
571 for component in nx.connected_components(G):
572 size = len(component)
573 if size > 2:
574 L = nx.laplacian_matrix(G, component)
575 x = None if method != "lobpcg" else _rcm_estimate(G, component)
576 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
577 sort_info = zip(fiedler, range(size), component)
578 order.extend(u for x, c, u in sorted(sort_info))
579 else:
580 order.extend(component)
582 return order
585@nx._dispatchable(edge_attrs="weight")
586def spectral_bisection(
587 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None
588):
589 """Bisect the graph using the Fiedler vector.
591 This method uses the Fiedler vector to bisect a graph.
592 The partition is defined by the nodes which are associated with
593 either positive or negative values in the vector.
595 Parameters
596 ----------
597 G : NetworkX Graph
599 weight : str, optional (default: weight)
600 The data key used to determine the weight of each edge. If None, then
601 each edge has unit weight.
603 normalized : bool, optional (default: False)
604 Whether the normalized Laplacian matrix is used.
606 tol : float, optional (default: 1e-8)
607 Tolerance of relative residual in eigenvalue computation.
609 method : string, optional (default: 'tracemin_pcg')
610 Method of eigenvalue computation. It must be one of the tracemin
611 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
612 or 'lobpcg' (LOBPCG).
614 The TraceMIN algorithm uses a linear system solver. The following
615 values allow specifying the solver to be used.
617 =============== ========================================
618 Value Solver
619 =============== ========================================
620 'tracemin_pcg' Preconditioned conjugate gradient method
621 'tracemin_lu' LU factorization
622 =============== ========================================
624 seed : integer, random_state, or None (default)
625 Indicator of random number generation state.
626 See :ref:`Randomness<randomness>`.
628 Returns
629 -------
630 bisection : tuple of sets
631 Sets with the bisection of nodes
633 Examples
634 --------
635 >>> G = nx.barbell_graph(3, 0)
636 >>> nx.spectral_bisection(G)
637 ({0, 1, 2}, {3, 4, 5})
639 References
640 ----------
641 .. [1] M. E. J Newman 'Networks: An Introduction', pages 364-370
642 Oxford University Press 2011.
643 """
644 import numpy as np
646 v = nx.fiedler_vector(G, weight, normalized, tol, method, seed)
647 nodes = np.array(list(G))
648 pos_vals = v >= 0
650 return set(nodes[~pos_vals].tolist()), set(nodes[pos_vals].tolist())