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

1"""Abstract linear algebra library. 

2 

3This module defines a class hierarchy that implements a kind of "lazy" 

4matrix representation, called the ``LinearOperator``. It can be used to do 

5linear algebra with extremely large sparse or structured matrices, without 

6representing those explicitly in memory. Such matrices can be added, 

7multiplied, transposed, etc. 

8 

9As a motivating example, suppose you want have a matrix where almost all of 

10the elements have the value one. The standard sparse matrix representation 

11skips the storage of zeros, but not ones. By contrast, a LinearOperator is 

12able to represent such matrices efficiently. First, we need a compact way to 

13represent an all-ones matrix:: 

14 

15 >>> import numpy as np 

16 >>> from scipy.sparse.linalg._interface import LinearOperator 

17 >>> class Ones(LinearOperator): 

18 ... def __init__(self, shape): 

19 ... super().__init__(dtype=None, shape=shape) 

20 ... def _matvec(self, x): 

21 ... return np.repeat(x.sum(), self.shape[0]) 

22 

23Instances of this class emulate ``np.ones(shape)``, but using a constant 

24amount of storage, independent of ``shape``. The ``_matvec`` method specifies 

25how this linear operator multiplies with (operates on) a vector. We can now 

26add this operator to a sparse matrix that stores only offsets from one:: 

27 

28 >>> from scipy.sparse.linalg._interface import aslinearoperator 

29 >>> from scipy.sparse import csr_matrix 

30 >>> offsets = csr_matrix([[1, 0, 2], [0, -1, 0], [0, 0, 3]]) 

31 >>> A = aslinearoperator(offsets) + Ones(offsets.shape) 

32 >>> A.dot([1, 2, 3]) 

33 array([13, 4, 15]) 

34 

35The result is the same as that given by its dense, explicitly-stored 

36counterpart:: 

37 

38 >>> (np.ones(A.shape, A.dtype) + offsets.toarray()).dot([1, 2, 3]) 

39 array([13, 4, 15]) 

40 

41Several algorithms in the ``scipy.sparse`` library are able to operate on 

42``LinearOperator`` instances. 

43""" 

44 

45import warnings 

46 

47import numpy as np 

48 

49from scipy.sparse import issparse 

50from scipy.sparse._sputils import isshape, isintlike, asmatrix, is_pydata_spmatrix 

51 

52__all__ = ['LinearOperator', 'aslinearoperator'] 

53 

54 

55class LinearOperator: 

56 """Common interface for performing matrix vector products 

57 

58 Many iterative methods (e.g. cg, gmres) do not need to know the 

59 individual entries of a matrix to solve a linear system A*x=b. 

60 Such solvers only require the computation of matrix vector 

61 products, A*v where v is a dense vector. This class serves as 

62 an abstract interface between iterative solvers and matrix-like 

63 objects. 

64 

65 To construct a concrete LinearOperator, either pass appropriate 

66 callables to the constructor of this class, or subclass it. 

67 

68 A subclass must implement either one of the methods ``_matvec`` 

69 and ``_matmat``, and the attributes/properties ``shape`` (pair of 

70 integers) and ``dtype`` (may be None). It may call the ``__init__`` 

71 on this class to have these attributes validated. Implementing 

72 ``_matvec`` automatically implements ``_matmat`` (using a naive 

73 algorithm) and vice-versa. 

74 

75 Optionally, a subclass may implement ``_rmatvec`` or ``_adjoint`` 

76 to implement the Hermitian adjoint (conjugate transpose). As with 

77 ``_matvec`` and ``_matmat``, implementing either ``_rmatvec`` or 

78 ``_adjoint`` implements the other automatically. Implementing 

79 ``_adjoint`` is preferable; ``_rmatvec`` is mostly there for 

80 backwards compatibility. 

81 

82 Parameters 

83 ---------- 

84 shape : tuple 

85 Matrix dimensions (M, N). 

86 matvec : callable f(v) 

87 Returns returns A * v. 

88 rmatvec : callable f(v) 

89 Returns A^H * v, where A^H is the conjugate transpose of A. 

90 matmat : callable f(V) 

91 Returns A * V, where V is a dense matrix with dimensions (N, K). 

92 dtype : dtype 

93 Data type of the matrix. 

94 rmatmat : callable f(V) 

95 Returns A^H * V, where V is a dense matrix with dimensions (M, K). 

96 

97 Attributes 

98 ---------- 

99 args : tuple 

100 For linear operators describing products etc. of other linear 

101 operators, the operands of the binary operation. 

102 ndim : int 

103 Number of dimensions (this is always 2) 

104 

105 See Also 

106 -------- 

107 aslinearoperator : Construct LinearOperators 

108 

109 Notes 

110 ----- 

111 The user-defined matvec() function must properly handle the case 

112 where v has shape (N,) as well as the (N,1) case. The shape of 

113 the return type is handled internally by LinearOperator. 

114 

115 LinearOperator instances can also be multiplied, added with each 

116 other and exponentiated, all lazily: the result of these operations 

117 is always a new, composite LinearOperator, that defers linear 

118 operations to the original operators and combines the results. 

119 

120 More details regarding how to subclass a LinearOperator and several 

121 examples of concrete LinearOperator instances can be found in the 

122 external project `PyLops <https://pylops.readthedocs.io>`_. 

123 

124 

125 Examples 

126 -------- 

127 >>> import numpy as np 

128 >>> from scipy.sparse.linalg import LinearOperator 

129 >>> def mv(v): 

130 ... return np.array([2*v[0], 3*v[1]]) 

131 ... 

132 >>> A = LinearOperator((2,2), matvec=mv) 

133 >>> A 

134 <2x2 _CustomLinearOperator with dtype=float64> 

135 >>> A.matvec(np.ones(2)) 

136 array([ 2., 3.]) 

137 >>> A * np.ones(2) 

138 array([ 2., 3.]) 

139 

140 """ 

141 

142 ndim = 2 

143 # Necessary for right matmul with numpy arrays. 

144 __array_ufunc__ = None 

145 

146 def __new__(cls, *args, **kwargs): 

147 if cls is LinearOperator: 

148 # Operate as _CustomLinearOperator factory. 

149 return super().__new__(_CustomLinearOperator) 

150 else: 

151 obj = super().__new__(cls) 

152 

153 if (type(obj)._matvec == LinearOperator._matvec 

154 and type(obj)._matmat == LinearOperator._matmat): 

155 warnings.warn("LinearOperator subclass should implement" 

156 " at least one of _matvec and _matmat.", 

157 category=RuntimeWarning, stacklevel=2) 

158 

159 return obj 

160 

161 def __init__(self, dtype, shape): 

162 """Initialize this LinearOperator. 

163 

164 To be called by subclasses. ``dtype`` may be None; ``shape`` should 

165 be convertible to a length-2 tuple. 

166 """ 

167 if dtype is not None: 

168 dtype = np.dtype(dtype) 

169 

170 shape = tuple(shape) 

171 if not isshape(shape): 

172 raise ValueError(f"invalid shape {shape!r} (must be 2-d)") 

173 

174 self.dtype = dtype 

175 self.shape = shape 

176 

177 def _init_dtype(self): 

178 """Called from subclasses at the end of the __init__ routine. 

179 """ 

180 if self.dtype is None: 

181 v = np.zeros(self.shape[-1]) 

182 self.dtype = np.asarray(self.matvec(v)).dtype 

183 

184 def _matmat(self, X): 

185 """Default matrix-matrix multiplication handler. 

186 

187 Falls back on the user-defined _matvec method, so defining that will 

188 define matrix multiplication (though in a very suboptimal way). 

189 """ 

190 

191 return np.hstack([self.matvec(col.reshape(-1,1)) for col in X.T]) 

192 

193 def _matvec(self, x): 

194 """Default matrix-vector multiplication handler. 

195 

196 If self is a linear operator of shape (M, N), then this method will 

197 be called on a shape (N,) or (N, 1) ndarray, and should return a 

198 shape (M,) or (M, 1) ndarray. 

199 

200 This default implementation falls back on _matmat, so defining that 

201 will define matrix-vector multiplication as well. 

202 """ 

203 return self.matmat(x.reshape(-1, 1)) 

204 

205 def matvec(self, x): 

206 """Matrix-vector multiplication. 

207 

208 Performs the operation y=A*x where A is an MxN linear 

209 operator and x is a column vector or 1-d array. 

210 

211 Parameters 

212 ---------- 

213 x : {matrix, ndarray} 

214 An array with shape (N,) or (N,1). 

215 

216 Returns 

217 ------- 

218 y : {matrix, ndarray} 

219 A matrix or ndarray with shape (M,) or (M,1) depending 

220 on the type and shape of the x argument. 

221 

222 Notes 

223 ----- 

224 This matvec wraps the user-specified matvec routine or overridden 

225 _matvec method to ensure that y has the correct shape and type. 

226 

227 """ 

228 

229 x = np.asanyarray(x) 

230 

231 M,N = self.shape 

232 

233 if x.shape != (N,) and x.shape != (N,1): 

234 raise ValueError('dimension mismatch') 

235 

236 y = self._matvec(x) 

237 

238 if isinstance(x, np.matrix): 

239 y = asmatrix(y) 

240 else: 

241 y = np.asarray(y) 

242 

243 if x.ndim == 1: 

244 y = y.reshape(M) 

245 elif x.ndim == 2: 

246 y = y.reshape(M,1) 

247 else: 

248 raise ValueError('invalid shape returned by user-defined matvec()') 

249 

250 return y 

251 

252 def rmatvec(self, x): 

253 """Adjoint matrix-vector multiplication. 

254 

255 Performs the operation y = A^H * x where A is an MxN linear 

256 operator and x is a column vector or 1-d array. 

257 

258 Parameters 

259 ---------- 

260 x : {matrix, ndarray} 

261 An array with shape (M,) or (M,1). 

262 

263 Returns 

264 ------- 

265 y : {matrix, ndarray} 

266 A matrix or ndarray with shape (N,) or (N,1) depending 

267 on the type and shape of the x argument. 

268 

269 Notes 

270 ----- 

271 This rmatvec wraps the user-specified rmatvec routine or overridden 

272 _rmatvec method to ensure that y has the correct shape and type. 

273 

274 """ 

275 

276 x = np.asanyarray(x) 

277 

278 M,N = self.shape 

279 

280 if x.shape != (M,) and x.shape != (M,1): 

281 raise ValueError('dimension mismatch') 

282 

283 y = self._rmatvec(x) 

284 

285 if isinstance(x, np.matrix): 

286 y = asmatrix(y) 

287 else: 

288 y = np.asarray(y) 

289 

290 if x.ndim == 1: 

291 y = y.reshape(N) 

292 elif x.ndim == 2: 

293 y = y.reshape(N,1) 

294 else: 

295 raise ValueError('invalid shape returned by user-defined rmatvec()') 

296 

297 return y 

298 

299 def _rmatvec(self, x): 

300 """Default implementation of _rmatvec; defers to adjoint.""" 

301 if type(self)._adjoint == LinearOperator._adjoint: 

302 # _adjoint not overridden, prevent infinite recursion 

303 raise NotImplementedError 

304 else: 

305 return self.H.matvec(x) 

306 

307 def matmat(self, X): 

308 """Matrix-matrix multiplication. 

309 

310 Performs the operation y=A*X where A is an MxN linear 

311 operator and X dense N*K matrix or ndarray. 

312 

313 Parameters 

314 ---------- 

315 X : {matrix, ndarray} 

316 An array with shape (N,K). 

317 

318 Returns 

319 ------- 

320 Y : {matrix, ndarray} 

321 A matrix or ndarray with shape (M,K) depending on 

322 the type of the X argument. 

323 

324 Notes 

325 ----- 

326 This matmat wraps any user-specified matmat routine or overridden 

327 _matmat method to ensure that y has the correct type. 

328 

329 """ 

330 if not (issparse(X) or is_pydata_spmatrix(X)): 

331 X = np.asanyarray(X) 

332 

333 if X.ndim != 2: 

334 raise ValueError(f'expected 2-d ndarray or matrix, not {X.ndim}-d') 

335 

336 if X.shape[0] != self.shape[1]: 

337 raise ValueError(f'dimension mismatch: {self.shape}, {X.shape}') 

338 

339 try: 

340 Y = self._matmat(X) 

341 except Exception as e: 

342 if issparse(X) or is_pydata_spmatrix(X): 

343 raise TypeError( 

344 "Unable to multiply a LinearOperator with a sparse matrix." 

345 " Wrap the matrix in aslinearoperator first." 

346 ) from e 

347 raise 

348 

349 if isinstance(Y, np.matrix): 

350 Y = asmatrix(Y) 

351 

352 return Y 

353 

354 def rmatmat(self, X): 

355 """Adjoint matrix-matrix multiplication. 

356 

357 Performs the operation y = A^H * x where A is an MxN linear 

358 operator and x is a column vector or 1-d array, or 2-d array. 

359 The default implementation defers to the adjoint. 

360 

361 Parameters 

362 ---------- 

363 X : {matrix, ndarray} 

364 A matrix or 2D array. 

365 

366 Returns 

367 ------- 

368 Y : {matrix, ndarray} 

369 A matrix or 2D array depending on the type of the input. 

370 

371 Notes 

372 ----- 

373 This rmatmat wraps the user-specified rmatmat routine. 

374 

375 """ 

376 if not (issparse(X) or is_pydata_spmatrix(X)): 

377 X = np.asanyarray(X) 

378 

379 if X.ndim != 2: 

380 raise ValueError('expected 2-d ndarray or matrix, not %d-d' 

381 % X.ndim) 

382 

383 if X.shape[0] != self.shape[0]: 

384 raise ValueError(f'dimension mismatch: {self.shape}, {X.shape}') 

385 

386 try: 

387 Y = self._rmatmat(X) 

388 except Exception as e: 

389 if issparse(X) or is_pydata_spmatrix(X): 

390 raise TypeError( 

391 "Unable to multiply a LinearOperator with a sparse matrix." 

392 " Wrap the matrix in aslinearoperator() first." 

393 ) from e 

394 raise 

395 

396 if isinstance(Y, np.matrix): 

397 Y = asmatrix(Y) 

398 return Y 

399 

400 def _rmatmat(self, X): 

401 """Default implementation of _rmatmat defers to rmatvec or adjoint.""" 

402 if type(self)._adjoint == LinearOperator._adjoint: 

403 return np.hstack([self.rmatvec(col.reshape(-1, 1)) for col in X.T]) 

404 else: 

405 return self.H.matmat(X) 

406 

407 def __call__(self, x): 

408 return self*x 

409 

410 def __mul__(self, x): 

411 return self.dot(x) 

412 

413 def __truediv__(self, other): 

414 if not np.isscalar(other): 

415 raise ValueError("Can only divide a linear operator by a scalar.") 

416 

417 return _ScaledLinearOperator(self, 1.0/other) 

418 

419 def dot(self, x): 

420 """Matrix-matrix or matrix-vector multiplication. 

421 

422 Parameters 

423 ---------- 

424 x : array_like 

425 1-d or 2-d array, representing a vector or matrix. 

426 

427 Returns 

428 ------- 

429 Ax : array 

430 1-d or 2-d array (depending on the shape of x) that represents 

431 the result of applying this linear operator on x. 

432 

433 """ 

434 if isinstance(x, LinearOperator): 

435 return _ProductLinearOperator(self, x) 

436 elif np.isscalar(x): 

437 return _ScaledLinearOperator(self, x) 

438 else: 

439 if not issparse(x) and not is_pydata_spmatrix(x): 

440 # Sparse matrices shouldn't be converted to numpy arrays. 

441 x = np.asarray(x) 

442 

443 if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1: 

444 return self.matvec(x) 

445 elif x.ndim == 2: 

446 return self.matmat(x) 

447 else: 

448 raise ValueError('expected 1-d or 2-d array or matrix, got %r' 

449 % x) 

450 

451 def __matmul__(self, other): 

452 if np.isscalar(other): 

453 raise ValueError("Scalar operands are not allowed, " 

454 "use '*' instead") 

455 return self.__mul__(other) 

456 

457 def __rmatmul__(self, other): 

458 if np.isscalar(other): 

459 raise ValueError("Scalar operands are not allowed, " 

460 "use '*' instead") 

461 return self.__rmul__(other) 

462 

463 def __rmul__(self, x): 

464 if np.isscalar(x): 

465 return _ScaledLinearOperator(self, x) 

466 else: 

467 return self._rdot(x) 

468 

469 def _rdot(self, x): 

470 """Matrix-matrix or matrix-vector multiplication from the right. 

471 

472 Parameters 

473 ---------- 

474 x : array_like 

475 1-d or 2-d array, representing a vector or matrix. 

476 

477 Returns 

478 ------- 

479 xA : array 

480 1-d or 2-d array (depending on the shape of x) that represents 

481 the result of applying this linear operator on x from the right. 

482 

483 Notes 

484 ----- 

485 This is copied from dot to implement right multiplication. 

486 """ 

487 if isinstance(x, LinearOperator): 

488 return _ProductLinearOperator(x, self) 

489 elif np.isscalar(x): 

490 return _ScaledLinearOperator(self, x) 

491 else: 

492 if not issparse(x) and not is_pydata_spmatrix(x): 

493 # Sparse matrices shouldn't be converted to numpy arrays. 

494 x = np.asarray(x) 

495 

496 # We use transpose instead of rmatvec/rmatmat to avoid 

497 # unnecessary complex conjugation if possible. 

498 if x.ndim == 1 or x.ndim == 2 and x.shape[0] == 1: 

499 return self.T.matvec(x.T).T 

500 elif x.ndim == 2: 

501 return self.T.matmat(x.T).T 

502 else: 

503 raise ValueError('expected 1-d or 2-d array or matrix, got %r' 

504 % x) 

505 

506 def __pow__(self, p): 

507 if np.isscalar(p): 

508 return _PowerLinearOperator(self, p) 

509 else: 

510 return NotImplemented 

511 

512 def __add__(self, x): 

513 if isinstance(x, LinearOperator): 

514 return _SumLinearOperator(self, x) 

515 else: 

516 return NotImplemented 

517 

518 def __neg__(self): 

519 return _ScaledLinearOperator(self, -1) 

520 

521 def __sub__(self, x): 

522 return self.__add__(-x) 

523 

524 def __repr__(self): 

525 M,N = self.shape 

526 if self.dtype is None: 

527 dt = 'unspecified dtype' 

528 else: 

529 dt = 'dtype=' + str(self.dtype) 

530 

531 return '<%dx%d %s with %s>' % (M, N, self.__class__.__name__, dt) 

532 

533 def adjoint(self): 

534 """Hermitian adjoint. 

535 

536 Returns the Hermitian adjoint of self, aka the Hermitian 

537 conjugate or Hermitian transpose. For a complex matrix, the 

538 Hermitian adjoint is equal to the conjugate transpose. 

539 

540 Can be abbreviated self.H instead of self.adjoint(). 

541 

542 Returns 

543 ------- 

544 A_H : LinearOperator 

545 Hermitian adjoint of self. 

546 """ 

547 return self._adjoint() 

548 

549 H = property(adjoint) 

550 

551 def transpose(self): 

552 """Transpose this linear operator. 

553 

554 Returns a LinearOperator that represents the transpose of this one. 

555 Can be abbreviated self.T instead of self.transpose(). 

556 """ 

557 return self._transpose() 

558 

559 T = property(transpose) 

560 

561 def _adjoint(self): 

562 """Default implementation of _adjoint; defers to rmatvec.""" 

563 return _AdjointLinearOperator(self) 

564 

565 def _transpose(self): 

566 """ Default implementation of _transpose; defers to rmatvec + conj""" 

567 return _TransposedLinearOperator(self) 

568 

569 

570class _CustomLinearOperator(LinearOperator): 

571 """Linear operator defined in terms of user-specified operations.""" 

572 

573 def __init__(self, shape, matvec, rmatvec=None, matmat=None, 

574 dtype=None, rmatmat=None): 

575 super().__init__(dtype, shape) 

576 

577 self.args = () 

578 

579 self.__matvec_impl = matvec 

580 self.__rmatvec_impl = rmatvec 

581 self.__rmatmat_impl = rmatmat 

582 self.__matmat_impl = matmat 

583 

584 self._init_dtype() 

585 

586 def _matmat(self, X): 

587 if self.__matmat_impl is not None: 

588 return self.__matmat_impl(X) 

589 else: 

590 return super()._matmat(X) 

591 

592 def _matvec(self, x): 

593 return self.__matvec_impl(x) 

594 

595 def _rmatvec(self, x): 

596 func = self.__rmatvec_impl 

597 if func is None: 

598 raise NotImplementedError("rmatvec is not defined") 

599 return self.__rmatvec_impl(x) 

600 

601 def _rmatmat(self, X): 

602 if self.__rmatmat_impl is not None: 

603 return self.__rmatmat_impl(X) 

604 else: 

605 return super()._rmatmat(X) 

606 

607 def _adjoint(self): 

608 return _CustomLinearOperator(shape=(self.shape[1], self.shape[0]), 

609 matvec=self.__rmatvec_impl, 

610 rmatvec=self.__matvec_impl, 

611 matmat=self.__rmatmat_impl, 

612 rmatmat=self.__matmat_impl, 

613 dtype=self.dtype) 

614 

615 

616class _AdjointLinearOperator(LinearOperator): 

617 """Adjoint of arbitrary Linear Operator""" 

618 

619 def __init__(self, A): 

620 shape = (A.shape[1], A.shape[0]) 

621 super().__init__(dtype=A.dtype, shape=shape) 

622 self.A = A 

623 self.args = (A,) 

624 

625 def _matvec(self, x): 

626 return self.A._rmatvec(x) 

627 

628 def _rmatvec(self, x): 

629 return self.A._matvec(x) 

630 

631 def _matmat(self, x): 

632 return self.A._rmatmat(x) 

633 

634 def _rmatmat(self, x): 

635 return self.A._matmat(x) 

636 

637class _TransposedLinearOperator(LinearOperator): 

638 """Transposition of arbitrary Linear Operator""" 

639 

640 def __init__(self, A): 

641 shape = (A.shape[1], A.shape[0]) 

642 super().__init__(dtype=A.dtype, shape=shape) 

643 self.A = A 

644 self.args = (A,) 

645 

646 def _matvec(self, x): 

647 # NB. np.conj works also on sparse matrices 

648 return np.conj(self.A._rmatvec(np.conj(x))) 

649 

650 def _rmatvec(self, x): 

651 return np.conj(self.A._matvec(np.conj(x))) 

652 

653 def _matmat(self, x): 

654 # NB. np.conj works also on sparse matrices 

655 return np.conj(self.A._rmatmat(np.conj(x))) 

656 

657 def _rmatmat(self, x): 

658 return np.conj(self.A._matmat(np.conj(x))) 

659 

660def _get_dtype(operators, dtypes=None): 

661 if dtypes is None: 

662 dtypes = [] 

663 for obj in operators: 

664 if obj is not None and hasattr(obj, 'dtype'): 

665 dtypes.append(obj.dtype) 

666 return np.result_type(*dtypes) 

667 

668 

669class _SumLinearOperator(LinearOperator): 

670 def __init__(self, A, B): 

671 if not isinstance(A, LinearOperator) or \ 

672 not isinstance(B, LinearOperator): 

673 raise ValueError('both operands have to be a LinearOperator') 

674 if A.shape != B.shape: 

675 raise ValueError(f'cannot add {A} and {B}: shape mismatch') 

676 self.args = (A, B) 

677 super().__init__(_get_dtype([A, B]), A.shape) 

678 

679 def _matvec(self, x): 

680 return self.args[0].matvec(x) + self.args[1].matvec(x) 

681 

682 def _rmatvec(self, x): 

683 return self.args[0].rmatvec(x) + self.args[1].rmatvec(x) 

684 

685 def _rmatmat(self, x): 

686 return self.args[0].rmatmat(x) + self.args[1].rmatmat(x) 

687 

688 def _matmat(self, x): 

689 return self.args[0].matmat(x) + self.args[1].matmat(x) 

690 

691 def _adjoint(self): 

692 A, B = self.args 

693 return A.H + B.H 

694 

695 

696class _ProductLinearOperator(LinearOperator): 

697 def __init__(self, A, B): 

698 if not isinstance(A, LinearOperator) or \ 

699 not isinstance(B, LinearOperator): 

700 raise ValueError('both operands have to be a LinearOperator') 

701 if A.shape[1] != B.shape[0]: 

702 raise ValueError(f'cannot multiply {A} and {B}: shape mismatch') 

703 super().__init__(_get_dtype([A, B]), 

704 (A.shape[0], B.shape[1])) 

705 self.args = (A, B) 

706 

707 def _matvec(self, x): 

708 return self.args[0].matvec(self.args[1].matvec(x)) 

709 

710 def _rmatvec(self, x): 

711 return self.args[1].rmatvec(self.args[0].rmatvec(x)) 

712 

713 def _rmatmat(self, x): 

714 return self.args[1].rmatmat(self.args[0].rmatmat(x)) 

715 

716 def _matmat(self, x): 

717 return self.args[0].matmat(self.args[1].matmat(x)) 

718 

719 def _adjoint(self): 

720 A, B = self.args 

721 return B.H * A.H 

722 

723 

724class _ScaledLinearOperator(LinearOperator): 

725 def __init__(self, A, alpha): 

726 if not isinstance(A, LinearOperator): 

727 raise ValueError('LinearOperator expected as A') 

728 if not np.isscalar(alpha): 

729 raise ValueError('scalar expected as alpha') 

730 if isinstance(A, _ScaledLinearOperator): 

731 A, alpha_original = A.args 

732 # Avoid in-place multiplication so that we don't accidentally mutate 

733 # the original prefactor. 

734 alpha = alpha * alpha_original 

735 

736 dtype = _get_dtype([A], [type(alpha)]) 

737 super().__init__(dtype, A.shape) 

738 self.args = (A, alpha) 

739 

740 def _matvec(self, x): 

741 return self.args[1] * self.args[0].matvec(x) 

742 

743 def _rmatvec(self, x): 

744 return np.conj(self.args[1]) * self.args[0].rmatvec(x) 

745 

746 def _rmatmat(self, x): 

747 return np.conj(self.args[1]) * self.args[0].rmatmat(x) 

748 

749 def _matmat(self, x): 

750 return self.args[1] * self.args[0].matmat(x) 

751 

752 def _adjoint(self): 

753 A, alpha = self.args 

754 return A.H * np.conj(alpha) 

755 

756 

757class _PowerLinearOperator(LinearOperator): 

758 def __init__(self, A, p): 

759 if not isinstance(A, LinearOperator): 

760 raise ValueError('LinearOperator expected as A') 

761 if A.shape[0] != A.shape[1]: 

762 raise ValueError('square LinearOperator expected, got %r' % A) 

763 if not isintlike(p) or p < 0: 

764 raise ValueError('non-negative integer expected as p') 

765 

766 super().__init__(_get_dtype([A]), A.shape) 

767 self.args = (A, p) 

768 

769 def _power(self, fun, x): 

770 res = np.array(x, copy=True) 

771 for i in range(self.args[1]): 

772 res = fun(res) 

773 return res 

774 

775 def _matvec(self, x): 

776 return self._power(self.args[0].matvec, x) 

777 

778 def _rmatvec(self, x): 

779 return self._power(self.args[0].rmatvec, x) 

780 

781 def _rmatmat(self, x): 

782 return self._power(self.args[0].rmatmat, x) 

783 

784 def _matmat(self, x): 

785 return self._power(self.args[0].matmat, x) 

786 

787 def _adjoint(self): 

788 A, p = self.args 

789 return A.H ** p 

790 

791 

792class MatrixLinearOperator(LinearOperator): 

793 def __init__(self, A): 

794 super().__init__(A.dtype, A.shape) 

795 self.A = A 

796 self.__adj = None 

797 self.args = (A,) 

798 

799 def _matmat(self, X): 

800 return self.A.dot(X) 

801 

802 def _adjoint(self): 

803 if self.__adj is None: 

804 self.__adj = _AdjointMatrixOperator(self) 

805 return self.__adj 

806 

807class _AdjointMatrixOperator(MatrixLinearOperator): 

808 def __init__(self, adjoint): 

809 self.A = adjoint.A.T.conj() 

810 self.__adjoint = adjoint 

811 self.args = (adjoint,) 

812 self.shape = adjoint.shape[1], adjoint.shape[0] 

813 

814 @property 

815 def dtype(self): 

816 return self.__adjoint.dtype 

817 

818 def _adjoint(self): 

819 return self.__adjoint 

820 

821 

822class IdentityOperator(LinearOperator): 

823 def __init__(self, shape, dtype=None): 

824 super().__init__(dtype, shape) 

825 

826 def _matvec(self, x): 

827 return x 

828 

829 def _rmatvec(self, x): 

830 return x 

831 

832 def _rmatmat(self, x): 

833 return x 

834 

835 def _matmat(self, x): 

836 return x 

837 

838 def _adjoint(self): 

839 return self 

840 

841 

842def aslinearoperator(A): 

843 """Return A as a LinearOperator. 

844 

845 'A' may be any of the following types: 

846 - ndarray 

847 - matrix 

848 - sparse matrix (e.g. csr_matrix, lil_matrix, etc.) 

849 - LinearOperator 

850 - An object with .shape and .matvec attributes 

851 

852 See the LinearOperator documentation for additional information. 

853 

854 Notes 

855 ----- 

856 If 'A' has no .dtype attribute, the data type is determined by calling 

857 :func:`LinearOperator.matvec()` - set the .dtype attribute to prevent this 

858 call upon the linear operator creation. 

859 

860 Examples 

861 -------- 

862 >>> import numpy as np 

863 >>> from scipy.sparse.linalg import aslinearoperator 

864 >>> M = np.array([[1,2,3],[4,5,6]], dtype=np.int32) 

865 >>> aslinearoperator(M) 

866 <2x3 MatrixLinearOperator with dtype=int32> 

867 """ 

868 if isinstance(A, LinearOperator): 

869 return A 

870 

871 elif isinstance(A, np.ndarray) or isinstance(A, np.matrix): 

872 if A.ndim > 2: 

873 raise ValueError('array must have ndim <= 2') 

874 A = np.atleast_2d(np.asarray(A)) 

875 return MatrixLinearOperator(A) 

876 

877 elif issparse(A) or is_pydata_spmatrix(A): 

878 return MatrixLinearOperator(A) 

879 

880 else: 

881 if hasattr(A, 'shape') and hasattr(A, 'matvec'):