Coverage for /usr/lib/python3/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 >>> class Ones(LinearOperator): 

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

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

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

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

21 

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

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

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

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

26 

27 >>> from scipy.sparse import csr_matrix 

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

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

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

31 array([13, 4, 15]) 

32 

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

34counterpart:: 

35 

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

37 array([13, 4, 15]) 

38 

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

40``LinearOperator`` instances. 

41""" 

42 

43import warnings 

44 

45import numpy as np 

46 

47from scipy.sparse import issparse 

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

49 

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

51 

52 

53class LinearOperator: 

54 """Common interface for performing matrix vector products 

55 

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

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

58 Such solvers only require the computation of matrix vector 

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

60 an abstract interface between iterative solvers and matrix-like 

61 objects. 

62 

63 To construct a concrete LinearOperator, either pass appropriate 

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

65 

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

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

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

69 on this class to have these attributes validated. Implementing 

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

71 algorithm) and vice-versa. 

72 

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

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

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

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

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

78 backwards compatibility. 

79 

80 Parameters 

81 ---------- 

82 shape : tuple 

83 Matrix dimensions (M, N). 

84 matvec : callable f(v) 

85 Returns returns A * v. 

86 rmatvec : callable f(v) 

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

88 matmat : callable f(V) 

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

90 dtype : dtype 

91 Data type of the matrix. 

92 rmatmat : callable f(V) 

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

94 

95 Attributes 

96 ---------- 

97 args : tuple 

98 For linear operators describing products etc. of other linear 

99 operators, the operands of the binary operation. 

100 ndim : int 

101 Number of dimensions (this is always 2) 

102 

103 See Also 

104 -------- 

105 aslinearoperator : Construct LinearOperators 

106 

107 Notes 

108 ----- 

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

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

111 the return type is handled internally by LinearOperator. 

112 

113 LinearOperator instances can also be multiplied, added with each 

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

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

116 operations to the original operators and combines the results. 

117 

118 More details regarding how to subclass a LinearOperator and several 

119 examples of concrete LinearOperator instances can be found in the 

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

121 

122 

123 Examples 

124 -------- 

125 >>> import numpy as np 

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

127 >>> def mv(v): 

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

129 ... 

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

131 >>> A 

132 <2x2 _CustomLinearOperator with dtype=float64> 

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

134 array([ 2., 3.]) 

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

136 array([ 2., 3.]) 

137 

138 """ 

139 

140 ndim = 2 

141 # Necessary for right matmul with numpy arrays. 

142 __array_ufunc__ = None 

143 

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

145 if cls is LinearOperator: 

146 # Operate as _CustomLinearOperator factory. 

147 return super().__new__(_CustomLinearOperator) 

148 else: 

149 obj = super().__new__(cls) 

150 

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

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

153 warnings.warn("LinearOperator subclass should implement" 

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

155 category=RuntimeWarning, stacklevel=2) 

156 

157 return obj 

158 

159 def __init__(self, dtype, shape): 

160 """Initialize this LinearOperator. 

161 

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

163 be convertible to a length-2 tuple. 

164 """ 

165 if dtype is not None: 

166 dtype = np.dtype(dtype) 

167 

168 shape = tuple(shape) 

169 if not isshape(shape): 

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

171 

172 self.dtype = dtype 

173 self.shape = shape 

174 

175 def _init_dtype(self): 

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

177 """ 

178 if self.dtype is None: 

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

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

181 

182 def _matmat(self, X): 

183 """Default matrix-matrix multiplication handler. 

184 

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

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

187 """ 

188 

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

190 

191 def _matvec(self, x): 

192 """Default matrix-vector multiplication handler. 

193 

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

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

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

197 

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

199 will define matrix-vector multiplication as well. 

200 """ 

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

202 

203 def matvec(self, x): 

204 """Matrix-vector multiplication. 

205 

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

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

208 

209 Parameters 

210 ---------- 

211 x : {matrix, ndarray} 

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

213 

214 Returns 

215 ------- 

216 y : {matrix, ndarray} 

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

218 on the type and shape of the x argument. 

219 

220 Notes 

221 ----- 

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

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

224 

225 """ 

226 

227 x = np.asanyarray(x) 

228 

229 M,N = self.shape 

230 

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

232 raise ValueError('dimension mismatch') 

233 

234 y = self._matvec(x) 

235 

236 if isinstance(x, np.matrix): 

237 y = asmatrix(y) 

238 else: 

239 y = np.asarray(y) 

240 

241 if x.ndim == 1: 

242 y = y.reshape(M) 

243 elif x.ndim == 2: 

244 y = y.reshape(M,1) 

245 else: 

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

247 

248 return y 

249 

250 def rmatvec(self, x): 

251 """Adjoint matrix-vector multiplication. 

252 

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

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

255 

256 Parameters 

257 ---------- 

258 x : {matrix, ndarray} 

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

260 

261 Returns 

262 ------- 

263 y : {matrix, ndarray} 

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

265 on the type and shape of the x argument. 

266 

267 Notes 

268 ----- 

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

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

271 

272 """ 

273 

274 x = np.asanyarray(x) 

275 

276 M,N = self.shape 

277 

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

279 raise ValueError('dimension mismatch') 

280 

281 y = self._rmatvec(x) 

282 

283 if isinstance(x, np.matrix): 

284 y = asmatrix(y) 

285 else: 

286 y = np.asarray(y) 

287 

288 if x.ndim == 1: 

289 y = y.reshape(N) 

290 elif x.ndim == 2: 

291 y = y.reshape(N,1) 

292 else: 

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

294 

295 return y 

296 

297 def _rmatvec(self, x): 

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

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

300 # _adjoint not overridden, prevent infinite recursion 

301 raise NotImplementedError 

302 else: 

303 return self.H.matvec(x) 

304 

305 def matmat(self, X): 

306 """Matrix-matrix multiplication. 

307 

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

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

310 

311 Parameters 

312 ---------- 

313 X : {matrix, ndarray} 

314 An array with shape (N,K). 

315 

316 Returns 

317 ------- 

318 Y : {matrix, ndarray} 

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

320 the type of the X argument. 

321 

322 Notes 

323 ----- 

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

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

326 

327 """ 

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

329 X = np.asanyarray(X) 

330 

331 if X.ndim != 2: 

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

333 

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

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

336 

337 try: 

338 Y = self._matmat(X) 

339 except Exception as e: 

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

341 raise TypeError( 

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

343 " Wrap the matrix in aslinearoperator first." 

344 ) from e 

345 raise 

346 

347 if isinstance(Y, np.matrix): 

348 Y = asmatrix(Y) 

349 

350 return Y 

351 

352 def rmatmat(self, X): 

353 """Adjoint matrix-matrix multiplication. 

354 

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

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

357 The default implementation defers to the adjoint. 

358 

359 Parameters 

360 ---------- 

361 X : {matrix, ndarray} 

362 A matrix or 2D array. 

363 

364 Returns 

365 ------- 

366 Y : {matrix, ndarray} 

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

368 

369 Notes 

370 ----- 

371 This rmatmat wraps the user-specified rmatmat routine. 

372 

373 """ 

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

375 X = np.asanyarray(X) 

376 

377 if X.ndim != 2: 

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

379 % X.ndim) 

380 

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

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

383 

384 try: 

385 Y = self._rmatmat(X) 

386 except Exception as e: 

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

388 raise TypeError( 

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

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

391 ) from e 

392 raise 

393 

394 if isinstance(Y, np.matrix): 

395 Y = asmatrix(Y) 

396 return Y 

397 

398 def _rmatmat(self, X): 

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

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

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

402 else: 

403 return self.H.matmat(X) 

404 

405 def __call__(self, x): 

406 return self*x 

407 

408 def __mul__(self, x): 

409 return self.dot(x) 

410 

411 def __truediv__(self, other): 

412 if not np.isscalar(other): 

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

414 

415 return _ScaledLinearOperator(self, 1.0/other) 

416 

417 def dot(self, x): 

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

419 

420 Parameters 

421 ---------- 

422 x : array_like 

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

424 

425 Returns 

426 ------- 

427 Ax : array 

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

429 the result of applying this linear operator on x. 

430 

431 """ 

432 if isinstance(x, LinearOperator): 

433 return _ProductLinearOperator(self, x) 

434 elif np.isscalar(x): 

435 return _ScaledLinearOperator(self, x) 

436 else: 

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

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

439 x = np.asarray(x) 

440 

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

442 return self.matvec(x) 

443 elif x.ndim == 2: 

444 return self.matmat(x) 

445 else: 

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

447 % x) 

448 

449 def __matmul__(self, other): 

450 if np.isscalar(other): 

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

452 "use '*' instead") 

453 return self.__mul__(other) 

454 

455 def __rmatmul__(self, other): 

456 if np.isscalar(other): 

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

458 "use '*' instead") 

459 return self.__rmul__(other) 

460 

461 def __rmul__(self, x): 

462 if np.isscalar(x): 

463 return _ScaledLinearOperator(self, x) 

464 else: 

465 return self._rdot(x) 

466 

467 def _rdot(self, x): 

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

469 

470 Parameters 

471 ---------- 

472 x : array_like 

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

474 

475 Returns 

476 ------- 

477 xA : array 

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

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

480 

481 Notes 

482 ----- 

483 This is copied from dot to implement right multiplication. 

484 """ 

485 if isinstance(x, LinearOperator): 

486 return _ProductLinearOperator(x, self) 

487 elif np.isscalar(x): 

488 return _ScaledLinearOperator(self, x) 

489 else: 

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

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

492 x = np.asarray(x) 

493 

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

495 # unnecessary complex conjugation if possible. 

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

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

498 elif x.ndim == 2: 

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

500 else: 

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

502 % x) 

503 

504 def __pow__(self, p): 

505 if np.isscalar(p): 

506 return _PowerLinearOperator(self, p) 

507 else: 

508 return NotImplemented 

509 

510 def __add__(self, x): 

511 if isinstance(x, LinearOperator): 

512 return _SumLinearOperator(self, x) 

513 else: 

514 return NotImplemented 

515 

516 def __neg__(self): 

517 return _ScaledLinearOperator(self, -1) 

518 

519 def __sub__(self, x): 

520 return self.__add__(-x) 

521 

522 def __repr__(self): 

523 M,N = self.shape 

524 if self.dtype is None: 

525 dt = 'unspecified dtype' 

526 else: 

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

528 

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

530 

531 def adjoint(self): 

532 """Hermitian adjoint. 

533 

534 Returns the Hermitian adjoint of self, aka the Hermitian 

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

536 Hermitian adjoint is equal to the conjugate transpose. 

537 

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

539 

540 Returns 

541 ------- 

542 A_H : LinearOperator 

543 Hermitian adjoint of self. 

544 """ 

545 return self._adjoint() 

546 

547 H = property(adjoint) 

548 

549 def transpose(self): 

550 """Transpose this linear operator. 

551 

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

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

554 """ 

555 return self._transpose() 

556 

557 T = property(transpose) 

558 

559 def _adjoint(self): 

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

561 return _AdjointLinearOperator(self) 

562 

563 def _transpose(self): 

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

565 return _TransposedLinearOperator(self) 

566 

567 

568class _CustomLinearOperator(LinearOperator): 

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

570 

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

572 dtype=None, rmatmat=None): 

573 super().__init__(dtype, shape) 

574 

575 self.args = () 

576 

577 self.__matvec_impl = matvec 

578 self.__rmatvec_impl = rmatvec 

579 self.__rmatmat_impl = rmatmat 

580 self.__matmat_impl = matmat 

581 

582 self._init_dtype() 

583 

584 def _matmat(self, X): 

585 if self.__matmat_impl is not None: 

586 return self.__matmat_impl(X) 

587 else: 

588 return super()._matmat(X) 

589 

590 def _matvec(self, x): 

591 return self.__matvec_impl(x) 

592 

593 def _rmatvec(self, x): 

594 func = self.__rmatvec_impl 

595 if func is None: 

596 raise NotImplementedError("rmatvec is not defined") 

597 return self.__rmatvec_impl(x) 

598 

599 def _rmatmat(self, X): 

600 if self.__rmatmat_impl is not None: 

601 return self.__rmatmat_impl(X) 

602 else: 

603 return super()._rmatmat(X) 

604 

605 def _adjoint(self): 

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

607 matvec=self.__rmatvec_impl, 

608 rmatvec=self.__matvec_impl, 

609 matmat=self.__rmatmat_impl, 

610 rmatmat=self.__matmat_impl, 

611 dtype=self.dtype) 

612 

613 

614class _AdjointLinearOperator(LinearOperator): 

615 """Adjoint of arbitrary Linear Operator""" 

616 

617 def __init__(self, A): 

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

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

620 self.A = A 

621 self.args = (A,) 

622 

623 def _matvec(self, x): 

624 return self.A._rmatvec(x) 

625 

626 def _rmatvec(self, x): 

627 return self.A._matvec(x) 

628 

629 def _matmat(self, x): 

630 return self.A._rmatmat(x) 

631 

632 def _rmatmat(self, x): 

633 return self.A._matmat(x) 

634 

635class _TransposedLinearOperator(LinearOperator): 

636 """Transposition of arbitrary Linear Operator""" 

637 

638 def __init__(self, A): 

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

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

641 self.A = A 

642 self.args = (A,) 

643 

644 def _matvec(self, x): 

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

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

647 

648 def _rmatvec(self, x): 

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

650 

651 def _matmat(self, x): 

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

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

654 

655 def _rmatmat(self, x): 

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

657 

658def _get_dtype(operators, dtypes=None): 

659 if dtypes is None: 

660 dtypes = [] 

661 for obj in operators: 

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

663 dtypes.append(obj.dtype) 

664 return np.result_type(*dtypes) 

665 

666 

667class _SumLinearOperator(LinearOperator): 

668 def __init__(self, A, B): 

669 if not isinstance(A, LinearOperator) or \ 

670 not isinstance(B, LinearOperator): 

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

672 if A.shape != B.shape: 

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

674 self.args = (A, B) 

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

676 

677 def _matvec(self, x): 

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

679 

680 def _rmatvec(self, x): 

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

682 

683 def _rmatmat(self, x): 

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

685 

686 def _matmat(self, x): 

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

688 

689 def _adjoint(self): 

690 A, B = self.args 

691 return A.H + B.H 

692 

693 

694class _ProductLinearOperator(LinearOperator): 

695 def __init__(self, A, B): 

696 if not isinstance(A, LinearOperator) or \ 

697 not isinstance(B, LinearOperator): 

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

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

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

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

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

703 self.args = (A, B) 

704 

705 def _matvec(self, x): 

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

707 

708 def _rmatvec(self, x): 

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

710 

711 def _rmatmat(self, x): 

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

713 

714 def _matmat(self, x): 

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

716 

717 def _adjoint(self): 

718 A, B = self.args 

719 return B.H * A.H 

720 

721 

722class _ScaledLinearOperator(LinearOperator): 

723 def __init__(self, A, alpha): 

724 if not isinstance(A, LinearOperator): 

725 raise ValueError('LinearOperator expected as A') 

726 if not np.isscalar(alpha): 

727 raise ValueError('scalar expected as alpha') 

728 if isinstance(A, _ScaledLinearOperator): 

729 A, alpha_original = A.args 

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

731 # the original prefactor. 

732 alpha = alpha * alpha_original 

733 

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

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

736 self.args = (A, alpha) 

737 

738 def _matvec(self, x): 

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

740 

741 def _rmatvec(self, x): 

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

743 

744 def _rmatmat(self, x): 

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

746 

747 def _matmat(self, x): 

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

749 

750 def _adjoint(self): 

751 A, alpha = self.args 

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

753 

754 

755class _PowerLinearOperator(LinearOperator): 

756 def __init__(self, A, p): 

757 if not isinstance(A, LinearOperator): 

758 raise ValueError('LinearOperator expected as A') 

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

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

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

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

763 

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

765 self.args = (A, p) 

766 

767 def _power(self, fun, x): 

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

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

770 res = fun(res) 

771 return res 

772 

773 def _matvec(self, x): 

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

775 

776 def _rmatvec(self, x): 

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

778 

779 def _rmatmat(self, x): 

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

781 

782 def _matmat(self, x): 

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

784 

785 def _adjoint(self): 

786 A, p = self.args 

787 return A.H ** p 

788 

789 

790class MatrixLinearOperator(LinearOperator): 

791 def __init__(self, A): 

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

793 self.A = A 

794 self.__adj = None 

795 self.args = (A,) 

796 

797 def _matmat(self, X): 

798 return self.A.dot(X) 

799 

800 def _adjoint(self): 

801 if self.__adj is None: 

802 self.__adj = _AdjointMatrixOperator(self) 

803 return self.__adj 

804 

805class _AdjointMatrixOperator(MatrixLinearOperator): 

806 def __init__(self, adjoint): 

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

808 self.__adjoint = adjoint 

809 self.args = (adjoint,) 

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

811 

812 @property 

813 def dtype(self): 

814 return self.__adjoint.dtype 

815 

816 def _adjoint(self): 

817 return self.__adjoint 

818 

819 

820class IdentityOperator(LinearOperator): 

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

822 super().__init__(dtype, shape) 

823 

824 def _matvec(self, x): 

825 return x 

826 

827 def _rmatvec(self, x): 

828 return x 

829 

830 def _rmatmat(self, x): 

831 return x 

832 

833 def _matmat(self, x): 

834 return x 

835 

836 def _adjoint(self): 

837 return self 

838 

839 

840def aslinearoperator(A): 

841 """Return A as a LinearOperator. 

842 

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

844 - ndarray 

845 - matrix 

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

847 - LinearOperator 

848 - An object with .shape and .matvec attributes 

849 

850 See the LinearOperator documentation for additional information. 

851 

852 Notes 

853 ----- 

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

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

856 call upon the linear operator creation. 

857 

858 Examples 

859 -------- 

860 >>> import numpy as np 

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

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

863 >>> aslinearoperator(M) 

864 <2x3 MatrixLinearOperator with dtype=int32> 

865 """ 

866 if isinstance(A, LinearOperator): 

867 return A 

868 

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

870 if A.ndim > 2: 

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

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

873 return MatrixLinearOperator(A) 

874 

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