Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/csgraph/_laplacian.py: 10%
153 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-03-22 06:44 +0000
1"""
2Laplacian of a compressed-sparse graph
3"""
5import numpy as np
6from scipy.sparse import issparse
7from scipy.sparse.linalg import LinearOperator
8from scipy.sparse._sputils import convert_pydata_sparse_to_scipy, is_pydata_spmatrix
11###############################################################################
12# Graph laplacian
13def laplacian(
14 csgraph,
15 normed=False,
16 return_diag=False,
17 use_out_degree=False,
18 *,
19 copy=True,
20 form="array",
21 dtype=None,
22 symmetrized=False,
23):
24 """
25 Return the Laplacian of a directed graph.
27 Parameters
28 ----------
29 csgraph : array_like or sparse matrix, 2 dimensions
30 compressed-sparse graph, with shape (N, N).
31 normed : bool, optional
32 If True, then compute symmetrically normalized Laplacian.
33 Default: False.
34 return_diag : bool, optional
35 If True, then also return an array related to vertex degrees.
36 Default: False.
37 use_out_degree : bool, optional
38 If True, then use out-degree instead of in-degree.
39 This distinction matters only if the graph is asymmetric.
40 Default: False.
41 copy: bool, optional
42 If False, then change `csgraph` in place if possible,
43 avoiding doubling the memory use.
44 Default: True, for backward compatibility.
45 form: 'array', or 'function', or 'lo'
46 Determines the format of the output Laplacian:
48 * 'array' is a numpy array;
49 * 'function' is a pointer to evaluating the Laplacian-vector
50 or Laplacian-matrix product;
51 * 'lo' results in the format of the `LinearOperator`.
53 Choosing 'function' or 'lo' always avoids doubling
54 the memory use, ignoring `copy` value.
55 Default: 'array', for backward compatibility.
56 dtype: None or one of numeric numpy dtypes, optional
57 The dtype of the output. If ``dtype=None``, the dtype of the
58 output matches the dtype of the input csgraph, except for
59 the case ``normed=True`` and integer-like csgraph, where
60 the output dtype is 'float' allowing accurate normalization,
61 but dramatically increasing the memory use.
62 Default: None, for backward compatibility.
63 symmetrized: bool, optional
64 If True, then the output Laplacian is symmetric/Hermitian.
65 The symmetrization is done by ``csgraph + csgraph.T.conj``
66 without dividing by 2 to preserve integer dtypes if possible
67 prior to the construction of the Laplacian.
68 The symmetrization will increase the memory footprint of
69 sparse matrices unless the sparsity pattern is symmetric or
70 `form` is 'function' or 'lo'.
71 Default: False, for backward compatibility.
73 Returns
74 -------
75 lap : ndarray, or sparse matrix, or `LinearOperator`
76 The N x N Laplacian of csgraph. It will be a NumPy array (dense)
77 if the input was dense, or a sparse matrix otherwise, or
78 the format of a function or `LinearOperator` if
79 `form` equals 'function' or 'lo', respectively.
80 diag : ndarray, optional
81 The length-N main diagonal of the Laplacian matrix.
82 For the normalized Laplacian, this is the array of square roots
83 of vertex degrees or 1 if the degree is zero.
85 Notes
86 -----
87 The Laplacian matrix of a graph is sometimes referred to as the
88 "Kirchhoff matrix" or just the "Laplacian", and is useful in many
89 parts of spectral graph theory.
90 In particular, the eigen-decomposition of the Laplacian can give
91 insight into many properties of the graph, e.g.,
92 is commonly used for spectral data embedding and clustering.
94 The constructed Laplacian doubles the memory use if ``copy=True`` and
95 ``form="array"`` which is the default.
96 Choosing ``copy=False`` has no effect unless ``form="array"``
97 or the matrix is sparse in the ``coo`` format, or dense array, except
98 for the integer input with ``normed=True`` that forces the float output.
100 Sparse input is reformatted into ``coo`` if ``form="array"``,
101 which is the default.
103 If the input adjacency matrix is not symmetric, the Laplacian is
104 also non-symmetric unless ``symmetrized=True`` is used.
106 Diagonal entries of the input adjacency matrix are ignored and
107 replaced with zeros for the purpose of normalization where ``normed=True``.
108 The normalization uses the inverse square roots of row-sums of the input
109 adjacency matrix, and thus may fail if the row-sums contain
110 negative or complex with a non-zero imaginary part values.
112 The normalization is symmetric, making the normalized Laplacian also
113 symmetric if the input csgraph was symmetric.
115 References
116 ----------
117 .. [1] Laplacian matrix. https://en.wikipedia.org/wiki/Laplacian_matrix
119 Examples
120 --------
121 >>> import numpy as np
122 >>> from scipy.sparse import csgraph
124 Our first illustration is the symmetric graph
126 >>> G = np.arange(4) * np.arange(4)[:, np.newaxis]
127 >>> G
128 array([[0, 0, 0, 0],
129 [0, 1, 2, 3],
130 [0, 2, 4, 6],
131 [0, 3, 6, 9]])
133 and its symmetric Laplacian matrix
135 >>> csgraph.laplacian(G)
136 array([[ 0, 0, 0, 0],
137 [ 0, 5, -2, -3],
138 [ 0, -2, 8, -6],
139 [ 0, -3, -6, 9]])
141 The non-symmetric graph
143 >>> G = np.arange(9).reshape(3, 3)
144 >>> G
145 array([[0, 1, 2],
146 [3, 4, 5],
147 [6, 7, 8]])
149 has different row- and column sums, resulting in two varieties
150 of the Laplacian matrix, using an in-degree, which is the default
152 >>> L_in_degree = csgraph.laplacian(G)
153 >>> L_in_degree
154 array([[ 9, -1, -2],
155 [-3, 8, -5],
156 [-6, -7, 7]])
158 or alternatively an out-degree
160 >>> L_out_degree = csgraph.laplacian(G, use_out_degree=True)
161 >>> L_out_degree
162 array([[ 3, -1, -2],
163 [-3, 8, -5],
164 [-6, -7, 13]])
166 Constructing a symmetric Laplacian matrix, one can add the two as
168 >>> L_in_degree + L_out_degree.T
169 array([[ 12, -4, -8],
170 [ -4, 16, -12],
171 [ -8, -12, 20]])
173 or use the ``symmetrized=True`` option
175 >>> csgraph.laplacian(G, symmetrized=True)
176 array([[ 12, -4, -8],
177 [ -4, 16, -12],
178 [ -8, -12, 20]])
180 that is equivalent to symmetrizing the original graph
182 >>> csgraph.laplacian(G + G.T)
183 array([[ 12, -4, -8],
184 [ -4, 16, -12],
185 [ -8, -12, 20]])
187 The goal of normalization is to make the non-zero diagonal entries
188 of the Laplacian matrix to be all unit, also scaling off-diagonal
189 entries correspondingly. The normalization can be done manually, e.g.,
191 >>> G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])
192 >>> L, d = csgraph.laplacian(G, return_diag=True)
193 >>> L
194 array([[ 2, -1, -1],
195 [-1, 2, -1],
196 [-1, -1, 2]])
197 >>> d
198 array([2, 2, 2])
199 >>> scaling = np.sqrt(d)
200 >>> scaling
201 array([1.41421356, 1.41421356, 1.41421356])
202 >>> (1/scaling)*L*(1/scaling)
203 array([[ 1. , -0.5, -0.5],
204 [-0.5, 1. , -0.5],
205 [-0.5, -0.5, 1. ]])
207 Or using ``normed=True`` option
209 >>> L, d = csgraph.laplacian(G, return_diag=True, normed=True)
210 >>> L
211 array([[ 1. , -0.5, -0.5],
212 [-0.5, 1. , -0.5],
213 [-0.5, -0.5, 1. ]])
215 which now instead of the diagonal returns the scaling coefficients
217 >>> d
218 array([1.41421356, 1.41421356, 1.41421356])
220 Zero scaling coefficients are substituted with 1s, where scaling
221 has thus no effect, e.g.,
223 >>> G = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]])
224 >>> G
225 array([[0, 0, 0],
226 [0, 0, 1],
227 [0, 1, 0]])
228 >>> L, d = csgraph.laplacian(G, return_diag=True, normed=True)
229 >>> L
230 array([[ 0., -0., -0.],
231 [-0., 1., -1.],
232 [-0., -1., 1.]])
233 >>> d
234 array([1., 1., 1.])
236 Only the symmetric normalization is implemented, resulting
237 in a symmetric Laplacian matrix if and only if its graph is symmetric
238 and has all non-negative degrees, like in the examples above.
240 The output Laplacian matrix is by default a dense array or a sparse matrix
241 inferring its shape, format, and dtype from the input graph matrix:
243 >>> G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]).astype(np.float32)
244 >>> G
245 array([[0., 1., 1.],
246 [1., 0., 1.],
247 [1., 1., 0.]], dtype=float32)
248 >>> csgraph.laplacian(G)
249 array([[ 2., -1., -1.],
250 [-1., 2., -1.],
251 [-1., -1., 2.]], dtype=float32)
253 but can alternatively be generated matrix-free as a LinearOperator:
255 >>> L = csgraph.laplacian(G, form="lo")
256 >>> L
257 <3x3 _CustomLinearOperator with dtype=float32>
258 >>> L(np.eye(3))
259 array([[ 2., -1., -1.],
260 [-1., 2., -1.],
261 [-1., -1., 2.]])
263 or as a lambda-function:
265 >>> L = csgraph.laplacian(G, form="function")
266 >>> L
267 <function _laplace.<locals>.<lambda> at 0x0000012AE6F5A598>
268 >>> L(np.eye(3))
269 array([[ 2., -1., -1.],
270 [-1., 2., -1.],
271 [-1., -1., 2.]])
273 The Laplacian matrix is used for
274 spectral data clustering and embedding
275 as well as for spectral graph partitioning.
276 Our final example illustrates the latter
277 for a noisy directed linear graph.
279 >>> from scipy.sparse import diags, random
280 >>> from scipy.sparse.linalg import lobpcg
282 Create a directed linear graph with ``N=35`` vertices
283 using a sparse adjacency matrix ``G``:
285 >>> N = 35
286 >>> G = diags(np.ones(N-1), 1, format="csr")
288 Fix a random seed ``rng`` and add a random sparse noise to the graph ``G``:
290 >>> rng = np.random.default_rng()
291 >>> G += 1e-2 * random(N, N, density=0.1, random_state=rng)
293 Set initial approximations for eigenvectors:
295 >>> X = rng.random((N, 2))
297 The constant vector of ones is always a trivial eigenvector
298 of the non-normalized Laplacian to be filtered out:
300 >>> Y = np.ones((N, 1))
302 Alternating (1) the sign of the graph weights allows determining
303 labels for spectral max- and min- cuts in a single loop.
304 Since the graph is undirected, the option ``symmetrized=True``
305 must be used in the construction of the Laplacian.
306 The option ``normed=True`` cannot be used in (2) for the negative weights
307 here as the symmetric normalization evaluates square roots.
308 The option ``form="lo"`` in (2) is matrix-free, i.e., guarantees
309 a fixed memory footprint and read-only access to the graph.
310 Calling the eigenvalue solver ``lobpcg`` (3) computes the Fiedler vector
311 that determines the labels as the signs of its components in (5).
312 Since the sign in an eigenvector is not deterministic and can flip,
313 we fix the sign of the first component to be always +1 in (4).
315 >>> for cut in ["max", "min"]:
316 ... G = -G # 1.
317 ... L = csgraph.laplacian(G, symmetrized=True, form="lo") # 2.
318 ... _, eves = lobpcg(L, X, Y=Y, largest=False, tol=1e-3) # 3.
319 ... eves *= np.sign(eves[0, 0]) # 4.
320 ... print(cut + "-cut labels:\\n", 1 * (eves[:, 0]>0)) # 5.
321 max-cut labels:
322 [1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]
323 min-cut labels:
324 [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
326 As anticipated for a (slightly noisy) linear graph,
327 the max-cut strips all the edges of the graph coloring all
328 odd vertices into one color and all even vertices into another one,
329 while the balanced min-cut partitions the graph
330 in the middle by deleting a single edge.
331 Both determined partitions are optimal.
332 """
333 is_pydata_sparse = is_pydata_spmatrix(csgraph)
334 if is_pydata_sparse:
335 pydata_sparse_cls = csgraph.__class__
336 csgraph = convert_pydata_sparse_to_scipy(csgraph)
337 if csgraph.ndim != 2 or csgraph.shape[0] != csgraph.shape[1]:
338 raise ValueError('csgraph must be a square matrix or array')
340 if normed and (
341 np.issubdtype(csgraph.dtype, np.signedinteger)
342 or np.issubdtype(csgraph.dtype, np.uint)
343 ):
344 csgraph = csgraph.astype(np.float64)
346 if form == "array":
347 create_lap = (
348 _laplacian_sparse if issparse(csgraph) else _laplacian_dense
349 )
350 else:
351 create_lap = (
352 _laplacian_sparse_flo
353 if issparse(csgraph)
354 else _laplacian_dense_flo
355 )
357 degree_axis = 1 if use_out_degree else 0
359 lap, d = create_lap(
360 csgraph,
361 normed=normed,
362 axis=degree_axis,
363 copy=copy,
364 form=form,
365 dtype=dtype,
366 symmetrized=symmetrized,
367 )
368 if is_pydata_sparse:
369 lap = pydata_sparse_cls.from_scipy_sparse(lap)
370 if return_diag:
371 return lap, d
372 return lap
375def _setdiag_dense(m, d):
376 step = len(d) + 1
377 m.flat[::step] = d
380def _laplace(m, d):
381 return lambda v: v * d[:, np.newaxis] - m @ v
384def _laplace_normed(m, d, nd):
385 laplace = _laplace(m, d)
386 return lambda v: nd[:, np.newaxis] * laplace(v * nd[:, np.newaxis])
389def _laplace_sym(m, d):
390 return (
391 lambda v: v * d[:, np.newaxis]
392 - m @ v
393 - np.transpose(np.conjugate(np.transpose(np.conjugate(v)) @ m))
394 )
397def _laplace_normed_sym(m, d, nd):
398 laplace_sym = _laplace_sym(m, d)
399 return lambda v: nd[:, np.newaxis] * laplace_sym(v * nd[:, np.newaxis])
402def _linearoperator(mv, shape, dtype):
403 return LinearOperator(matvec=mv, matmat=mv, shape=shape, dtype=dtype)
406def _laplacian_sparse_flo(graph, normed, axis, copy, form, dtype, symmetrized):
407 # The keyword argument `copy` is unused and has no effect here.
408 del copy
410 if dtype is None:
411 dtype = graph.dtype
413 graph_sum = np.asarray(graph.sum(axis=axis)).ravel()
414 graph_diagonal = graph.diagonal()
415 diag = graph_sum - graph_diagonal
416 if symmetrized:
417 graph_sum += np.asarray(graph.sum(axis=1 - axis)).ravel()
418 diag = graph_sum - graph_diagonal - graph_diagonal
420 if normed:
421 isolated_node_mask = diag == 0
422 w = np.where(isolated_node_mask, 1, np.sqrt(diag))
423 if symmetrized:
424 md = _laplace_normed_sym(graph, graph_sum, 1.0 / w)
425 else:
426 md = _laplace_normed(graph, graph_sum, 1.0 / w)
427 if form == "function":
428 return md, w.astype(dtype, copy=False)
429 elif form == "lo":
430 m = _linearoperator(md, shape=graph.shape, dtype=dtype)
431 return m, w.astype(dtype, copy=False)
432 else:
433 raise ValueError(f"Invalid form: {form!r}")
434 else:
435 if symmetrized:
436 md = _laplace_sym(graph, graph_sum)
437 else:
438 md = _laplace(graph, graph_sum)
439 if form == "function":
440 return md, diag.astype(dtype, copy=False)
441 elif form == "lo":
442 m = _linearoperator(md, shape=graph.shape, dtype=dtype)
443 return m, diag.astype(dtype, copy=False)
444 else:
445 raise ValueError(f"Invalid form: {form!r}")
448def _laplacian_sparse(graph, normed, axis, copy, form, dtype, symmetrized):
449 # The keyword argument `form` is unused and has no effect here.
450 del form
452 if dtype is None:
453 dtype = graph.dtype
455 needs_copy = False
456 if graph.format in ('lil', 'dok'):
457 m = graph.tocoo()
458 else:
459 m = graph
460 if copy:
461 needs_copy = True
463 if symmetrized:
464 m += m.T.conj()
466 w = np.asarray(m.sum(axis=axis)).ravel() - m.diagonal()
467 if normed:
468 m = m.tocoo(copy=needs_copy)
469 isolated_node_mask = (w == 0)
470 w = np.where(isolated_node_mask, 1, np.sqrt(w))
471 m.data /= w[m.row]
472 m.data /= w[m.col]
473 m.data *= -1
474 m.setdiag(1 - isolated_node_mask)
475 else:
476 if m.format == 'dia':
477 m = m.copy()
478 else:
479 m = m.tocoo(copy=needs_copy)
480 m.data *= -1
481 m.setdiag(w)
483 return m.astype(dtype, copy=False), w.astype(dtype)
486def _laplacian_dense_flo(graph, normed, axis, copy, form, dtype, symmetrized):
488 if copy:
489 m = np.array(graph)
490 else:
491 m = np.asarray(graph)
493 if dtype is None:
494 dtype = m.dtype
496 graph_sum = m.sum(axis=axis)
497 graph_diagonal = m.diagonal()
498 diag = graph_sum - graph_diagonal
499 if symmetrized:
500 graph_sum += m.sum(axis=1 - axis)
501 diag = graph_sum - graph_diagonal - graph_diagonal
503 if normed:
504 isolated_node_mask = diag == 0
505 w = np.where(isolated_node_mask, 1, np.sqrt(diag))
506 if symmetrized:
507 md = _laplace_normed_sym(m, graph_sum, 1.0 / w)
508 else:
509 md = _laplace_normed(m, graph_sum, 1.0 / w)
510 if form == "function":
511 return md, w.astype(dtype, copy=False)
512 elif form == "lo":
513 m = _linearoperator(md, shape=graph.shape, dtype=dtype)
514 return m, w.astype(dtype, copy=False)
515 else:
516 raise ValueError(f"Invalid form: {form!r}")
517 else:
518 if symmetrized:
519 md = _laplace_sym(m, graph_sum)
520 else:
521 md = _laplace(m, graph_sum)
522 if form == "function":
523 return md, diag.astype(dtype, copy=False)
524 elif form == "lo":
525 m = _linearoperator(md, shape=graph.shape, dtype=dtype)
526 return m, diag.astype(dtype, copy=False)
527 else:
528 raise ValueError(f"Invalid form: {form!r}")
531def _laplacian_dense(graph, normed, axis, copy, form, dtype, symmetrized):
533 if form != "array":
534 raise ValueError(f'{form!r} must be "array"')
536 if dtype is None:
537 dtype = graph.dtype
539 if copy:
540 m = np.array(graph)
541 else:
542 m = np.asarray(graph)
544 if dtype is None:
545 dtype = m.dtype
547 if symmetrized:
548 m += m.T.conj()
549 np.fill_diagonal(m, 0)
550 w = m.sum(axis=axis)
551 if normed:
552 isolated_node_mask = (w == 0)
553 w = np.where(isolated_node_mask, 1, np.sqrt(w))
554 m /= w
555 m /= w[:, np.newaxis]
556 m *= -1
557 _setdiag_dense(m, 1 - isolated_node_mask)
558 else:
559 m *= -1
560 _setdiag_dense(m, w)
562 return m.astype(dtype, copy=False), w.astype(dtype, copy=False)