Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/sparse/linalg/_special_sparse_arrays.py: 20%

1import numpy as np 

2from scipy.sparse.linalg import LinearOperator 

3from scipy.sparse import kron, eye, dia_array 

4 

5__all__ = ['LaplacianNd'] 

6# Sakurai and Mikota classes are intended for tests and benchmarks 

7# and explicitly not included in the public API of this module. 

8 

9 

10class LaplacianNd(LinearOperator): 

11 """ 

12 The grid Laplacian in ``N`` dimensions and its eigenvalues/eigenvectors. 

13 

14 Construct Laplacian on a uniform rectangular grid in `N` dimensions 

15 and output its eigenvalues and eigenvectors. 

16 The Laplacian ``L`` is square, negative definite, real symmetric array 

17 with signed integer entries and zeros otherwise. 

18 

19 Parameters 

20 ---------- 

21 grid_shape : tuple 

22 A tuple of integers of length ``N`` (corresponding to the dimension of 

23 the Lapacian), where each entry gives the size of that dimension. The 

24 Laplacian matrix is square of the size ``np.prod(grid_shape)``. 

25 boundary_conditions : {'neumann', 'dirichlet', 'periodic'}, optional 

26 The type of the boundary conditions on the boundaries of the grid. 

27 Valid values are ``'dirichlet'`` or ``'neumann'``(default) or 

28 ``'periodic'``. 

29 dtype : dtype 

30 Numerical type of the array. Default is ``np.int8``. 

31 

32 Methods 

33 ------- 

34 toarray() 

35 Construct a dense array from Laplacian data 

36 tosparse() 

37 Construct a sparse array from Laplacian data 

38 eigenvalues(m=None) 

39 Construct a 1D array of `m` largest (smallest in absolute value) 

40 eigenvalues of the Laplacian matrix in ascending order. 

41 eigenvectors(m=None): 

42 Construct the array with columns made of `m` eigenvectors (``float``) 

43 of the ``Nd`` Laplacian corresponding to the `m` ordered eigenvalues. 

44 

45 .. versionadded:: 1.12.0 

46 

47 Notes 

48 ----- 

49 Compared to the MATLAB/Octave implementation [1] of 1-, 2-, and 3-D 

50 Laplacian, this code allows the arbitrary N-D case and the matrix-free 

51 callable option, but is currently limited to pure Dirichlet, Neumann or 

52 Periodic boundary conditions only. 

53 

54 The Laplacian matrix of a graph (`scipy.sparse.csgraph.laplacian`) of a 

55 rectangular grid corresponds to the negative Laplacian with the Neumann 

56 conditions, i.e., ``boundary_conditions = 'neumann'``. 

57 

58 All eigenvalues and eigenvectors of the discrete Laplacian operator for 

59 an ``N``-dimensional regular grid of shape `grid_shape` with the grid 

60 step size ``h=1`` are analytically known [2]. 

61 

62 References 

63 ---------- 

64 .. [1] https://github.com/lobpcg/blopex/blob/master/blopex_\ 

65tools/matlab/laplacian/laplacian.m 

66 .. [2] "Eigenvalues and eigenvectors of the second derivative", Wikipedia 

67 https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors_\ 

68of_the_second_derivative 

69 

70 Examples 

71 -------- 

72 >>> import numpy as np 

73 >>> from scipy.sparse.linalg import LaplacianNd 

74 >>> from scipy.sparse import diags, csgraph 

75 >>> from scipy.linalg import eigvalsh 

76 

77 The one-dimensional Laplacian demonstrated below for pure Neumann boundary 

78 conditions on a regular grid with ``n=6`` grid points is exactly the 

79 negative graph Laplacian for the undirected linear graph with ``n`` 

80 vertices using the sparse adjacency matrix ``G`` represented by the 

81 famous tri-diagonal matrix: 

82 

83 >>> n = 6 

84 >>> G = diags(np.ones(n - 1), 1, format='csr') 

85 >>> Lf = csgraph.laplacian(G, symmetrized=True, form='function') 

86 >>> grid_shape = (n, ) 

87 >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann') 

88 >>> np.array_equal(lap.matmat(np.eye(n)), -Lf(np.eye(n))) 

89 True 

90 

91 Since all matrix entries of the Laplacian are integers, ``'int8'`` is 

92 the default dtype for storing matrix representations. 

93 

94 >>> lap.tosparse() 

95 <6x6 sparse array of type '<class 'numpy.int8'>' 

96 with 16 stored elements (3 diagonals) in DIAgonal format> 

97 >>> lap.toarray() 

98 array([[-1, 1, 0, 0, 0, 0], 

99 [ 1, -2, 1, 0, 0, 0], 

100 [ 0, 1, -2, 1, 0, 0], 

101 [ 0, 0, 1, -2, 1, 0], 

102 [ 0, 0, 0, 1, -2, 1], 

103 [ 0, 0, 0, 0, 1, -1]], dtype=int8) 

104 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) 

105 True 

106 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) 

107 True 

108 

109 Any number of extreme eigenvalues and/or eigenvectors can be computed. 

110  

111 >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic') 

112 >>> lap.eigenvalues() 

113 array([-4., -3., -3., -1., -1., 0.]) 

114 >>> lap.eigenvalues()[-2:] 

115 array([-1., 0.]) 

116 >>> lap.eigenvalues(2) 

117 array([-1., 0.]) 

118 >>> lap.eigenvectors(1) 

119 array([[0.40824829], 

120 [0.40824829], 

121 [0.40824829], 

122 [0.40824829], 

123 [0.40824829], 

124 [0.40824829]]) 

125 >>> lap.eigenvectors(2) 

126 array([[ 0.5 , 0.40824829], 

127 [ 0. , 0.40824829], 

128 [-0.5 , 0.40824829], 

129 [-0.5 , 0.40824829], 

130 [ 0. , 0.40824829], 

131 [ 0.5 , 0.40824829]]) 

132 >>> lap.eigenvectors() 

133 array([[ 0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 , 

134 0.40824829], 

135 [-0.40824829, -0.57735027, -0.57735027, 0. , 0. , 

136 0.40824829], 

137 [ 0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 , 

138 0.40824829], 

139 [-0.40824829, 0.28867513, 0.28867513, -0.5 , -0.5 , 

140 0.40824829], 

141 [ 0.40824829, -0.57735027, -0.57735027, 0. , 0. , 

142 0.40824829], 

143 [-0.40824829, 0.28867513, 0.28867513, 0.5 , 0.5 , 

144 0.40824829]]) 

145 

146 The two-dimensional Laplacian is illustrated on a regular grid with 

147 ``grid_shape = (2, 3)`` points in each dimension. 

148 

149 >>> grid_shape = (2, 3) 

150 >>> n = np.prod(grid_shape) 

151 

152 Numeration of grid points is as follows: 

153 

154 >>> np.arange(n).reshape(grid_shape + (-1,)) 

155 array([[[0], 

156 [1], 

157 [2]], 

158 <BLANKLINE> 

159 [[3], 

160 [4], 

161 [5]]]) 

162 

163 Each of the boundary conditions ``'dirichlet'``, ``'periodic'``, and 

164 ``'neumann'`` is illustrated separately; with ``'dirichlet'`` 

165 

166 >>> lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet') 

167 >>> lap.tosparse() 

168 <6x6 sparse array of type '<class 'numpy.int8'>' 

169 with 20 stored elements in Compressed Sparse Row format> 

170 >>> lap.toarray() 

171 array([[-4, 1, 0, 1, 0, 0], 

172 [ 1, -4, 1, 0, 1, 0], 

173 [ 0, 1, -4, 0, 0, 1], 

174 [ 1, 0, 0, -4, 1, 0], 

175 [ 0, 1, 0, 1, -4, 1], 

176 [ 0, 0, 1, 0, 1, -4]], dtype=int8) 

177 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) 

178 True 

179 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) 

180 True 

181 >>> lap.eigenvalues() 

182 array([-6.41421356, -5. , -4.41421356, -3.58578644, -3. , 

183 -1.58578644]) 

184 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) 

185 >>> np.allclose(lap.eigenvalues(), eigvals) 

186 True 

187 >>> np.allclose(lap.toarray() @ lap.eigenvectors(), 

188 ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) 

189 True 

190 

191 with ``'periodic'`` 

192 

193 >>> lap = LaplacianNd(grid_shape, boundary_conditions='periodic') 

194 >>> lap.tosparse() 

195 <6x6 sparse array of type '<class 'numpy.int8'>' 

196 with 24 stored elements in Compressed Sparse Row format> 

197 >>> lap.toarray() 

198 array([[-4, 1, 1, 2, 0, 0], 

199 [ 1, -4, 1, 0, 2, 0], 

200 [ 1, 1, -4, 0, 0, 2], 

201 [ 2, 0, 0, -4, 1, 1], 

202 [ 0, 2, 0, 1, -4, 1], 

203 [ 0, 0, 2, 1, 1, -4]], dtype=int8) 

204 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) 

205 True 

206 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) 

207 True 

208 >>> lap.eigenvalues() 

209 array([-7., -7., -4., -3., -3., 0.]) 

210 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) 

211 >>> np.allclose(lap.eigenvalues(), eigvals) 

212 True 

213 >>> np.allclose(lap.toarray() @ lap.eigenvectors(), 

214 ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) 

215 True 

216 

217 and with ``'neumann'`` 

218 

219 >>> lap = LaplacianNd(grid_shape, boundary_conditions='neumann') 

220 >>> lap.tosparse() 

221 <6x6 sparse array of type '<class 'numpy.int8'>' 

222 with 20 stored elements in Compressed Sparse Row format> 

223 >>> lap.toarray() 

224 array([[-2, 1, 0, 1, 0, 0], 

225 [ 1, -3, 1, 0, 1, 0], 

226 [ 0, 1, -2, 0, 0, 1], 

227 [ 1, 0, 0, -2, 1, 0], 

228 [ 0, 1, 0, 1, -3, 1], 

229 [ 0, 0, 1, 0, 1, -2]]) 

230 >>> np.array_equal(lap.matmat(np.eye(n)), lap.toarray()) 

231 True 

232 >>> np.array_equal(lap.tosparse().toarray(), lap.toarray()) 

233 True 

234 >>> lap.eigenvalues() 

235 array([-5., -3., -3., -2., -1., 0.]) 

236 >>> eigvals = eigvalsh(lap.toarray().astype(np.float64)) 

237 >>> np.allclose(lap.eigenvalues(), eigvals) 

238 True 

239 >>> np.allclose(lap.toarray() @ lap.eigenvectors(), 

240 ... lap.eigenvectors() @ np.diag(lap.eigenvalues())) 

241 True 

242 

243 """ 

244 

245 def __init__(self, grid_shape, *, 

246 boundary_conditions='neumann', 

247 dtype=np.int8): 

248 

249 if boundary_conditions not in ('dirichlet', 'neumann', 'periodic'): 

250 raise ValueError( 

251 f"Unknown value {boundary_conditions!r} is given for " 

252 "'boundary_conditions' parameter. The valid options are " 

253 "'dirichlet', 'periodic', and 'neumann' (default)." 

254 ) 

255 

256 self.grid_shape = grid_shape 

257 self.boundary_conditions = boundary_conditions 

258 # LaplacianNd folds all dimensions in `grid_shape` into a single one 

259 N = np.prod(grid_shape) 

260 super().__init__(dtype=dtype, shape=(N, N)) 

261 

262 def _eigenvalue_ordering(self, m): 

263 """Compute `m` largest eigenvalues in each of the ``N`` directions, 

264 i.e., up to ``m * N`` total, order them and return `m` largest. 

265 """ 

266 grid_shape = self.grid_shape 

267 if m is None: 

268 indices = np.indices(grid_shape) 

269 Leig = np.zeros(grid_shape) 

270 else: 

271 grid_shape_min = min(grid_shape, 

272 tuple(np.ones_like(grid_shape) * m)) 

273 indices = np.indices(grid_shape_min) 

274 Leig = np.zeros(grid_shape_min) 

275 

276 for j, n in zip(indices, grid_shape): 

277 if self.boundary_conditions == 'dirichlet': 

278 Leig += -4 * np.sin(np.pi * (j + 1) / (2 * (n + 1))) ** 2 

279 elif self.boundary_conditions == 'neumann': 

280 Leig += -4 * np.sin(np.pi * j / (2 * n)) ** 2 

281 else: # boundary_conditions == 'periodic' 

282 Leig += -4 * np.sin(np.pi * np.floor((j + 1) / 2) / n) ** 2 

283 

284 Leig_ravel = Leig.ravel() 

285 ind = np.argsort(Leig_ravel) 

286 eigenvalues = Leig_ravel[ind] 

287 if m is not None: 

288 eigenvalues = eigenvalues[-m:] 

289 ind = ind[-m:] 

290 

291 return eigenvalues, ind 

292 

293 def eigenvalues(self, m=None): 

294 """Return the requested number of eigenvalues. 

295  

296 Parameters 

297 ---------- 

298 m : int, optional 

299 The positive number of smallest eigenvalues to return. 

300 If not provided, then all eigenvalues will be returned. 

301  

302 Returns 

303 ------- 

304 eigenvalues : float array 

305 The requested `m` smallest or all eigenvalues, in ascending order. 

306 """ 

307 eigenvalues, _ = self._eigenvalue_ordering(m) 

308 return eigenvalues 

309 

310 def _ev1d(self, j, n): 

311 """Return 1 eigenvector in 1d with index `j` 

312 and number of grid points `n` where ``j < n``.  

313 """ 

314 if self.boundary_conditions == 'dirichlet': 

315 i = np.pi * (np.arange(n) + 1) / (n + 1) 

316 ev = np.sqrt(2. / (n + 1.)) * np.sin(i * (j + 1)) 

317 elif self.boundary_conditions == 'neumann': 

318 i = np.pi * (np.arange(n) + 0.5) / n 

319 ev = np.sqrt((1. if j == 0 else 2.) / n) * np.cos(i * j) 

320 else: # boundary_conditions == 'periodic' 

321 if j == 0: 

322 ev = np.sqrt(1. / n) * np.ones(n) 

323 elif j + 1 == n and n % 2 == 0: 

324 ev = np.sqrt(1. / n) * np.tile([1, -1], n//2) 

325 else: 

326 i = 2. * np.pi * (np.arange(n) + 0.5) / n 

327 ev = np.sqrt(2. / n) * np.cos(i * np.floor((j + 1) / 2)) 

328 # make small values exact zeros correcting round-off errors 

329 # due to symmetry of eigenvectors the exact 0. is correct  

330 ev[np.abs(ev) < np.finfo(np.float64).eps] = 0. 

331 return ev 

332 

333 def _one_eve(self, k): 

334 """Return 1 eigenvector in Nd with multi-index `j` 

335 as a tensor product of the corresponding 1d eigenvectors.  

336 """ 

337 phi = [self._ev1d(j, n) for j, n in zip(k, self.grid_shape)] 

338 result = phi[0] 

339 for phi in phi[1:]: 

340 result = np.tensordot(result, phi, axes=0) 

341 return np.asarray(result).ravel() 

342 

343 def eigenvectors(self, m=None): 

344 """Return the requested number of eigenvectors for ordered eigenvalues. 

345  

346 Parameters 

347 ---------- 

348 m : int, optional 

349 The positive number of eigenvectors to return. If not provided, 

350 then all eigenvectors will be returned. 

351  

352 Returns 

353 ------- 

354 eigenvectors : float array 

355 An array with columns made of the requested `m` or all eigenvectors. 

356 The columns are ordered according to the `m` ordered eigenvalues.  

357 """ 

358 _, ind = self._eigenvalue_ordering(m) 

359 if m is None: 

360 grid_shape_min = self.grid_shape 

361 else: 

362 grid_shape_min = min(self.grid_shape, 

363 tuple(np.ones_like(self.grid_shape) * m)) 

364 

365 N_indices = np.unravel_index(ind, grid_shape_min) 

366 N_indices = [tuple(x) for x in zip(*N_indices)] 

367 eigenvectors_list = [self._one_eve(k) for k in N_indices] 

368 return np.column_stack(eigenvectors_list) 

369 

370 def toarray(self): 

371 """ 

372 Converts the Laplacian data to a dense array. 

373 

374 Returns 

375 ------- 

376 L : ndarray 

377 The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``. 

378 

379 """ 

380 grid_shape = self.grid_shape 

381 n = np.prod(grid_shape) 

382 L = np.zeros([n, n], dtype=np.int8) 

383 # Scratch arrays 

384 L_i = np.empty_like(L) 

385 Ltemp = np.empty_like(L) 

386 

387 for ind, dim in enumerate(grid_shape): 

388 # Start zeroing out L_i 

389 L_i[:] = 0 

390 # Allocate the top left corner with the kernel of L_i 

391 # Einsum returns writable view of arrays 

392 np.einsum("ii->i", L_i[:dim, :dim])[:] = -2 

393 np.einsum("ii->i", L_i[: dim - 1, 1:dim])[:] = 1 

394 np.einsum("ii->i", L_i[1:dim, : dim - 1])[:] = 1 

395 

396 if self.boundary_conditions == 'neumann': 

397 L_i[0, 0] = -1 

398 L_i[dim - 1, dim - 1] = -1 

399 elif self.boundary_conditions == 'periodic': 

400 if dim > 1: 

401 L_i[0, dim - 1] += 1 

402 L_i[dim - 1, 0] += 1 

403 else: 

404 L_i[0, 0] += 1 

405 

406 # kron is too slow for large matrices hence the next two tricks 

407 # 1- kron(eye, mat) is block_diag(mat, mat, ...) 

408 # 2- kron(mat, eye) can be performed by 4d stride trick 

409 

410 # 1- 

411 new_dim = dim 

412 # for block_diag we tile the top left portion on the diagonal 

413 if ind > 0: 

414 tiles = np.prod(grid_shape[:ind]) 

415 for j in range(1, tiles): 

416 L_i[j*dim:(j+1)*dim, j*dim:(j+1)*dim] = L_i[:dim, :dim] 

417 new_dim += dim 

418 # 2- 

419 # we need the keep L_i, but reset the array 

420 Ltemp[:new_dim, :new_dim] = L_i[:new_dim, :new_dim] 

421 tiles = int(np.prod(grid_shape[ind+1:])) 

422 # Zero out the top left, the rest is already 0 

423 L_i[:new_dim, :new_dim] = 0 

424 idx = [x for x in range(tiles)] 

425 L_i.reshape( 

426 (new_dim, tiles, 

427 new_dim, tiles) 

428 )[:, idx, :, idx] = Ltemp[:new_dim, :new_dim] 

429 

430 L += L_i 

431 

432 return L.astype(self.dtype) 

433 

434 def tosparse(self): 

435 """ 

436 Constructs a sparse array from the Laplacian data. The returned sparse 

437 array format is dependent on the selected boundary conditions. 

438 

439 Returns 

440 ------- 

441 L : scipy.sparse.sparray 

442 The shape is ``(N, N)`` where ``N = np.prod(grid_shape)``. 

443 

444 """ 

445 N = len(self.grid_shape) 

446 p = np.prod(self.grid_shape) 

447 L = dia_array((p, p), dtype=np.int8) 

448 

449 for i in range(N): 

450 dim = self.grid_shape[i] 

451 data = np.ones([3, dim], dtype=np.int8) 

452 data[1, :] *= -2 

453 

454 if self.boundary_conditions == 'neumann': 

455 data[1, 0] = -1 

456 data[1, -1] = -1 

457 

458 L_i = dia_array((data, [-1, 0, 1]), shape=(dim, dim), 

459 dtype=np.int8 

460 ) 

461 

462 if self.boundary_conditions == 'periodic': 

463 t = dia_array((dim, dim), dtype=np.int8) 

464 t.setdiag([1], k=-dim+1) 

465 t.setdiag([1], k=dim-1) 

466 L_i += t 

467 

468 for j in range(i): 

469 L_i = kron(eye(self.grid_shape[j], dtype=np.int8), L_i) 

470 for j in range(i + 1, N): 

471 L_i = kron(L_i, eye(self.grid_shape[j], dtype=np.int8)) 

472 L += L_i 

473 return L.astype(self.dtype) 

474 

475 def _matvec(self, x): 

476 grid_shape = self.grid_shape 

477 N = len(grid_shape) 

478 X = x.reshape(grid_shape + (-1,)) 

479 Y = -2 * N * X 

480 for i in range(N): 

481 Y += np.roll(X, 1, axis=i) 

482 Y += np.roll(X, -1, axis=i) 

483 if self.boundary_conditions in ('neumann', 'dirichlet'): 

484 Y[(slice(None),)*i + (0,) + (slice(None),)*(N-i-1) 

485 ] -= np.roll(X, 1, axis=i)[ 

486 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1) 

487 ] 

488 Y[ 

489 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1) 

490 ] -= np.roll(X, -1, axis=i)[ 

491 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1) 

492 ] 

493 

494 if self.boundary_conditions == 'neumann': 

495 Y[ 

496 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1) 

497 ] += np.roll(X, 0, axis=i)[ 

498 (slice(None),) * i + (0,) + (slice(None),) * (N-i-1) 

499 ] 

500 Y[ 

501 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1) 

502 ] += np.roll(X, 0, axis=i)[ 

503 (slice(None),) * i + (-1,) + (slice(None),) * (N-i-1) 

504 ] 

505 

506 return Y.reshape(-1, X.shape[-1]) 

507 

508 def _matmat(self, x): 

509 return self._matvec(x) 

510 

511 def _adjoint(self): 

512 return self 

513 

514 def _transpose(self): 

515 return self 

516 

517 

518class Sakurai(LinearOperator): 

519 """ 

520 Construct a Sakurai matrix in various formats and its eigenvalues. 

521 

522 Constructs the "Sakurai" matrix motivated by reference [1]_: 

523 square real symmetric positive definite and 5-diagonal 

524 with the main digonal ``[5, 6, 6, ..., 6, 6, 5], the ``+1`` and ``-1`` 

525 diagonals filled with ``-4``, and the ``+2`` and ``-2`` diagonals 

526 made of ``1``. Its eigenvalues are analytically known to be 

527 ``16. * np.power(np.cos(0.5 * k * np.pi / (n + 1)), 4)``. 

528 The matrix gets ill-conditioned with its size growing. 

529 It is useful for testing and benchmarking sparse eigenvalue solvers 

530 especially those taking advantage of its banded 5-diagonal structure. 

531 See the notes below for details. 

532 

533 Parameters 

534 ---------- 

535 n : int 

536 The size of the matrix. 

537 dtype : dtype 

538 Numerical type of the array. Default is ``np.int8``. 

539 

540 Methods 

541 ------- 

542 toarray() 

543 Construct a dense array from Laplacian data 

544 tosparse() 

545 Construct a sparse array from Laplacian data 

546 tobanded() 

547 The Sakurai matrix in the format for banded symmetric matrices, 

548 i.e., (3, n) ndarray with 3 upper diagonals 

549 placing the main diagonal at the bottom. 

550 eigenvalues 

551 All eigenvalues of the Sakurai matrix ordered ascending. 

552 

553 Notes 

554 ----- 

555 Reference [1]_ introduces a generalized eigenproblem for the matrix pair 

556 `A` and `B` where `A` is the identity so we turn it into an eigenproblem 

557 just for the matrix `B` that this function outputs in various formats 

558 together with its eigenvalues. 

559  

560 .. versionadded:: 1.12.0 

561 

562 References 

563 ---------- 

564 .. [1] T. Sakurai, H. Tadano, Y. Inadomi, and U. Nagashima, 

565 "A moment-based method for large-scale generalized 

566 eigenvalue problems", 

567 Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004). 

568 

569 Examples 

570 -------- 

571 >>> import numpy as np 

572 >>> from scipy.sparse.linalg._special_sparse_arrays import Sakurai 

573 >>> from scipy.linalg import eig_banded 

574 >>> n = 6 

575 >>> sak = Sakurai(n) 

576 

577 Since all matrix entries are small integers, ``'int8'`` is 

578 the default dtype for storing matrix representations. 

579 

580 >>> sak.toarray() 

581 array([[ 5, -4, 1, 0, 0, 0], 

582 [-4, 6, -4, 1, 0, 0], 

583 [ 1, -4, 6, -4, 1, 0], 

584 [ 0, 1, -4, 6, -4, 1], 

585 [ 0, 0, 1, -4, 6, -4], 

586 [ 0, 0, 0, 1, -4, 5]], dtype=int8) 

587 >>> sak.tobanded() 

588 array([[ 1, 1, 1, 1, 1, 1], 

589 [-4, -4, -4, -4, -4, -4], 

590 [ 5, 6, 6, 6, 6, 5]], dtype=int8) 

591 >>> sak.tosparse() 

592 <6x6 sparse matrix of type '<class 'numpy.int8'>' 

593 with 24 stored elements (5 diagonals) in DIAgonal format> 

594 >>> np.array_equal(sak.dot(np.eye(n)), sak.tosparse().toarray()) 

595 True 

596 >>> sak.eigenvalues() 

597 array([0.03922866, 0.56703972, 2.41789479, 5.97822974, 

598 10.54287655, 14.45473055]) 

599 >>> sak.eigenvalues(2) 

600 array([0.03922866, 0.56703972]) 

601 

602 The banded form can be used in scipy functions for banded matrices, e.g., 

603 

604 >>> e = eig_banded(sak.tobanded(), eigvals_only=True) 

605 >>> np.allclose(sak.eigenvalues, e, atol= n * n * n * np.finfo(float).eps) 

606 True 

607 

608 """ 

609 def __init__(self, n, dtype=np.int8): 

610 self.n = n 

611 self.dtype = dtype 

612 shape = (n, n) 

613 super().__init__(dtype, shape) 

614 

615 def eigenvalues(self, m=None): 

616 """Return the requested number of eigenvalues. 

617  

618 Parameters 

619 ---------- 

620 m : int, optional 

621 The positive number of smallest eigenvalues to return. 

622 If not provided, then all eigenvalues will be returned. 

623  

624 Returns 

625 ------- 

626 eigenvalues : `np.float64` array 

627 The requested `m` smallest or all eigenvalues, in ascending order. 

628 """ 

629 if m is None: 

630 m = self.n 

631 k = np.arange(self.n + 1 -m, self.n + 1) 

632 return np.flip(16. * np.power(np.cos(0.5 * k * np.pi / (self.n + 1)), 4)) 

633 

634 def tobanded(self): 

635 """ 

636 Construct the Sakurai matrix as a banded array. 

637 """ 

638 d0 = np.r_[5, 6 * np.ones(self.n - 2, dtype=self.dtype), 5] 

639 d1 = -4 * np.ones(self.n, dtype=self.dtype) 

640 d2 = np.ones(self.n, dtype=self.dtype) 

641 return np.array([d2, d1, d0]).astype(self.dtype) 

642 

643 def tosparse(self): 

644 """ 

645 Construct the Sakurai matrix is a sparse format. 

646 """ 

647 from scipy.sparse import spdiags 

648 d = self.tobanded() 

649 # the banded format has the main diagonal at the bottom 

650 # `spdiags` has no `dtype` parameter so inherits dtype from banded 

651 return spdiags([d[0], d[1], d[2], d[1], d[0]], [-2, -1, 0, 1, 2], 

652 self.n, self.n) 

653 

654 def toarray(self): 

655 return self.tosparse().toarray() 

656 

657 def _matvec(self, x): 

658 """ 

659 Construct matrix-free callable banded-matrix-vector multiplication by 

660 the Sakurai matrix without constructing or storing the matrix itself 

661 using the knowledge of its entries and the 5-diagonal format. 

662 """ 

663 x = x.reshape(self.n, -1) 

664 result_dtype = np.promote_types(x.dtype, self.dtype) 

665 sx = np.zeros_like(x, dtype=result_dtype) 

666 sx[0, :] = 5 * x[0, :] - 4 * x[1, :] + x[2, :] 

667 sx[-1, :] = 5 * x[-1, :] - 4 * x[-2, :] + x[-3, :] 

668 sx[1: -1, :] = (6 * x[1: -1, :] - 4 * (x[:-2, :] + x[2:, :]) 

669 + np.pad(x[:-3, :], ((1, 0), (0, 0))) 

670 + np.pad(x[3:, :], ((0, 1), (0, 0)))) 

671 return sx 

672 

673 def _matmat(self, x): 

674 """ 

675 Construct matrix-free callable matrix-matrix multiplication by 

676 the Sakurai matrix without constructing or storing the matrix itself 

677 by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. 

678 """ 

679 return self._matvec(x) 

680 

681 def _adjoint(self): 

682 return self 

683 

684 def _transpose(self): 

685 return self 

686 

687 

688class MikotaM(LinearOperator): 

689 """ 

690 Construct a mass matrix in various formats of Mikota pair. 

691 

692 The mass matrix `M` is square real diagonal 

693 positive definite with entries that are reciprocal to integers. 

694 

695 Parameters 

696 ---------- 

697 shape : tuple of int 

698 The shape of the matrix. 

699 dtype : dtype 

700 Numerical type of the array. Default is ``np.float64``. 

701 

702 Methods 

703 ------- 

704 toarray() 

705 Construct a dense array from Mikota data 

706 tosparse() 

707 Construct a sparse array from Mikota data 

708 tobanded() 

709 The format for banded symmetric matrices, 

710 i.e., (1, n) ndarray with the main diagonal. 

711 """ 

712 def __init__(self, shape, dtype=np.float64): 

713 self.shape = shape 

714 self.dtype = dtype 

715 super().__init__(dtype, shape) 

716 

717 def _diag(self): 

718 # The matrix is constructed from its diagonal 1 / [1, ..., N+1]; 

719 # compute in a function to avoid duplicated code & storage footprint 

720 return (1. / np.arange(1, self.shape[0] + 1)).astype(self.dtype) 

721 

722 def tobanded(self): 

723 return self._diag() 

724 

725 def tosparse(self): 

726 from scipy.sparse import diags 

727 return diags([self._diag()], [0], shape=self.shape, dtype=self.dtype) 

728 

729 def toarray(self): 

730 return np.diag(self._diag()).astype(self.dtype) 

731 

732 def _matvec(self, x): 

733 """ 

734 Construct matrix-free callable banded-matrix-vector multiplication by 

735 the Mikota mass matrix without constructing or storing the matrix itself 

736 using the knowledge of its entries and the diagonal format. 

737 """ 

738 x = x.reshape(self.shape[0], -1) 

739 return self._diag()[:, np.newaxis] * x 

740 

741 def _matmat(self, x): 

742 """ 

743 Construct matrix-free callable matrix-matrix multiplication by 

744 the Mikota mass matrix without constructing or storing the matrix itself 

745 by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. 

746 """ 

747 return self._matvec(x) 

748 

749 def _adjoint(self): 

750 return self 

751 

752 def _transpose(self): 

753 return self 

754 

755 

756class MikotaK(LinearOperator): 

757 """ 

758 Construct a stiffness matrix in various formats of Mikota pair. 

759 

760 The stiffness matrix `K` is square real tri-diagonal symmetric 

761 positive definite with integer entries.  

762 

763 Parameters 

764 ---------- 

765 shape : tuple of int 

766 The shape of the matrix. 

767 dtype : dtype 

768 Numerical type of the array. Default is ``np.int32``. 

769 

770 Methods 

771 ------- 

772 toarray() 

773 Construct a dense array from Mikota data 

774 tosparse() 

775 Construct a sparse array from Mikota data 

776 tobanded() 

777 The format for banded symmetric matrices, 

778 i.e., (2, n) ndarray with 2 upper diagonals 

779 placing the main diagonal at the bottom. 

780 """ 

781 def __init__(self, shape, dtype=np.int32): 

782 self.shape = shape 

783 self.dtype = dtype 

784 super().__init__(dtype, shape) 

785 # The matrix is constructed from its diagonals; 

786 # we precompute these to avoid duplicating the computation 

787 n = shape[0] 

788 self._diag0 = np.arange(2 * n - 1, 0, -2, dtype=self.dtype) 

789 self._diag1 = - np.arange(n - 1, 0, -1, dtype=self.dtype) 

790 

791 def tobanded(self): 

792 return np.array([np.pad(self._diag1, (1, 0), 'constant'), self._diag0]) 

793 

794 def tosparse(self): 

795 from scipy.sparse import diags 

796 return diags([self._diag1, self._diag0, self._diag1], [-1, 0, 1], 

797 shape=self.shape, dtype=self.dtype) 

798 

799 def toarray(self): 

800 return self.tosparse().toarray() 

801 

802 def _matvec(self, x): 

803 """ 

804 Construct matrix-free callable banded-matrix-vector multiplication by 

805 the Mikota stiffness matrix without constructing or storing the matrix 

806 itself using the knowledge of its entries and the 3-diagonal format. 

807 """ 

808 x = x.reshape(self.shape[0], -1) 

809 result_dtype = np.promote_types(x.dtype, self.dtype) 

810 kx = np.zeros_like(x, dtype=result_dtype) 

811 d1 = self._diag1 

812 d0 = self._diag0 

813 kx[0, :] = d0[0] * x[0, :] + d1[0] * x[1, :] 

814 kx[-1, :] = d1[-1] * x[-2, :] + d0[-1] * x[-1, :] 

815 kx[1: -1, :] = (d1[:-1, None] * x[: -2, :] 

816 + d0[1: -1, None] * x[1: -1, :] 

817 + d1[1:, None] * x[2:, :]) 

818 return kx 

819 

820 def _matmat(self, x): 

821 """ 

822 Construct matrix-free callable matrix-matrix multiplication by 

823 the Stiffness mass matrix without constructing or storing the matrix itself 

824 by reusing the ``_matvec(x)`` that supports both 1D and 2D arrays ``x``. 

825 """ 

826 return self._matvec(x) 

827 

828 def _adjoint(self): 

829 return self 

830 

831 def _transpose(self): 

832 return self 

833 

834 

835class MikotaPair: 

836 """ 

837 Construct the Mikota pair of matrices in various formats and 

838 eigenvalues of the generalized eigenproblem with them. 

839 

840 The Mikota pair of matrices [1, 2]_ models a vibration problem 

841 of a linear mass-spring system with the ends attached where 

842 the stiffness of the springs and the masses increase along 

843 the system length such that vibration frequencies are subsequent 

844 integers 1, 2, ..., `n` where `n` is the number of the masses. Thus, 

845 eigenvalues of the generalized eigenvalue problem for 

846 the matrix pair `K` and `M` where `K` is he system stiffness matrix 

847 and `M` is the system mass matrix are the squares of the integers, 

848 i.e., 1, 4, 9, ..., ``n * n``. 

849 

850 The stiffness matrix `K` is square real tri-diagonal symmetric 

851 positive definite. The mass matrix `M` is diagonal with diagonal 

852 entries 1, 1/2, 1/3, ...., ``1/n``. Both matrices get 

853 ill-conditioned with `n` growing. 

854 

855 Parameters 

856 ---------- 

857 n : int 

858 The size of the matrices of the Mikota pair. 

859 dtype : dtype 

860 Numerical type of the array. Default is ``np.float64``. 

861 

862 Attributes 

863 ---------- 

864 eigenvalues : 1D ndarray, ``np.uint64`` 

865 All eigenvalues of the Mikota pair ordered ascending. 

866 

867 Methods 

868 ------- 

869 MikotaK() 

870 A `LinearOperator` custom object for the stiffness matrix. 

871 MikotaM() 

872 A `LinearOperator` custom object for the mass matrix. 

873  

874 .. versionadded:: 1.12.0 

875 

876 References 

877 ---------- 

878 .. [1] J. Mikota, "Frequency tuning of chain structure multibody oscillators 

879 to place the natural frequencies at omega1 and N-1 integer multiples 

880 omega2,..., omegaN", Z. Angew. Math. Mech. 81 (2001), S2, S201-S202. 

881 Appl. Num. Anal. Comp. Math. Vol. 1 No. 2 (2004). 

882 .. [2] Peter C. Muller and Metin Gurgoze, 

883 "Natural frequencies of a multi-degree-of-freedom vibration system", 

884 Proc. Appl. Math. Mech. 6, 319-320 (2006). 

885 http://dx.doi.org/10.1002/pamm.200610141.