Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/networkx/linalg/algebraicconnectivity.py: 14%

1"""

2Algebraic connectivity and Fiedler vectors of undirected graphs.

3"""

4from functools import partial

6import networkx as nx

7from networkx.utils import (

8 not_implemented_for,

9 np_random_state,

10 reverse_cuthill_mckee_ordering,

11)

13__all__ = [

14 "algebraic_connectivity",

15 "fiedler_vector",

16 "spectral_ordering",

17 "spectral_bisection",

18]

21class _PCGSolver:

22 """Preconditioned conjugate gradient method.

24 To solve Ax = b:

25 M = A.diagonal() # or some other preconditioner

26 solver = _PCGSolver(lambda x: A * x, lambda x: M * x)

27 x = solver.solve(b)

29 The inputs A and M are functions which compute

30 matrix multiplication on the argument.

31 A - multiply by the matrix A in Ax=b

32 M - multiply by M, the preconditioner surrogate for A

34 Warning: There is no limit on number of iterations.

35 """

37 def __init__(self, A, M):

38 self._A = A

39 self._M = M

41 def solve(self, B, tol):

42 import numpy as np

44 # Densifying step - can this be kept sparse?

45 B = np.asarray(B)

46 X = np.ndarray(B.shape, order="F")

47 for j in range(B.shape[1]):

48 X[:, j] = self._solve(B[:, j], tol)

49 return X

51 def _solve(self, b, tol):

52 import numpy as np

53 import scipy as sp

55 A = self._A

56 M = self._M

57 tol *= sp.linalg.blas.dasum(b)

58 # Initialize.

59 x = np.zeros(b.shape)

60 r = b.copy()

61 z = M(r)

62 rz = sp.linalg.blas.ddot(r, z)

63 p = z.copy()

64 # Iterate.

65 while True:

66 Ap = A(p)

67 alpha = rz / sp.linalg.blas.ddot(p, Ap)

68 x = sp.linalg.blas.daxpy(p, x, a=alpha)

69 r = sp.linalg.blas.daxpy(Ap, r, a=-alpha)

70 if sp.linalg.blas.dasum(r) < tol:

71 return x

72 z = M(r)

73 beta = sp.linalg.blas.ddot(r, z)

74 beta, rz = beta / rz, beta

75 p = sp.linalg.blas.daxpy(p, z, a=beta)

78class _LUSolver:

79 """LU factorization.

81 To solve Ax = b:

82 solver = _LUSolver(A)

83 x = solver.solve(b)

85 optional argument `tol` on solve method is ignored but included

86 to match _PCGsolver API.

87 """

89 def __init__(self, A):

90 import scipy as sp

92 self._LU = sp.sparse.linalg.splu(

93 A,

94 permc_spec="MMD_AT_PLUS_A",

95 diag_pivot_thresh=0.0,

96 options={"Equil": True, "SymmetricMode": True},

97 )

99 def solve(self, B, tol=None):

100 import numpy as np

101

102 B = np.asarray(B)

103 X = np.ndarray(B.shape, order="F")

104 for j in range(B.shape[1]):

105 X[:, j] = self._LU.solve(B[:, j])

106 return X

107

108

109def _preprocess_graph(G, weight):

110 """Compute edge weights and eliminate zero-weight edges."""

111 if G.is_directed():

112 H = nx.MultiGraph()

113 H.add_nodes_from(G)

114 H.add_weighted_edges_from(

115 ((u, v, e.get(weight, 1.0)) for u, v, e in G.edges(data=True) if u != v),

116 weight=weight,

117 )

118 G = H

119 if not G.is_multigraph():

120 edges = (

121 (u, v, abs(e.get(weight, 1.0))) for u, v, e in G.edges(data=True) if u != v

122 )

123 else:

124 edges = (

125 (u, v, sum(abs(e.get(weight, 1.0)) for e in G[u][v].values()))

126 for u, v in G.edges()

127 if u != v

128 )

129 H = nx.Graph()

130 H.add_nodes_from(G)

131 H.add_weighted_edges_from((u, v, e) for u, v, e in edges if e != 0)

132 return H

133

134

135def _rcm_estimate(G, nodelist):

136 """Estimate the Fiedler vector using the reverse Cuthill-McKee ordering."""

137 import numpy as np

138

139 G = G.subgraph(nodelist)

140 order = reverse_cuthill_mckee_ordering(G)

141 n = len(nodelist)

142 index = dict(zip(nodelist, range(n)))

143 x = np.ndarray(n, dtype=float)

144 for i, u in enumerate(order):

145 x[index[u]] = i

146 x -= (n - 1) / 2.0

147 return x

148

149

150def _tracemin_fiedler(L, X, normalized, tol, method):

151 """Compute the Fiedler vector of L using the TraceMIN-Fiedler algorithm.

152

153 The Fiedler vector of a connected undirected graph is the eigenvector

154 corresponding to the second smallest eigenvalue of the Laplacian matrix

155 of the graph. This function starts with the Laplacian L, not the Graph.

156

157 Parameters

158 ----------

159 L : Laplacian of a possibly weighted or normalized, but undirected graph

160

161 X : Initial guess for a solution. Usually a matrix of random numbers.

162 This function allows more than one column in X to identify more than

163 one eigenvector if desired.

164

165 normalized : bool

166 Whether the normalized Laplacian matrix is used.

167

168 tol : float

169 Tolerance of relative residual in eigenvalue computation.

170 Warning: There is no limit on number of iterations.

171

172 method : string

173 Should be 'tracemin_pcg' or 'tracemin_lu'.

174 Otherwise exception is raised.

175

176 Returns

177 -------

178 sigma, X : Two NumPy arrays of floats.

179 The lowest eigenvalues and corresponding eigenvectors of L.

180 The size of input X determines the size of these outputs.

181 As this is for Fiedler vectors, the zero eigenvalue (and

182 constant eigenvector) are avoided.

183 """

184 import numpy as np

185 import scipy as sp

186

187 n = X.shape[0]

188

189 if normalized:

190 # Form the normalized Laplacian matrix and determine the eigenvector of

191 # its nullspace.

192 e = np.sqrt(L.diagonal())

193 # TODO: rm csr_array wrapper when spdiags array creation becomes available

194 D = sp.sparse.csr_array(sp.sparse.spdiags(1 / e, 0, n, n, format="csr"))

195 L = D @ L @ D

196 e *= 1.0 / np.linalg.norm(e, 2)

197

198 if normalized:

199

200 def project(X):

201 """Make X orthogonal to the nullspace of L."""

202 X = np.asarray(X)

203 for j in range(X.shape[1]):

204 X[:, j] -= (X[:, j] @ e) * e

205

206 else:

207

208 def project(X):

209 """Make X orthogonal to the nullspace of L."""

210 X = np.asarray(X)

211 for j in range(X.shape[1]):

212 X[:, j] -= X[:, j].sum() / n

213

214 if method == "tracemin_pcg":

215 D = L.diagonal().astype(float)

216 solver = _PCGSolver(lambda x: L @ x, lambda x: D * x)

217 elif method == "tracemin_lu":

218 # Convert A to CSC to suppress SparseEfficiencyWarning.

219 A = sp.sparse.csc_array(L, dtype=float, copy=True)

220 # Force A to be nonsingular. Since A is the Laplacian matrix of a

221 # connected graph, its rank deficiency is one, and thus one diagonal

222 # element needs to modified. Changing to infinity forces a zero in the

223 # corresponding element in the solution.

224 i = (A.indptr[1:] - A.indptr[:-1]).argmax()

225 A[i, i] = float("inf")

226 solver = _LUSolver(A)

227 else:

228 raise nx.NetworkXError(f"Unknown linear system solver: {method}")

229

230 # Initialize.

231 Lnorm = abs(L).sum(axis=1).flatten().max()

232 project(X)

233 W = np.ndarray(X.shape, order="F")

234

235 while True:

236 # Orthonormalize X.

237 X = np.linalg.qr(X)[0]

238 # Compute iteration matrix H.

239 W[:, :] = L @ X

240 H = X.T @ W

241 sigma, Y = sp.linalg.eigh(H, overwrite_a=True)

242 # Compute the Ritz vectors.

243 X = X @ Y

244 # Test for convergence exploiting the fact that L * X == W * Y.

245 res = sp.linalg.blas.dasum(W @ Y[:, 0] - sigma[0] * X[:, 0]) / Lnorm

246 if res < tol:

247 break

248 # Compute X = L \ X / (X' * (L \ X)).

249 # L \ X can have an arbitrary projection on the nullspace of L,

250 # which will be eliminated.

251 W[:, :] = solver.solve(X, tol)

252 X = (sp.linalg.inv(W.T @ X) @ W.T).T # Preserves Fortran storage order.

253 project(X)

254

255 return sigma, np.asarray(X)

256

257

258def _get_fiedler_func(method):

259 """Returns a function that solves the Fiedler eigenvalue problem."""

260 import numpy as np

261

262 if method == "tracemin": # old style keyword <v2.1

263 method = "tracemin_pcg"

264 if method in ("tracemin_pcg", "tracemin_lu"):

265

266 def find_fiedler(L, x, normalized, tol, seed):

267 q = 1 if method == "tracemin_pcg" else min(4, L.shape[0] - 1)

268 X = np.asarray(seed.normal(size=(q, L.shape[0]))).T

269 sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)

270 return sigma[0], X[:, 0]

271

272 elif method == "lanczos" or method == "lobpcg":

273

274 def find_fiedler(L, x, normalized, tol, seed):

275 import scipy as sp

276

277 L = sp.sparse.csc_array(L, dtype=float)

278 n = L.shape[0]

279 if normalized:

280 # TODO: rm csc_array wrapping when spdiags array becomes available

281 D = sp.sparse.csc_array(

282 sp.sparse.spdiags(

283 1.0 / np.sqrt(L.diagonal()), [0], n, n, format="csc"

284 )

285 )

286 L = D @ L @ D

287 if method == "lanczos" or n < 10:

288 # Avoid LOBPCG when n < 10 due to

289 # https://github.com/scipy/scipy/issues/3592

290 # https://github.com/scipy/scipy/pull/3594

291 sigma, X = sp.sparse.linalg.eigsh(

292 L, 2, which="SM", tol=tol, return_eigenvectors=True

293 )

294 return sigma[1], X[:, 1]

295 else:

296 X = np.asarray(np.atleast_2d(x).T)

297 # TODO: rm csr_array wrapping when spdiags array becomes available

298 M = sp.sparse.csr_array(sp.sparse.spdiags(1.0 / L.diagonal(), 0, n, n))

299 Y = np.ones(n)

300 if normalized:

301 Y /= D.diagonal()

302 sigma, X = sp.sparse.linalg.lobpcg(

303 L, X, M=M, Y=np.atleast_2d(Y).T, tol=tol, maxiter=n, largest=False

304 )

305 return sigma[0], X[:, 0]

306

307 else:

308 raise nx.NetworkXError(f"unknown method {method!r}.")

309

310 return find_fiedler

311

312

313@not_implemented_for("directed")

314@np_random_state(5)

315@nx._dispatch(edge_attrs="weight")

316def algebraic_connectivity(

317 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None

318):

319 r"""Returns the algebraic connectivity of an undirected graph.

320

321 The algebraic connectivity of a connected undirected graph is the second

322 smallest eigenvalue of its Laplacian matrix.

323

324 Parameters

325 ----------

326 G : NetworkX graph

327 An undirected graph.

328

329 weight : object, optional (default: None)

330 The data key used to determine the weight of each edge. If None, then

331 each edge has unit weight.

332

333 normalized : bool, optional (default: False)

334 Whether the normalized Laplacian matrix is used.

335

336 tol : float, optional (default: 1e-8)

337 Tolerance of relative residual in eigenvalue computation.

338

339 method : string, optional (default: 'tracemin_pcg')

340 Method of eigenvalue computation. It must be one of the tracemin

341 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)

342 or 'lobpcg' (LOBPCG).

343

344 The TraceMIN algorithm uses a linear system solver. The following

345 values allow specifying the solver to be used.

346

347 =============== ========================================

348 Value Solver

349 =============== ========================================

350 'tracemin_pcg' Preconditioned conjugate gradient method

351 'tracemin_lu' LU factorization

352 =============== ========================================

353

354 seed : integer, random_state, or None (default)

355 Indicator of random number generation state.

356 See :ref:`Randomness<randomness>`.

357

358 Returns

359 -------

360 algebraic_connectivity : float

361 Algebraic connectivity.

362

363 Raises

364 ------

365 NetworkXNotImplemented

366 If G is directed.

367

368 NetworkXError

369 If G has less than two nodes.

370

371 Notes

372 -----

373 Edge weights are interpreted by their absolute values. For MultiGraph's,

374 weights of parallel edges are summed. Zero-weighted edges are ignored.

375

376 See Also

377 --------

378 laplacian_matrix

379

380 Examples

381 --------

382 For undirected graphs algebraic connectivity can tell us if a graph is connected or not

383 `G` is connected iff ``algebraic_connectivity(G) > 0``:

384

385 >>> G = nx.complete_graph(5)

386 >>> nx.algebraic_connectivity(G) > 0

387 True

388 >>> G.add_node(10) # G is no longer connected

389 >>> nx.algebraic_connectivity(G) > 0

390 False

391

392 """

393 if len(G) < 2:

394 raise nx.NetworkXError("graph has less than two nodes.")

395 G = _preprocess_graph(G, weight)

396 if not nx.is_connected(G):

397 return 0.0

398

399 L = nx.laplacian_matrix(G)

400 if L.shape[0] == 2:

401 return 2.0 * L[0, 0] if not normalized else 2.0

402

403 find_fiedler = _get_fiedler_func(method)

404 x = None if method != "lobpcg" else _rcm_estimate(G, G)

405 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)

406 return sigma

407

408

409@not_implemented_for("directed")

410@np_random_state(5)

411@nx._dispatch(edge_attrs="weight")

412def fiedler_vector(

413 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None

414):

415 """Returns the Fiedler vector of a connected undirected graph.

416

417 The Fiedler vector of a connected undirected graph is the eigenvector

418 corresponding to the second smallest eigenvalue of the Laplacian matrix

419 of the graph.

420

421 Parameters

422 ----------

423 G : NetworkX graph

424 An undirected graph.

425

426 weight : object, optional (default: None)

427 The data key used to determine the weight of each edge. If None, then

428 each edge has unit weight.

429

430 normalized : bool, optional (default: False)

431 Whether the normalized Laplacian matrix is used.

432

433 tol : float, optional (default: 1e-8)

434 Tolerance of relative residual in eigenvalue computation.

435

436 method : string, optional (default: 'tracemin_pcg')

437 Method of eigenvalue computation. It must be one of the tracemin

438 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)

439 or 'lobpcg' (LOBPCG).

440

441 The TraceMIN algorithm uses a linear system solver. The following

442 values allow specifying the solver to be used.

443

444 =============== ========================================

445 Value Solver

446 =============== ========================================

447 'tracemin_pcg' Preconditioned conjugate gradient method

448 'tracemin_lu' LU factorization

449 =============== ========================================

450

451 seed : integer, random_state, or None (default)

452 Indicator of random number generation state.

453 See :ref:`Randomness<randomness>`.

454

455 Returns

456 -------

457 fiedler_vector : NumPy array of floats.

458 Fiedler vector.

459

460 Raises

461 ------

462 NetworkXNotImplemented

463 If G is directed.

464

465 NetworkXError

466 If G has less than two nodes or is not connected.

467

468 Notes

469 -----

470 Edge weights are interpreted by their absolute values. For MultiGraph's,

471 weights of parallel edges are summed. Zero-weighted edges are ignored.

472

473 See Also

474 --------

475 laplacian_matrix

476

477 Examples

478 --------

479 Given a connected graph the signs of the values in the Fiedler vector can be

480 used to partition the graph into two components.

481

482 >>> G = nx.barbell_graph(5, 0)

483 >>> nx.fiedler_vector(G, normalized=True, seed=1)

484 array([-0.32864129, -0.32864129, -0.32864129, -0.32864129, -0.26072899,

485 0.26072899, 0.32864129, 0.32864129, 0.32864129, 0.32864129])

486

487 The connected components are the two 5-node cliques of the barbell graph.

488 """

489 import numpy as np

490

491 if len(G) < 2:

492 raise nx.NetworkXError("graph has less than two nodes.")

493 G = _preprocess_graph(G, weight)

494 if not nx.is_connected(G):

495 raise nx.NetworkXError("graph is not connected.")

496

497 if len(G) == 2:

498 return np.array([1.0, -1.0])

499

500 find_fiedler = _get_fiedler_func(method)

501 L = nx.laplacian_matrix(G)

502 x = None if method != "lobpcg" else _rcm_estimate(G, G)

503 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)

504 return fiedler

505

506

507@np_random_state(5)

508@nx._dispatch(edge_attrs="weight")

509def spectral_ordering(

510 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None

511):

512 """Compute the spectral_ordering of a graph.

513

514 The spectral ordering of a graph is an ordering of its nodes where nodes

515 in the same weakly connected components appear contiguous and ordered by

516 their corresponding elements in the Fiedler vector of the component.

517

518 Parameters

519 ----------

520 G : NetworkX graph

521 A graph.

522

523 weight : object, optional (default: None)

524 The data key used to determine the weight of each edge. If None, then

525 each edge has unit weight.

526

527 normalized : bool, optional (default: False)

528 Whether the normalized Laplacian matrix is used.

529

530 tol : float, optional (default: 1e-8)

531 Tolerance of relative residual in eigenvalue computation.

532

533 method : string, optional (default: 'tracemin_pcg')

534 Method of eigenvalue computation. It must be one of the tracemin

535 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)

536 or 'lobpcg' (LOBPCG).

537

538 The TraceMIN algorithm uses a linear system solver. The following

539 values allow specifying the solver to be used.

540

541 =============== ========================================

542 Value Solver

543 =============== ========================================

544 'tracemin_pcg' Preconditioned conjugate gradient method

545 'tracemin_lu' LU factorization

546 =============== ========================================

547

548 seed : integer, random_state, or None (default)

549 Indicator of random number generation state.

550 See :ref:`Randomness<randomness>`.

551

552 Returns

553 -------

554 spectral_ordering : NumPy array of floats.

555 Spectral ordering of nodes.

556

557 Raises

558 ------

559 NetworkXError

560 If G is empty.

561

562 Notes

563 -----

564 Edge weights are interpreted by their absolute values. For MultiGraph's,

565 weights of parallel edges are summed. Zero-weighted edges are ignored.

566

567 See Also

568 --------

569 laplacian_matrix

570 """

571 if len(G) == 0:

572 raise nx.NetworkXError("graph is empty.")

573 G = _preprocess_graph(G, weight)

574

575 find_fiedler = _get_fiedler_func(method)

576 order = []

577 for component in nx.connected_components(G):

578 size = len(component)

579 if size > 2:

580 L = nx.laplacian_matrix(G, component)

581 x = None if method != "lobpcg" else _rcm_estimate(G, component)

582 sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)

583 sort_info = zip(fiedler, range(size), component)

584 order.extend(u for x, c, u in sorted(sort_info))

585 else:

586 order.extend(component)

587

588 return order

589

590

591@nx._dispatch(edge_attrs="weight")

592def spectral_bisection(

593 G, weight="weight", normalized=False, tol=1e-8, method="tracemin_pcg", seed=None

594):

595 """Bisect the graph using the Fiedler vector.

596

597 This method uses the Fiedler vector to bisect a graph.

598 The partition is defined by the nodes which are associated with

599 either positive or negative values in the vector.

600

601 Parameters

602 ----------

603 G : NetworkX Graph

604

605 weight : str, optional (default: weight)

606 The data key used to determine the weight of each edge. If None, then

607 each edge has unit weight.

608

609 normalized : bool, optional (default: False)

610 Whether the normalized Laplacian matrix is used.

611

612 tol : float, optional (default: 1e-8)

613 Tolerance of relative residual in eigenvalue computation.

614

615 method : string, optional (default: 'tracemin_pcg')

616 Method of eigenvalue computation. It must be one of the tracemin

617 options shown below (TraceMIN), 'lanczos' (Lanczos iteration)

618 or 'lobpcg' (LOBPCG).

619

620 The TraceMIN algorithm uses a linear system solver. The following

621 values allow specifying the solver to be used.

622

623 =============== ========================================

624 Value Solver

625 =============== ========================================

626 'tracemin_pcg' Preconditioned conjugate gradient method

627 'tracemin_lu' LU factorization

628 =============== ========================================

629

630 seed : integer, random_state, or None (default)

631 Indicator of random number generation state.

632 See :ref:`Randomness<randomness>`.

633

634 Returns

635 -------

636 bisection : tuple of sets

637 Sets with the bisection of nodes

638

639 Examples

640 --------

641 >>> G = nx.barbell_graph(3, 0)

642 >>> nx.spectral_bisection(G)

643 ({0, 1, 2}, {3, 4, 5})

644

645 References

646 ----------

647 .. [1] M. E. J Newman 'Networks: An Introduction', pages 364-370

648 Oxford University Press 2011.

649 """

650 import numpy as np

651

652 v = nx.fiedler_vector(G, weight, normalized, tol, method, seed)

653 nodes = np.array(list(G))

654 pos_vals = v >= 0

655

656 return set(nodes[~pos_vals]), set(nodes[pos_vals])