Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/csgraph/_laplacian.py: 10%

1""" 

2Laplacian of a compressed-sparse graph 

3""" 

4 

5import numpy as np 

6from scipy.sparse import isspmatrix 

7from scipy.sparse.linalg import LinearOperator 

8 

9 

10############################################################################### 

11# Graph laplacian 

12def laplacian( 

13 csgraph, 

14 normed=False, 

15 return_diag=False, 

16 use_out_degree=False, 

17 *, 

18 copy=True, 

19 form="array", 

20 dtype=None, 

21 symmetrized=False, 

22): 

23 """ 

24 Return the Laplacian of a directed graph. 

25 

26 Parameters 

27 ---------- 

28 csgraph : array_like or sparse matrix, 2 dimensions 

29 compressed-sparse graph, with shape (N, N). 

30 normed : bool, optional 

31 If True, then compute symmetrically normalized Laplacian. 

32 Default: False. 

33 return_diag : bool, optional 

34 If True, then also return an array related to vertex degrees. 

35 Default: False. 

36 use_out_degree : bool, optional 

37 If True, then use out-degree instead of in-degree. 

38 This distinction matters only if the graph is asymmetric. 

39 Default: False. 

40 copy: bool, optional 

41 If False, then change `csgraph` in place if possible, 

42 avoiding doubling the memory use. 

43 Default: True, for backward compatibility. 

44 form: 'array', or 'function', or 'lo' 

45 Determines the format of the output Laplacian: 

46 

47 * 'array' is a numpy array; 

48 * 'function' is a pointer to evaluating the Laplacian-vector 

49 or Laplacian-matrix product; 

50 * 'lo' results in the format of the `LinearOperator`. 

51 

52 Choosing 'function' or 'lo' always avoids doubling 

53 the memory use, ignoring `copy` value. 

54 Default: 'array', for backward compatibility. 

55 dtype: None or one of numeric numpy dtypes, optional 

56 The dtype of the output. If ``dtype=None``, the dtype of the 

57 output matches the dtype of the input csgraph, except for 

58 the case ``normed=True`` and integer-like csgraph, where 

59 the output dtype is 'float' allowing accurate normalization, 

60 but dramatically increasing the memory use. 

61 Default: None, for backward compatibility. 

62 symmetrized: bool, optional 

63 If True, then the output Laplacian is symmetric/Hermitian. 

64 The symmetrization is done by ``csgraph + csgraph.T.conj`` 

65 without dividing by 2 to preserve integer dtypes if possible 

66 prior to the construction of the Laplacian. 

67 The symmetrization will increase the memory footprint of 

68 sparse matrices unless the sparsity pattern is symmetric or 

69 `form` is 'function' or 'lo'. 

70 Default: False, for backward compatibility. 

71 

72 Returns 

73 ------- 

74 lap : ndarray, or sparse matrix, or `LinearOperator` 

75 The N x N Laplacian of csgraph. It will be a NumPy array (dense) 

76 if the input was dense, or a sparse matrix otherwise, or 

77 the format of a function or `LinearOperator` if 

78 `form` equals 'function' or 'lo', respectively. 

79 diag : ndarray, optional 

80 The length-N main diagonal of the Laplacian matrix. 

81 For the normalized Laplacian, this is the array of square roots 

82 of vertex degrees or 1 if the degree is zero. 

83 

84 Notes 

85 ----- 

86 The Laplacian matrix of a graph is sometimes referred to as the 

87 "Kirchhoff matrix" or just the "Laplacian", and is useful in many 

88 parts of spectral graph theory. 

89 In particular, the eigen-decomposition of the Laplacian can give 

90 insight into many properties of the graph, e.g., 

91 is commonly used for spectral data embedding and clustering. 

92 

93 The constructed Laplacian doubles the memory use if ``copy=True`` and 

94 ``form="array"`` which is the default. 

95 Choosing ``copy=False`` has no effect unless ``form="array"`` 

96 or the matrix is sparse in the ``coo`` format, or dense array, except 

97 for the integer input with ``normed=True`` that forces the float output. 

98 

99 Sparse input is reformatted into ``coo`` if ``form="array"``, 

100 which is the default. 

101 

102 If the input adjacency matrix is not symmetic, the Laplacian is 

103 also non-symmetric unless ``symmetrized=True`` is used. 

104 

105 Diagonal entries of the input adjacency matrix are ignored and 

106 replaced with zeros for the purpose of normalization where ``normed=True``. 

107 The normalization uses the inverse square roots of row-sums of the input 

108 adjacency matrix, and thus may fail if the row-sums contain 

109 negative or complex with a non-zero imaginary part values. 

110 

111 The normalization is symmetric, making the normalized Laplacian also 

112 symmetric if the input csgraph was symmetric. 

113 

114 References 

115 ---------- 

116 .. [1] Laplacian matrix. https://en.wikipedia.org/wiki/Laplacian_matrix 

117 

118 Examples 

119 -------- 

120 >>> import numpy as np 

121 >>> from scipy.sparse import csgraph 

122 

123 Our first illustration is the symmetric graph 

124 

125 >>> G = np.arange(4) * np.arange(4)[:, np.newaxis] 

126 >>> G 

127 array([[0, 0, 0, 0], 

128 [0, 1, 2, 3], 

129 [0, 2, 4, 6], 

130 [0, 3, 6, 9]]) 

131 

132 and its symmetric Laplacian matrix 

133 

134 >>> csgraph.laplacian(G) 

135 array([[ 0, 0, 0, 0], 

136 [ 0, 5, -2, -3], 

137 [ 0, -2, 8, -6], 

138 [ 0, -3, -6, 9]]) 

139 

140 The non-symmetric graph 

141 

142 >>> G = np.arange(9).reshape(3, 3) 

143 >>> G 

144 array([[0, 1, 2], 

145 [3, 4, 5], 

146 [6, 7, 8]]) 

147 

148 has different row- and column sums, resulting in two varieties 

149 of the Laplacian matrix, using an in-degree, which is the default 

150 

151 >>> L_in_degree = csgraph.laplacian(G) 

152 >>> L_in_degree 

153 array([[ 9, -1, -2], 

154 [-3, 8, -5], 

155 [-6, -7, 7]]) 

156 

157 or alternatively an out-degree 

158 

159 >>> L_out_degree = csgraph.laplacian(G, use_out_degree=True) 

160 >>> L_out_degree 

161 array([[ 3, -1, -2], 

162 [-3, 8, -5], 

163 [-6, -7, 13]]) 

164 

165 Constructing a symmetric Laplacian matrix, one can add the two as 

166 

167 >>> L_in_degree + L_out_degree.T 

168 array([[ 12, -4, -8], 

169 [ -4, 16, -12], 

170 [ -8, -12, 20]]) 

171 

172 or use the ``symmetrized=True`` option 

173 

174 >>> csgraph.laplacian(G, symmetrized=True) 

175 array([[ 12, -4, -8], 

176 [ -4, 16, -12], 

177 [ -8, -12, 20]]) 

178 

179 that is equivalent to symmetrizing the original graph 

180 

181 >>> csgraph.laplacian(G + G.T) 

182 array([[ 12, -4, -8], 

183 [ -4, 16, -12], 

184 [ -8, -12, 20]]) 

185 

186 The goal of normalization is to make the non-zero diagonal entries 

187 of the Laplacian matrix to be all unit, also scaling off-diagonal 

188 entries correspondingly. The normalization can be done manually, e.g., 

189 

190 >>> G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]) 

191 >>> L, d = csgraph.laplacian(G, return_diag=True) 

192 >>> L 

193 array([[ 2, -1, -1], 

194 [-1, 2, -1], 

195 [-1, -1, 2]]) 

196 >>> d 

197 array([2, 2, 2]) 

198 >>> scaling = np.sqrt(d) 

199 >>> scaling 

200 array([1.41421356, 1.41421356, 1.41421356]) 

201 >>> (1/scaling)*L*(1/scaling) 

202 array([[ 1. , -0.5, -0.5], 

203 [-0.5, 1. , -0.5], 

204 [-0.5, -0.5, 1. ]]) 

205 

206 Or using ``normed=True`` option 

207 

208 >>> L, d = csgraph.laplacian(G, return_diag=True, normed=True) 

209 >>> L 

210 array([[ 1. , -0.5, -0.5], 

211 [-0.5, 1. , -0.5], 

212 [-0.5, -0.5, 1. ]]) 

213 

214 which now instead of the diagonal returns the scaling coefficients 

215 

216 >>> d 

217 array([1.41421356, 1.41421356, 1.41421356]) 

218 

219 Zero scaling coefficients are substituted with 1s, where scaling 

220 has thus no effect, e.g., 

221 

222 >>> G = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]) 

223 >>> G 

224 array([[0, 0, 0], 

225 [0, 0, 1], 

226 [0, 1, 0]]) 

227 >>> L, d = csgraph.laplacian(G, return_diag=True, normed=True) 

228 >>> L 

229 array([[ 0., -0., -0.], 

230 [-0., 1., -1.], 

231 [-0., -1., 1.]]) 

232 >>> d 

233 array([1., 1., 1.]) 

234 

235 Only the symmetric normalization is implemented, resulting 

236 in a symmetric Laplacian matrix if and only if its graph is symmetric 

237 and has all non-negative degrees, like in the examples above. 

238 

239 The output Laplacian matrix is by default a dense array or a sparse matrix 

240 inferring its shape, format, and dtype from the input graph matrix: 

241 

242 >>> G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]).astype(np.float32) 

243 >>> G 

244 array([[0., 1., 1.], 

245 [1., 0., 1.], 

246 [1., 1., 0.]], dtype=float32) 

247 >>> csgraph.laplacian(G) 

248 array([[ 2., -1., -1.], 

249 [-1., 2., -1.], 

250 [-1., -1., 2.]], dtype=float32) 

251 

252 but can alternatively be generated matrix-free as a LinearOperator: 

253 

254 >>> L = csgraph.laplacian(G, form="lo") 

255 >>> L 

256 <3x3 _CustomLinearOperator with dtype=float32> 

257 >>> L(np.eye(3)) 

258 array([[ 2., -1., -1.], 

259 [-1., 2., -1.], 

260 [-1., -1., 2.]]) 

261 

262 or as a lambda-function: 

263 

264 >>> L = csgraph.laplacian(G, form="function") 

265 >>> L 

266 <function _laplace.<locals>.<lambda> at 0x0000012AE6F5A598> 

267 >>> L(np.eye(3)) 

268 array([[ 2., -1., -1.], 

269 [-1., 2., -1.], 

270 [-1., -1., 2.]]) 

271 

272 The Laplacian matrix is used for 

273 spectral data clustering and embedding 

274 as well as for spectral graph partitioning. 

275 Our final example illustrates the latter 

276 for a noisy directed linear graph. 

277 

278 >>> from scipy.sparse import diags, random 

279 >>> from scipy.sparse.linalg import lobpcg 

280 

281 Create a directed linear graph with ``N=35`` vertices 

282 using a sparse adjacency matrix ``G``: 

283 

284 >>> N = 35 

285 >>> G = diags(np.ones(N-1), 1, format="csr") 

286 

287 Fix a random seed ``rng`` and add a random sparse noise to the graph ``G``: 

288 

289 >>> rng = np.random.default_rng() 

290 >>> G += 1e-2 * random(N, N, density=0.1, random_state=rng) 

291 

292 Set initial approximations for eigenvectors: 

293 

294 >>> X = rng.random((N, 2)) 

295 

296 The constant vector of ones is always a trivial eigenvector 

297 of the non-normalized Laplacian to be filtered out: 

298 

299 >>> Y = np.ones((N, 1)) 

300 

301 Alternating (1) the sign of the graph weights allows determining 

302 labels for spectral max- and min- cuts in a single loop. 

303 Since the graph is undirected, the option ``symmetrized=True`` 

304 must be used in the construction of the Laplacian. 

305 The option ``normed=True`` cannot be used in (2) for the negative weights 

306 here as the symmetric normalization evaluates square roots. 

307 The option ``form="lo"`` in (2) is matrix-free, i.e., guarantees 

308 a fixed memory footprint and read-only access to the graph. 

309 Calling the eigenvalue solver ``lobpcg`` (3) computes the Fiedler vector 

310 that determines the labels as the signs of its components in (5). 

311 Since the sign in an eigenvector is not deterministic and can flip, 

312 we fix the sign of the first component to be always +1 in (4). 

313 

314 >>> for cut in ["max", "min"]: 

315 ... G = -G # 1. 

316 ... L = csgraph.laplacian(G, symmetrized=True, form="lo") # 2. 

317 ... _, eves = lobpcg(L, X, Y=Y, largest=False, tol=1e-3) # 3. 

318 ... eves *= np.sign(eves[0, 0]) # 4. 

319 ... print(cut + "-cut labels:\\n", 1 * (eves[:, 0]>0)) # 5. 

320 max-cut labels: 

321 [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] 

322 min-cut labels: 

323 [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] 

324 

325 As anticipated for a (slightly noisy) linear graph, 

326 the max-cut strips all the edges of the graph coloring all 

327 odd vertices into one color and all even vertices into another one, 

328 while the balanced min-cut partitions the graph 

329 in the middle by deleting a single edge. 

330 Both determined partitions are optimal. 

331 """ 

332 if csgraph.ndim != 2 or csgraph.shape[0] != csgraph.shape[1]: 

333 raise ValueError('csgraph must be a square matrix or array') 

334 

335 if normed and ( 

336 np.issubdtype(csgraph.dtype, np.signedinteger) 

337 or np.issubdtype(csgraph.dtype, np.uint) 

338 ): 

339 csgraph = csgraph.astype(np.float64) 

340 

341 if form == "array": 

342 create_lap = ( 

343 _laplacian_sparse if isspmatrix(csgraph) else _laplacian_dense 

344 ) 

345 else: 

346 create_lap = ( 

347 _laplacian_sparse_flo 

348 if isspmatrix(csgraph) 

349 else _laplacian_dense_flo 

350 ) 

351 

352 degree_axis = 1 if use_out_degree else 0 

353 

354 lap, d = create_lap( 

355 csgraph, 

356 normed=normed, 

357 axis=degree_axis, 

358 copy=copy, 

359 form=form, 

360 dtype=dtype, 

361 symmetrized=symmetrized, 

362 ) 

363 if return_diag: 

364 return lap, d 

365 return lap 

366 

367 

368def _setdiag_dense(m, d): 

369 step = len(d) + 1 

370 m.flat[::step] = d 

371 

372 

373def _laplace(m, d): 

374 return lambda v: v * d[:, np.newaxis] - m @ v 

375 

376 

377def _laplace_normed(m, d, nd): 

378 laplace = _laplace(m, d) 

379 return lambda v: nd[:, np.newaxis] * laplace(v * nd[:, np.newaxis]) 

380 

381 

382def _laplace_sym(m, d): 

383 return ( 

384 lambda v: v * d[:, np.newaxis] 

385 - m @ v 

386 - np.transpose(np.conjugate(np.transpose(np.conjugate(v)) @ m)) 

387 ) 

388 

389 

390def _laplace_normed_sym(m, d, nd): 

391 laplace_sym = _laplace_sym(m, d) 

392 return lambda v: nd[:, np.newaxis] * laplace_sym(v * nd[:, np.newaxis]) 

393 

394 

395def _linearoperator(mv, shape, dtype): 

396 return LinearOperator(matvec=mv, matmat=mv, shape=shape, dtype=dtype) 

397 

398 

399def _laplacian_sparse_flo(graph, normed, axis, copy, form, dtype, symmetrized): 

400 # The keyword argument `copy` is unused and has no effect here. 

401 del copy 

402 

403 if dtype is None: 

404 dtype = graph.dtype 

405 

406 graph_sum = graph.sum(axis=axis).getA1() 

407 graph_diagonal = graph.diagonal() 

408 diag = graph_sum - graph_diagonal 

409 if symmetrized: 

410 graph_sum += graph.sum(axis=1 - axis).getA1() 

411 diag = graph_sum - graph_diagonal - graph_diagonal 

412 

413 if normed: 

414 isolated_node_mask = diag == 0 

415 w = np.where(isolated_node_mask, 1, np.sqrt(diag)) 

416 if symmetrized: 

417 md = _laplace_normed_sym(graph, graph_sum, 1.0 / w) 

418 else: 

419 md = _laplace_normed(graph, graph_sum, 1.0 / w) 

420 if form == "function": 

421 return md, w.astype(dtype, copy=False) 

422 elif form == "lo": 

423 m = _linearoperator(md, shape=graph.shape, dtype=dtype) 

424 return m, w.astype(dtype, copy=False) 

425 else: 

426 raise ValueError(f"Invalid form: {form!r}") 

427 else: 

428 if symmetrized: 

429 md = _laplace_sym(graph, graph_sum) 

430 else: 

431 md = _laplace(graph, graph_sum) 

432 if form == "function": 

433 return md, diag.astype(dtype, copy=False) 

434 elif form == "lo": 

435 m = _linearoperator(md, shape=graph.shape, dtype=dtype) 

436 return m, diag.astype(dtype, copy=False) 

437 else: 

438 raise ValueError(f"Invalid form: {form!r}") 

439 

440 

441def _laplacian_sparse(graph, normed, axis, copy, form, dtype, symmetrized): 

442 # The keyword argument `form` is unused and has no effect here. 

443 del form 

444 

445 if dtype is None: 

446 dtype = graph.dtype 

447 

448 needs_copy = False 

449 if graph.format in ('lil', 'dok'): 

450 m = graph.tocoo() 

451 else: 

452 m = graph 

453 if copy: 

454 needs_copy = True 

455 

456 if symmetrized: 

457 m += m.T.conj() 

458 

459 w = m.sum(axis=axis).getA1() - m.diagonal() 

460 if normed: 

461 m = m.tocoo(copy=needs_copy) 

462 isolated_node_mask = (w == 0) 

463 w = np.where(isolated_node_mask, 1, np.sqrt(w)) 

464 m.data /= w[m.row] 

465 m.data /= w[m.col] 

466 m.data *= -1 

467 m.setdiag(1 - isolated_node_mask) 

468 else: 

469 if m.format == 'dia': 

470 m = m.copy() 

471 else: 

472 m = m.tocoo(copy=needs_copy) 

473 m.data *= -1 

474 m.setdiag(w) 

475 

476 return m.astype(dtype, copy=False), w.astype(dtype) 

477 

478 

479def _laplacian_dense_flo(graph, normed, axis, copy, form, dtype, symmetrized): 

480 

481 if copy: 

482 m = np.array(graph) 

483 else: 

484 m = np.asarray(graph) 

485 

486 if dtype is None: 

487 dtype = m.dtype 

488 

489 graph_sum = m.sum(axis=axis) 

490 graph_diagonal = m.diagonal() 

491 diag = graph_sum - graph_diagonal 

492 if symmetrized: 

493 graph_sum += m.sum(axis=1 - axis) 

494 diag = graph_sum - graph_diagonal - graph_diagonal 

495 

496 if normed: 

497 isolated_node_mask = diag == 0 

498 w = np.where(isolated_node_mask, 1, np.sqrt(diag)) 

499 if symmetrized: 

500 md = _laplace_normed_sym(m, graph_sum, 1.0 / w) 

501 else: 

502 md = _laplace_normed(m, graph_sum, 1.0 / w) 

503 if form == "function": 

504 return md, w.astype(dtype, copy=False) 

505 elif form == "lo": 

506 m = _linearoperator(md, shape=graph.shape, dtype=dtype) 

507 return m, w.astype(dtype, copy=False) 

508 else: 

509 raise ValueError(f"Invalid form: {form!r}") 

510 else: 

511 if symmetrized: 

512 md = _laplace_sym(m, graph_sum) 

513 else: 

514 md = _laplace(m, graph_sum) 

515 if form == "function": 

516 return md, diag.astype(dtype, copy=False) 

517 elif form == "lo": 

518 m = _linearoperator(md, shape=graph.shape, dtype=dtype) 

519 return m, diag.astype(dtype, copy=False) 

520 else: 

521 raise ValueError(f"Invalid form: {form!r}") 

522 

523 

524def _laplacian_dense(graph, normed, axis, copy, form, dtype, symmetrized): 

525 

526 if form != "array": 

527 raise ValueError(f'{form!r} must be "array"') 

528 

529 if dtype is None: 

530 dtype = graph.dtype 

531 

532 if copy: 

533 m = np.array(graph) 

534 else: 

535 m = np.asarray(graph) 

536 

537 if dtype is None: 

538 dtype = m.dtype 

539 

540 if symmetrized: 

541 m += m.T.conj() 

542 np.fill_diagonal(m, 0) 

543 w = m.sum(axis=axis) 

544 if normed: 

545 isolated_node_mask = (w == 0) 

546 w = np.where(isolated_node_mask, 1, np.sqrt(w)) 

547 m /= w 

548 m /= w[:, np.newaxis] 

549 m *= -1 

550 _setdiag_dense(m, 1 - isolated_node_mask) 

551 else: 

552 m *= -1 

553 _setdiag_dense(m, w) 

554 

555 return m.astype(dtype, copy=False), w.astype(dtype, copy=False)