Statistics
| Branch: | Revision:

iof-tools / networkxMiCe / networkx-master / networkx / linalg / algebraicconnectivity.py @ 5cef0f13

History | View | Annotate | Download (19.1 KB)

1
# -*- coding: utf-8 -*-
2
# Copyright (C) 2014 ysitu <ysitu@users.noreply.github.com>
3
# All rights reserved.
4
# BSD license.
5
#
6
# Author: ysitu <ysitu@users.noreply.github.com>
7
"""
8
Algebraic connectivity and Fiedler vectors of undirected graphs.
9
"""
10
from functools import partial
11
import networkx as nx
12
from networkx.utils import not_implemented_for
13
from networkx.utils import reverse_cuthill_mckee_ordering
14
from networkx.utils import random_state
15

    
16
try:
17
    from numpy import array, asmatrix, asarray, dot, ndarray, ones, sqrt, zeros
18
    from numpy.linalg import norm, qr
19
    from numpy.random import normal
20
    from scipy.linalg import eigh, inv
21
    from scipy.sparse import csc_matrix, spdiags
22
    from scipy.sparse.linalg import eigsh, lobpcg
23
    __all__ = ['algebraic_connectivity', 'fiedler_vector', 'spectral_ordering']
24
except ImportError:
25
    __all__ = []
26

    
27
try:
28
    from scipy.linalg.blas import dasum, daxpy, ddot
29
except ImportError:
30
    if __all__:
31
        # Make sure the imports succeeded.
32
        # Use minimal replacements if BLAS is unavailable from SciPy.
33
        dasum = partial(norm, ord=1)
34
        ddot = dot
35

    
36
        def daxpy(x, y, a):
37
            y += a * x
38
            return y
39

    
40

    
41
class _PCGSolver(object):
42
    """Preconditioned conjugate gradient method.
43

44
    To solve Ax = b:
45
        M = A.diagonal() # or some other preconditioner
46
        solver = _PCGSolver(lambda x: A * x, lambda x: M * x)
47
        x = solver.solve(b)
48

49
    The inputs A and M are functions which compute
50
    matrix multiplication on the argument.
51
    A - multiply by the matrix A in Ax=b
52
    M - multiply by M, the preconditioner surragate for A
53

54
    Warning: There is no limit on number of iterations.
55
    """
56

    
57
    def __init__(self, A, M):
58
        self._A = A
59
        self._M = M or (lambda x: x.copy())
60

    
61
    def solve(self, B, tol):
62
        B = asarray(B)
63
        X = ndarray(B.shape, order='F')
64
        for j in range(B.shape[1]):
65
            X[:, j] = self._solve(B[:, j], tol)
66
        return X
67

    
68
    def _solve(self, b, tol):
69
        A = self._A
70
        M = self._M
71
        tol *= dasum(b)
72
        # Initialize.
73
        x = zeros(b.shape)
74
        r = b.copy()
75
        z = M(r)
76
        rz = ddot(r, z)
77
        p = z.copy()
78
        # Iterate.
79
        while True:
80
            Ap = A(p)
81
            alpha = rz / ddot(p, Ap)
82
            x = daxpy(p, x, a=alpha)
83
            r = daxpy(Ap, r, a=-alpha)
84
            if dasum(r) < tol:
85
                return x
86
            z = M(r)
87
            beta = ddot(r, z)
88
            beta, rz = beta / rz, beta
89
            p = daxpy(p, z, a=beta)
90

    
91

    
92
class _CholeskySolver(object):
93
    """Cholesky factorization.
94

95
    To solve Ax = b:
96
        solver = _CholeskySolver(A)
97
        x = solver.solve(b)
98

99
    optional argument `tol` on solve method is ignored but included
100
    to match _PCGsolver API.
101
    """
102

    
103
    def __init__(self, A):
104
        if not self._cholesky:
105
            raise nx.NetworkXError('Cholesky solver unavailable.')
106
        self._chol = self._cholesky(A)
107

    
108
    def solve(self, B, tol=None):
109
        return self._chol(B)
110

    
111
    try:
112
        from scikits.sparse.cholmod import cholesky
113
        _cholesky = cholesky
114
    except ImportError:
115
        _cholesky = None
116

    
117

    
118
class _LUSolver(object):
119
    """LU factorization.
120

121
    To solve Ax = b:
122
        solver = _LUSolver(A)
123
        x = solver.solve(b)
124

125
    optional argument `tol` on solve method is ignored but included
126
    to match _PCGsolver API.
127
    """
128

    
129
    def __init__(self, A):
130
        if not self._splu:
131
            raise nx.NetworkXError('LU solver unavailable.')
132
        self._LU = self._splu(A)
133

    
134
    def solve(self, B, tol=None):
135
        B = asarray(B)
136
        X = ndarray(B.shape, order='F')
137
        for j in range(B.shape[1]):
138
            X[:, j] = self._LU.solve(B[:, j])
139
        return X
140

    
141
    try:
142
        from scipy.sparse.linalg import splu
143
        _splu = partial(splu, permc_spec='MMD_AT_PLUS_A', diag_pivot_thresh=0.,
144
                        options={'Equil': True, 'SymmetricMode': True})
145
    except ImportError:
146
        _splu = None
147

    
148

    
149
def _preprocess_graph(G, weight):
150
    """Compute edge weights and eliminate zero-weight edges.
151
    """
152
    if G.is_directed():
153
        H = nx.MultiGraph()
154
        H.add_nodes_from(G)
155
        H.add_weighted_edges_from(((u, v, e.get(weight, 1.))
156
                                   for u, v, e in G.edges(data=True)
157
                                   if u != v), weight=weight)
158
        G = H
159
    if not G.is_multigraph():
160
        edges = ((u, v, abs(e.get(weight, 1.)))
161
                 for u, v, e in G.edges(data=True) if u != v)
162
    else:
163
        edges = ((u, v, sum(abs(e.get(weight, 1.)) for e in G[u][v].values()))
164
                 for u, v in G.edges() if u != v)
165
    H = nx.Graph()
166
    H.add_nodes_from(G)
167
    H.add_weighted_edges_from((u, v, e) for u, v, e in edges if e != 0)
168
    return H
169

    
170

    
171
def _rcm_estimate(G, nodelist):
172
    """Estimate the Fiedler vector using the reverse Cuthill-McKee ordering.
173
    """
174
    G = G.subgraph(nodelist)
175
    order = reverse_cuthill_mckee_ordering(G)
176
    n = len(nodelist)
177
    index = dict(zip(nodelist, range(n)))
178
    x = ndarray(n, dtype=float)
179
    for i, u in enumerate(order):
180
        x[index[u]] = i
181
    x -= (n - 1) / 2.
182
    return x
183

    
184

    
185
def _tracemin_fiedler(L, X, normalized, tol, method):
186
    """Compute the Fiedler vector of L using the TraceMIN-Fiedler algorithm.
187

188
    The Fiedler vector of a connected undirected graph is the eigenvector
189
    corresponding to the second smallest eigenvalue of the Laplacian matrix of
190
    of the graph. This function starts with the Laplacian L, not the Graph.
191

192
    Parameters
193
    ----------
194
    L : Laplacian of a possibly weighted or normalized, but undirected graph
195

196
    X : Initial guess for a solution. Usually a matrix of random numbers.
197
        This function allows more than one column in X to identify more than
198
        one eigenvector if desired.
199

200
    normalized : bool
201
        Whether the normalized Laplacian matrix is used.
202

203
    tol : float
204
        Tolerance of relative residual in eigenvalue computation.
205
        Warning: There is no limit on number of iterations.
206

207
    method : string
208
        Should be 'tracemin_pcg', 'tracemin_chol' or 'tracemin_lu'.
209
        Otherwise exception is raised.
210

211
    Returns
212
    -------
213
    sigma, X : Two NumPy arrays of floats.
214
        The lowest eigenvalues and corresponding eigenvectors of L.
215
        The size of input X determines the size of these outputs.
216
        As this is for Fiedler vectors, the zero eigenvalue (and
217
        constant eigenvector) are avoided.
218
    """
219
    n = X.shape[0]
220

    
221
    if normalized:
222
        # Form the normalized Laplacian matrix and determine the eigenvector of
223
        # its nullspace.
224
        e = sqrt(L.diagonal())
225
        D = spdiags(1. / e, [0], n, n, format='csr')
226
        L = D * L * D
227
        e *= 1. / norm(e, 2)
228

    
229
    if normalized:
230
        def project(X):
231
            """Make X orthogonal to the nullspace of L.
232
            """
233
            X = asarray(X)
234
            for j in range(X.shape[1]):
235
                X[:, j] -= dot(X[:, j], e) * e
236
    else:
237
        def project(X):
238
            """Make X orthogonal to the nullspace of L.
239
            """
240
            X = asarray(X)
241
            for j in range(X.shape[1]):
242
                X[:, j] -= X[:, j].sum() / n
243

    
244
    if method == 'tracemin_pcg':
245
        D = L.diagonal().astype(float)
246
        solver = _PCGSolver(lambda x: L * x, lambda x: D * x)
247
    elif method == 'tracemin_chol' or method == 'tracemin_lu':
248
        # Convert A to CSC to suppress SparseEfficiencyWarning.
249
        A = csc_matrix(L, dtype=float, copy=True)
250
        # Force A to be nonsingular. Since A is the Laplacian matrix of a
251
        # connected graph, its rank deficiency is one, and thus one diagonal
252
        # element needs to modified. Changing to infinity forces a zero in the
253
        # corresponding element in the solution.
254
        i = (A.indptr[1:] - A.indptr[:-1]).argmax()
255
        A[i, i] = float('inf')
256
        if method == 'tracemin_chol':
257
            solver = _CholeskySolver(A)
258
        else:
259
            solver = _LUSolver(A)
260
    else:
261
        raise nx.NetworkXError('Unknown linear system solver: ' + method)
262

    
263
    # Initialize.
264
    Lnorm = abs(L).sum(axis=1).flatten().max()
265
    project(X)
266
    W = asmatrix(ndarray(X.shape, order='F'))
267

    
268
    while True:
269
        # Orthonormalize X.
270
        X = qr(X)[0]
271
        # Compute iteration matrix H.
272
        W[:, :] = L * X
273
        H = X.T * W
274
        sigma, Y = eigh(H, overwrite_a=True)
275
        # Compute the Ritz vectors.
276
        X *= Y
277
        # Test for convergence exploiting the fact that L * X == W * Y.
278
        res = dasum(W * asmatrix(Y)[:, 0] - sigma[0] * X[:, 0]) / Lnorm
279
        if res < tol:
280
            break
281
        # Compute X = L \ X / (X' * (L \ X)).
282
        # L \ X can have an arbitrary projection on the nullspace of L,
283
        # which will be eliminated.
284
        W[:, :] = solver.solve(X, tol)
285
        X = (inv(W.T * X) * W.T).T  # Preserves Fortran storage order.
286
        project(X)
287

    
288
    return sigma, asarray(X)
289

    
290

    
291
def _get_fiedler_func(method):
292
    """Returns a function that solves the Fiedler eigenvalue problem.
293
    """
294
    if method == "tracemin":  # old style keyword <v2.1
295
        method = "tracemin_pcg"
296
    if method in ("tracemin_pcg", "tracemin_chol", "tracemin_lu"):
297
        def find_fiedler(L, x, normalized, tol, seed):
298
            q = 1 if method == 'tracemin_pcg' else min(4, L.shape[0] - 1)
299
            X = asmatrix(seed.normal(size=(q, L.shape[0]))).T
300
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
301
            return sigma[0], X[:, 0]
302
    elif method == 'lanczos' or method == 'lobpcg':
303
        def find_fiedler(L, x, normalized, tol, seed):
304
            L = csc_matrix(L, dtype=float)
305
            n = L.shape[0]
306
            if normalized:
307
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
308
                L = D * L * D
309
            if method == 'lanczos' or n < 10:
310
                # Avoid LOBPCG when n < 10 due to
311
                # https://github.com/scipy/scipy/issues/3592
312
                # https://github.com/scipy/scipy/pull/3594
313
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
314
                                 return_eigenvectors=True)
315
                return sigma[1], X[:, 1]
316
            else:
317
                X = asarray(asmatrix(x).T)
318
                M = spdiags(1. / L.diagonal(), [0], n, n)
319
                Y = ones(n)
320
                if normalized:
321
                    Y /= D.diagonal()
322
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
323
                                  maxiter=n, largest=False)
324
                return sigma[0], X[:, 0]
325
    else:
326
        raise nx.NetworkXError("unknown method '%s'." % method)
327

    
328
    return find_fiedler
329

    
330

    
331
@random_state(5)
332
@not_implemented_for('directed')
333
def algebraic_connectivity(G, weight='weight', normalized=False, tol=1e-8,
334
                           method='tracemin_pcg', seed=None):
335
    """Returns the algebraic connectivity of an undirected graph.
336

337
    The algebraic connectivity of a connected undirected graph is the second
338
    smallest eigenvalue of its Laplacian matrix.
339

340
    Parameters
341
    ----------
342
    G : NetworkX graph
343
        An undirected graph.
344

345
    weight : object, optional (default: None)
346
        The data key used to determine the weight of each edge. If None, then
347
        each edge has unit weight.
348

349
    normalized : bool, optional (default: False)
350
        Whether the normalized Laplacian matrix is used.
351

352
    tol : float, optional (default: 1e-8)
353
        Tolerance of relative residual in eigenvalue computation.
354

355
    method : string, optional (default: 'tracemin_pcg')
356
        Method of eigenvalue computation. It must be one of the tracemin
357
        options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
358
        or 'lobpcg' (LOBPCG).
359

360
        The TraceMIN algorithm uses a linear system solver. The following
361
        values allow specifying the solver to be used.
362

363
        =============== ========================================
364
        Value           Solver
365
        =============== ========================================
366
        'tracemin_pcg'  Preconditioned conjugate gradient method
367
        'tracemin_chol' Cholesky factorization
368
        'tracemin_lu'   LU factorization
369
        =============== ========================================
370

371
    seed : integer, random_state, or None (default)
372
        Indicator of random number generation state.
373
        See :ref:`Randomness<randomness>`.
374

375
    Returns
376
    -------
377
    algebraic_connectivity : float
378
        Algebraic connectivity.
379

380
    Raises
381
    ------
382
    NetworkXNotImplemented
383
        If G is directed.
384

385
    NetworkXError
386
        If G has less than two nodes.
387

388
    Notes
389
    -----
390
    Edge weights are interpreted by their absolute values. For MultiGraph's,
391
    weights of parallel edges are summed. Zero-weighted edges are ignored.
392

393
    To use Cholesky factorization in the TraceMIN algorithm, the
394
    :samp:`scikits.sparse` package must be installed.
395

396
    See Also
397
    --------
398
    laplacian_matrix
399
    """
400
    if len(G) < 2:
401
        raise nx.NetworkXError('graph has less than two nodes.')
402
    G = _preprocess_graph(G, weight)
403
    if not nx.is_connected(G):
404
        return 0.
405

    
406
    L = nx.laplacian_matrix(G)
407
    if L.shape[0] == 2:
408
        return 2. * L[0, 0] if not normalized else 2.
409

    
410
    find_fiedler = _get_fiedler_func(method)
411
    x = None if method != 'lobpcg' else _rcm_estimate(G, G)
412
    sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
413
    return sigma
414

    
415

    
416
@random_state(5)
417
@not_implemented_for('directed')
418
def fiedler_vector(G, weight='weight', normalized=False, tol=1e-8,
419
                   method='tracemin_pcg', seed=None):
420
    """Returns the Fiedler vector of a connected undirected graph.
421

422
    The Fiedler vector of a connected undirected graph is the eigenvector
423
    corresponding to the second smallest eigenvalue of the Laplacian matrix of
424
    of the graph.
425

426
    Parameters
427
    ----------
428
    G : NetworkX graph
429
        An undirected graph.
430

431
    weight : object, optional (default: None)
432
        The data key used to determine the weight of each edge. If None, then
433
        each edge has unit weight.
434

435
    normalized : bool, optional (default: False)
436
        Whether the normalized Laplacian matrix is used.
437

438
    tol : float, optional (default: 1e-8)
439
        Tolerance of relative residual in eigenvalue computation.
440

441
    method : string, optional (default: 'tracemin_pcg')
442
        Method of eigenvalue computation. It must be one of the tracemin
443
        options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
444
        or 'lobpcg' (LOBPCG).
445

446
        The TraceMIN algorithm uses a linear system solver. The following
447
        values allow specifying the solver to be used.
448

449
        =============== ========================================
450
        Value           Solver
451
        =============== ========================================
452
        'tracemin_pcg'  Preconditioned conjugate gradient method
453
        'tracemin_chol' Cholesky factorization
454
        'tracemin_lu'   LU factorization
455
        =============== ========================================
456

457
    seed : integer, random_state, or None (default)
458
        Indicator of random number generation state.
459
        See :ref:`Randomness<randomness>`.
460

461
    Returns
462
    -------
463
    fiedler_vector : NumPy array of floats.
464
        Fiedler vector.
465

466
    Raises
467
    ------
468
    NetworkXNotImplemented
469
        If G is directed.
470

471
    NetworkXError
472
        If G has less than two nodes or is not connected.
473

474
    Notes
475
    -----
476
    Edge weights are interpreted by their absolute values. For MultiGraph's,
477
    weights of parallel edges are summed. Zero-weighted edges are ignored.
478

479
    To use Cholesky factorization in the TraceMIN algorithm, the
480
    :samp:`scikits.sparse` package must be installed.
481

482
    See Also
483
    --------
484
    laplacian_matrix
485
    """
486
    if len(G) < 2:
487
        raise nx.NetworkXError('graph has less than two nodes.')
488
    G = _preprocess_graph(G, weight)
489
    if not nx.is_connected(G):
490
        raise nx.NetworkXError('graph is not connected.')
491

    
492
    if len(G) == 2:
493
        return array([1., -1.])
494

    
495
    find_fiedler = _get_fiedler_func(method)
496
    L = nx.laplacian_matrix(G)
497
    x = None if method != 'lobpcg' else _rcm_estimate(G, G)
498
    sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
499
    return fiedler
500

    
501

    
502
@random_state(5)
503
def spectral_ordering(G, weight='weight', normalized=False, tol=1e-8,
504
                      method='tracemin_pcg', seed=None):
505
    """Compute the spectral_ordering of a graph.
506

507
    The spectral ordering of a graph is an ordering of its nodes where nodes
508
    in the same weakly connected components appear contiguous and ordered by
509
    their corresponding elements in the Fiedler vector of the component.
510

511
    Parameters
512
    ----------
513
    G : NetworkX graph
514
        A graph.
515

516
    weight : object, optional (default: None)
517
        The data key used to determine the weight of each edge. If None, then
518
        each edge has unit weight.
519

520
    normalized : bool, optional (default: False)
521
        Whether the normalized Laplacian matrix is used.
522

523
    tol : float, optional (default: 1e-8)
524
        Tolerance of relative residual in eigenvalue computation.
525

526
    method : string, optional (default: 'tracemin_pcg')
527
        Method of eigenvalue computation. It must be one of the tracemin
528
        options shown below (TraceMIN), 'lanczos' (Lanczos iteration)
529
        or 'lobpcg' (LOBPCG).
530

531
        The TraceMIN algorithm uses a linear system solver. The following
532
        values allow specifying the solver to be used.
533

534
        =============== ========================================
535
        Value           Solver
536
        =============== ========================================
537
        'tracemin_pcg'  Preconditioned conjugate gradient method
538
        'tracemin_chol' Cholesky factorization
539
        'tracemin_lu'   LU factorization
540
        =============== ========================================
541

542
    seed : integer, random_state, or None (default)
543
        Indicator of random number generation state.
544
        See :ref:`Randomness<randomness>`.
545

546
    Returns
547
    -------
548
    spectral_ordering : NumPy array of floats.
549
        Spectral ordering of nodes.
550

551
    Raises
552
    ------
553
    NetworkXError
554
        If G is empty.
555

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

561
    To use Cholesky factorization in the TraceMIN algorithm, the
562
    :samp:`scikits.sparse` package must be installed.
563

564
    See Also
565
    --------
566
    laplacian_matrix
567
    """
568
    if len(G) == 0:
569
        raise nx.NetworkXError('graph is empty.')
570
    G = _preprocess_graph(G, weight)
571

    
572
    find_fiedler = _get_fiedler_func(method)
573
    order = []
574
    for component in nx.connected_components(G):
575
        size = len(component)
576
        if size > 2:
577
            L = nx.laplacian_matrix(G, component)
578
            x = None if method != 'lobpcg' else _rcm_estimate(G, component)
579
            sigma, fiedler = find_fiedler(L, x, normalized, tol, seed)
580
            sort_info = zip(fiedler, range(size), component)
581
            order.extend(u for x, c, u in sorted(sort_info))
582
        else:
583
            order.extend(component)
584

    
585
    return order
586

    
587

    
588
# fixture for nose tests
589
def setup_module(module):
590
    from nose import SkipTest
591
    try:
592
        import numpy
593
        import scipy.sparse
594
    except ImportError:
595
        raise SkipTest('SciPy not available.')