Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_special_sparse_arrays.py: 12%
142 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +0000
1import numpy as np
2from scipy.sparse.linalg import LinearOperator
3from scipy.sparse import kron, eye, dia_array
5__all__ = ['LaplacianNd']
8class LaplacianNd(LinearOperator):
9 """
10 The grid Laplacian in ``N`` dimensions and its eigenvalues/eigenvectors.
12 Construct Laplacian on a uniform rectangular grid in `N` dimensions
13 and output its eigenvalues and eigenvectors.
14 The Laplacian ``L`` is square, negative definite, real symmetric array
15 with signed integer entries and zeros otherwise.
17 Parameters
18 ----------
19 grid_shape : tuple
20 A tuple of integers of length ``N`` (corresponding to the dimension of
21 the Lapacian), where each entry gives the size of that dimension. The
22 Laplacian matrix is square of the size ``np.prod(grid_shape)``.
23 boundary_conditions : {'neumann', 'dirichlet', 'periodic'}, optional
24 The type of the boundary conditions on the boundaries of the grid.
25 Valid values are ``'dirichlet'`` or ``'neumann'``(default) or
26 ``'periodic'``.
27 dtype : dtype
28 Numerical type of the array. Default is ``np.int8``.
30 Methods
31 -------
32 toarray()
33 Construct a dense array from Laplacian data
34 tosparse()
35 Construct a sparse array from Laplacian data
36 eigenvalues(m=None)
37 Construct a 1D array of `m` largest (smallest in absolute value)
38 eigenvalues of the Laplacian matrix in ascending order.
39 eigenvectors(m=None):
40 Construct the array with columns made of `m` eigenvectors (``float``)
41 of the ``Nd`` Laplacian corresponding to the `m` ordered eigenvalues.
43 .. versionadded:: 1.12.0
45 Notes
46 -----
47 Compared to the MATLAB/Octave implementation [1] of 1-, 2-, and 3-D
48 Laplacian, this code allows the arbitrary N-D case and the matrix-free
49 callable option, but is currently limited to pure Dirichlet, Neumann or
50 Periodic boundary conditions only.
52 The Laplacian matrix of a graph (`scipy.sparse.csgraph.laplacian`) of a
53 rectangular grid corresponds to the negative Laplacian with the Neumann
54 conditions, i.e., ``boundary_conditions = 'neumann'``.
56 All eigenvalues and eigenvectors of the discrete Laplacian operator for
57 an ``N``-dimensional regular grid of shape `grid_shape` with the grid
58 step size ``h=1`` are analytically known [2].
60 References
61 ----------
62 .. [1] https://github.com/lobpcg/blopex/blob/master/blopex_\
63tools/matlab/laplacian/laplacian.m
64 .. [2] "Eigenvalues and eigenvectors of the second derivative", Wikipedia
65 https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_\
66of_the_second_derivative
68 Examples
69 --------
70 >>> import numpy as np
71 >>> from scipy.sparse.linalg import LaplacianNd
72 >>> from scipy.sparse import diags, csgraph
73 >>> from scipy.linalg import eigvalsh
75 The one-dimensional Laplacian demonstrated below for pure Neumann boundary
76 conditions on a regular grid with ``n=6`` grid points is exactly the
77 negative graph Laplacian for the undirected linear graph with ``n``
78 vertices using the sparse adjacency matrix ``G`` represented by the
79 famous tri-diagonal matrix:
81 >>> n = 6
82 >>> G = diags(np.ones(n - 1), 1, format='csr')
83 >>> Lf = csgraph.laplacian(G, symmetrized=True, form='function')
84 >>> grid_shape = (n, )
85 >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
86 >>> np.array_equal(lap.matmat(np.eye(n)), -Lf(np.eye(n)))
87 True
89 Since all matrix entries of the Laplacian are integers, ``'int8'`` is
90 the default dtype for storing matrix representations.
92 >>> lap.tosparse()
93 <6x6 sparse matrix of type '<class 'numpy.int8'>'
94 with 16 stored elements (3 diagonals) in DIAgonal format>
95 >>> lap.toarray()
96 array([[-1, 1, 0, 0, 0, 0],
97 [ 1, -2, 1, 0, 0, 0],
98 [ 0, 1, -2, 1, 0, 0],
99 [ 0, 0, 1, -2, 1, 0],
100 [ 0, 0, 0, 1, -2, 1],
101 [ 0, 0, 0, 0, 1, -1]], dtype=int8)
102 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
103 True
104 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
105 True
107 Any number of extreme eigenvalues and/or eigenvectors can be computed.
109 >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
110 >>> lap.eigenvalues()
111 array([-4., -3., -3., -1., -1., 0.])
112 >>> lap.eigenvalues()[-2:]
113 array([-1., 0.])
114 >>> lap.eigenvalues(2)
115 array([-1., 0.])
116 >>> lap.eigenvectors(1)
117 array([[0.40824829],
118 [0.40824829],
119 [0.40824829],
120 [0.40824829],
121 [0.40824829],
122 [0.40824829]])
123 >>> lap.eigenvectors(2)
124 array([[ 0.5 , 0.40824829],
125 [ 0. , 0.40824829],
126 [-0.5 , 0.40824829],
127 [-0.5 , 0.40824829],
128 [ 0. , 0.40824829],
129 [ 0.5 , 0.40824829]])
130 >>> lap.eigenvectors()
131 array([[ 0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 ,
132 0.40824829],
133 [-0.40824829, -0.57735027, -0.57735027, 0. , 0. ,
134 0.40824829],
135 [ 0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 ,
136 0.40824829],
137 [-0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 ,
138 0.40824829],
139 [ 0.40824829, -0.57735027, -0.57735027, 0. , 0. ,
140 0.40824829],
141 [-0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 ,
142 0.40824829]])
144 The two-dimensional Laplacian is illustrated on a regular grid with
145 ``grid_shape = (2, 3)`` points in each dimension.
147 >>> grid_shape = (2, 3)
148 >>> n = np.prod(grid_shape)
150 Numeration of grid points is as follows:
152 >>> np.arange(n).reshape(grid_shape + (-1,))
153 array([[[0],
154 [1],
155 [2]],
156 <BLANKLINE>
157 [[3],
158 [4],
159 [5]]])
161 Each of the boundary conditions ``'dirichlet'``, ``'periodic'``, and
162 ``'neumann'`` is illustrated separately; with ``'dirichlet'``
164 >>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet')
165 >>> lap.tosparse()
166 <6x6 sparse array of type '<class 'numpy.int8'>'
167 with 20 stored elements in Compressed Sparse Row format>
168 >>> lap.toarray()
169 array([[-4, 1, 0, 1, 0, 0],
170 [ 1, -4, 1, 0, 1, 0],
171 [ 0, 1, -4, 0, 0, 1],
172 [ 1, 0, 0, -4, 1, 0],
173 [ 0, 1, 0, 1, -4, 1],
174 [ 0, 0, 1, 0, 1, -4]], dtype=int8)
175 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
176 True
177 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
178 True
179 >>> lap.eigenvalues()
180 array([-6.41421356, -5. , -4.41421356, -3.58578644, -3. ,
181 -1.58578644])
182 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
183 >>> np.allclose(lap.eigenvalues(), eigvals)
184 True
185 >>> np.allclose(lap.toarray() @ lap.eigenvectors(),
186 ... lap.eigenvectors() @ np.diag(lap.eigenvalues()))
187 True
189 with ``'periodic'``
191 >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic')
192 >>> lap.tosparse()
193 <6x6 sparse array of type '<class 'numpy.int8'>'
194 with 24 stored elements in Compressed Sparse Row format>
195 >>> lap.toarray()
196 array([[-4, 1, 1, 2, 0, 0],
197 [ 1, -4, 1, 0, 2, 0],
198 [ 1, 1, -4, 0, 0, 2],
199 [ 2, 0, 0, -4, 1, 1],
200 [ 0, 2, 0, 1, -4, 1],
201 [ 0, 0, 2, 1, 1, -4]], dtype=int8)
202 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
203 True
204 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
205 True
206 >>> lap.eigenvalues()
207 array([-7., -7., -4., -3., -3., 0.])
208 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
209 >>> np.allclose(lap.eigenvalues(), eigvals)
210 True
211 >>> np.allclose(lap.toarray() @ lap.eigenvectors(),
212 ... lap.eigenvectors() @ np.diag(lap.eigenvalues()))
213 True
215 and with ``'neumann'``
217 >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann')
218 >>> lap.tosparse()
219 <6x6 sparse array of type '<class 'numpy.int8'>'
220 with 20 stored elements in Compressed Sparse Row format>
221 >>> lap.toarray()
222 array([[-2, 1, 0, 1, 0, 0],
223 [ 1, -3, 1, 0, 1, 0],
224 [ 0, 1, -2, 0, 0, 1],
225 [ 1, 0, 0, -2, 1, 0],
226 [ 0, 1, 0, 1, -3, 1],
227 [ 0, 0, 1, 0, 1, -2]])
228 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray())
229 True
230 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray())
231 True
232 >>> lap.eigenvalues()
233 array([-5., -3., -3., -2., -1., 0.])
234 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64))
235 >>> np.allclose(lap.eigenvalues(), eigvals)
236 True
237 >>> np.allclose(lap.toarray() @ lap.eigenvectors(),
238 ... lap.eigenvectors() @ np.diag(lap.eigenvalues()))
239 True
241 """
243 def __init__(self, grid_shape, *,
244 boundary_conditions='neumann',
245 dtype=np.int8):
247 if boundary_conditions not in ('dirichlet', 'neumann', 'periodic'):
248 raise ValueError(
249 f"Unknown value {boundary_conditions!r} is given for "
250 "'boundary_conditions' parameter. The valid options are "
251 "'dirichlet', 'periodic', and 'neumann' (default)."
252 )
254 self.grid_shape = grid_shape
255 self.boundary_conditions = boundary_conditions
256 # LaplacianNd folds all dimensions in `grid_shape` into a single one
257 N = np.prod(grid_shape)
258 super().__init__(dtype=dtype, shape=(N, N))
260 def _eigenvalue_ordering(self, m):
261 """Compute `m` largest eigenvalues in each of the ``N`` directions,
262 i.e., up to ``m * N`` total, order them and return `m` largest.
263 """
264 grid_shape = self.grid_shape
265 if m is None:
266 indices = np.indices(grid_shape)
267 Leig = np.zeros(grid_shape)
268 else:
269 grid_shape_min = min(grid_shape,
270 tuple(np.ones_like(grid_shape) * m))
271 indices = np.indices(grid_shape_min)
272 Leig = np.zeros(grid_shape_min)
274 for j, n in zip(indices, grid_shape):
275 if self.boundary_conditions == 'dirichlet':
276 Leig += -4 * np.sin(np.pi * (j + 1) / (2 * (n + 1))) ** 2
277 elif self.boundary_conditions == 'neumann':
278 Leig += -4 * np.sin(np.pi * j / (2 * n)) ** 2
279 else: # boundary_conditions == 'periodic'
280 Leig += -4 * np.sin(np.pi * np.floor((j + 1) / 2) / n) ** 2
282 Leig_ravel = Leig.ravel()
283 ind = np.argsort(Leig_ravel)
284 eigenvalues = Leig_ravel[ind]
285 if m is not None:
286 eigenvalues = eigenvalues[-m:]
287 ind = ind[-m:]
289 return eigenvalues, ind
291 def eigenvalues(self, m=None):
292 """Return the requested number of eigenvalues.
294 Parameters
295 ----------
296 m : int, optional
297 The positive number of eigenvalues to return. If not provided,
298 then all eigenvalues will be returned.
300 Returns
301 -------
302 eigenvalues : float array
303 The requested `m` or all eigenvalues, in ascending order.
304 """
305 eigenvalues, _ = self._eigenvalue_ordering(m)
306 return eigenvalues
308 def _ev1d(self, j, n):
309 """Return 1 eigenvector in 1d with index `j`
310 and number of grid points `n` where ``j < n``.
311 """
312 if self.boundary_conditions == 'dirichlet':
313 i = np.pi * (np.arange(n) + 1) / (n + 1)
314 ev = np.sqrt(2. / (n + 1.)) * np.sin(i * (j + 1))
315 elif self.boundary_conditions == 'neumann':
316 i = np.pi * (np.arange(n) + 0.5) / n
317 ev = np.sqrt((1. if j == 0 else 2.) / n) * np.cos(i * j)
318 else: # boundary_conditions == 'periodic'
319 if j == 0:
320 ev = np.sqrt(1. / n) * np.ones(n)
321 elif j + 1 == n and n % 2 == 0:
322 ev = np.sqrt(1. / n) * np.tile([1, -1], n//2)
323 else:
324 i = 2. * np.pi * (np.arange(n) + 0.5) / n
325 ev = np.sqrt(2. / n) * np.cos(i * np.floor((j + 1) / 2))
326 # make small values exact zeros correcting round-off errors
327 # due to symmetry of eigenvectors the exact 0. is correct
328 ev[np.abs(ev) < np.finfo(np.float64).eps] = 0.
329 return ev
331 def _one_eve(self, k):
332 """Return 1 eigenvector in Nd with multi-index `j`
333 as a tensor product of the corresponding 1d eigenvectors.
334 """
335 phi = [self._ev1d(j, n) for j, n in zip(k, self.grid_shape)]
336 result = phi[0]
337 for phi in phi[1:]:
338 result = np.tensordot(result, phi, axes=0)
339 return np.asarray(result).ravel()
341 def eigenvectors(self, m=None):
342 """Return the requested number of eigenvectors for ordered eigenvalues.
344 Parameters
345 ----------
346 m : int, optional
347 The positive number of eigenvectors to return. If not provided,
348 then all eigenvectors will be returned.
350 Returns
351 -------
352 eigenvectors : float array
353 An array with columns made of the requested `m` or all eigenvectors.
354 The columns are ordered according to the `m` ordered eigenvalues.
355 """
356 _, ind = self._eigenvalue_ordering(m)
357 if m is None:
358 grid_shape_min = self.grid_shape
359 else:
360 grid_shape_min = min(self.grid_shape,
361 tuple(np.ones_like(self.grid_shape) * m))
363 N_indices = np.unravel_index(ind, grid_shape_min)
364 N_indices = [tuple(x) for x in zip(*N_indices)]
365 eigenvectors_list = [self._one_eve(k) for k in N_indices]
366 return np.column_stack(eigenvectors_list)
368 def toarray(self):
369 """
370 Converts the Laplacian data to a dense array.
372 Returns
373 -------
374 L : ndarray
375 The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``.
377 """
378 grid_shape = self.grid_shape
379 n = np.prod(grid_shape)
380 L = np.zeros([n, n], dtype=np.int8)
381 # Scratch arrays
382 L_i = np.empty_like(L)
383 Ltemp = np.empty_like(L)
385 for ind, dim in enumerate(grid_shape):
386 # Start zeroing out L_i
387 L_i[:] = 0
388 # Allocate the top left corner with the kernel of L_i
389 # Einsum returns writable view of arrays
390 np.einsum("ii->i", L_i[:dim, :dim])[:] = -2
391 np.einsum("ii->i", L_i[: dim - 1, 1:dim])[:] = 1
392 np.einsum("ii->i", L_i[1:dim, : dim - 1])[:] = 1
394 if self.boundary_conditions == 'neumann':
395 L_i[0, 0] = -1
396 L_i[dim - 1, dim - 1] = -1
397 elif self.boundary_conditions == 'periodic':
398 if dim > 1:
399 L_i[0, dim - 1] += 1
400 L_i[dim - 1, 0] += 1
401 else:
402 L_i[0, 0] += 1
404 # kron is too slow for large matrices hence the next two tricks
405 # 1- kron(eye, mat) is block_diag(mat, mat, ...)
406 # 2- kron(mat, eye) can be performed by 4d stride trick
408 # 1-
409 new_dim = dim
410 # for block_diag we tile the top left portion on the diagonal
411 if ind > 0:
412 tiles = np.prod(grid_shape[:ind])
413 for j in range(1, tiles):
414 L_i[j*dim:(j+1)*dim, j*dim:(j+1)*dim] = L_i[:dim, :dim]
415 new_dim += dim
416 # 2-
417 # we need the keep L_i, but reset the array
418 Ltemp[:new_dim, :new_dim] = L_i[:new_dim, :new_dim]
419 tiles = int(np.prod(grid_shape[ind+1:]))
420 # Zero out the top left, the rest is already 0
421 L_i[:new_dim, :new_dim] = 0
422 idx = [x for x in range(tiles)]
423 L_i.reshape(
424 (new_dim, tiles,
425 new_dim, tiles)
426 )[:, idx, :, idx] = Ltemp[:new_dim, :new_dim]
428 L += L_i
430 return L.astype(self.dtype)
432 def tosparse(self):
433 """
434 Constructs a sparse array from the Laplacian data. The returned sparse
435 array format is dependent on the selected boundary conditions.
437 Returns
438 -------
439 L : scipy.sparse.sparray
440 The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``.
442 """
443 N = len(self.grid_shape)
444 p = np.prod(self.grid_shape)
445 L = dia_array((p, p), dtype=np.int8)
447 for i in range(N):
448 dim = self.grid_shape[i]
449 data = np.ones([3, dim], dtype=np.int8)
450 data[1, :] *= -2
452 if self.boundary_conditions == 'neumann':
453 data[1, 0] = -1
454 data[1, -1] = -1
456 L_i = dia_array((data, [-1, 0, 1]), shape=(dim, dim),
457 dtype=np.int8
458 )
460 if self.boundary_conditions == 'periodic':
461 t = dia_array((dim, dim), dtype=np.int8)
462 t.setdiag([1], k=-dim+1)
463 t.setdiag([1], k=dim-1)
464 L_i += t
466 for j in range(i):
467 L_i = kron(eye(self.grid_shape[j], dtype=np.int8), L_i)
468 for j in range(i + 1, N):
469 L_i = kron(L_i, eye(self.grid_shape[j], dtype=np.int8))
470 L += L_i
471 return L.astype(self.dtype)
473 def _matvec(self, x):
474 grid_shape = self.grid_shape
475 N = len(grid_shape)
476 X = x.reshape(grid_shape + (-1,))
477 Y = -2 * N * X
478 for i in range(N):
479 Y += np.roll(X, 1, axis=i)
480 Y += np.roll(X, -1, axis=i)
481 if self.boundary_conditions in ('neumann', 'dirichlet'):
482 Y[(slice(None),)*i + (0,) + (slice(None),)*(N-i-1)
483 ] -= np.roll(X, 1, axis=i)[
484 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1)
485 ]
486 Y[
487 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1)
488 ] -= np.roll(X, -1, axis=i)[
489 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1)
490 ]
492 if self.boundary_conditions == 'neumann':
493 Y[
494 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1)
495 ] += np.roll(X, 0, axis=i)[
496 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1)
497 ]
498 Y[
499 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1)
500 ] += np.roll(X, 0, axis=i)[
501 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1)
502 ]
504 return Y.reshape(-1, X.shape[-1])
506 def _matmat(self, x):
507 return self._matvec(x)
509 def _adjoint(self):
510 return self
512 def _transpose(self):
513 return self