1 """
2
3 Classes
4 -------
5
6 LinearOperator : Provides the LinearOperator class, which is the base
7 class for all operators. This class is shamelessly taken from the
8 scipy.sparse.linalg package.
9
10 Functions
11 ---------
12
13 aslinearoperator : transforms a class with the matvec method into a
14 LinearOperator
15
16 concatenate : Transform a list of LinearOperators into a
17 LinearOperator by concatenation.
18
19 block_diagonal : Transform a list of LinearOperators into a
20 block-diagonal LinearOperator with given LinearOperators on the
21 diagonal.
22 """
23
24 import numpy as np
25 import warnings
32 """Common interface for performing matrix vector products
33
34 Many iterative methods (e.g. cg, gmres) do not need to know the
35 individual entries of a matrix to solve a linear system A*x=b.
36 Such solvers only require the computation of matrix vector
37 products, A*v where v is a dense vector. This class serves as
38 an abstract interface between iterative solvers and matrix-like
39 objects.
40
41 Parameters
42 ----------
43 shape : tuple
44 Matrix dimensions (M,N)
45 matvec : callable f(v)
46 Returns returns A * v.
47
48 Optional Parameters
49 -------------------
50 rmatvec : callable f(v)
51 Returns A^H * v, where A^H is the conjugate transpose of A.
52 matmat : callable f(V)
53 Returns A * V, where V is a dense matrix with dimensions (N,K).
54 rmatmat : callable f(V)
55 Returns A^H * V, where V is a dense matrix with dimensions (N,K).
56 dtype : dtype
57 Data type of the matrix.
58
59 See Also
60 --------
61 aslinearoperator : Construct LinearOperators
62
63 Notes
64 -----
65 The user-defined matvec() function must properly handle the case
66 where v has shape (N,) as well as the (N,1) case. The shape of
67 the return type is handled internally by LinearOperator.
68
69 Examples
70 --------
71 >>> from scipy.sparse.linalg import LinearOperator
72 >>> from scipy import *
73 >>> def mv(v):
74 ... return array([ 2*v[0], 3*v[1]])
75 ...
76 >>> A = LinearOperator( (2,2), matvec=mv, rmatvec=mv)
77 >>> A
78 <2x2 LinearOperator with unspecified dtype>
79 >>> A.matvec( ones(2) )
80 array([ 2., 3.])
81 >>> A * ones(2)
82 array([ 2., 3.])
83 >>> A.dense()
84 array([[ 2., 0.],
85 [ 0., 3.]])
86 >>> (2 * A.T * A + 1) * ones(2)
87 array([ 9., 19.])
88
89 """
90 - def __init__(self, shape, matvec, rmatvec=None, matmat=None, rmatmat=None,
91 dtypein=None, dtypeout=None, dtype=np.float64):
92
93 shape = tuple(shape)
94
95
96
97
98 self.shape = shape
99 self._matvec = matvec
100
101 if rmatvec is None:
102 def rmatvec(v):
103 raise NotImplementedError('rmatvec is not defined')
104 self.rmatvec = rmatvec
105 else:
106 self.rmatvec = rmatvec
107
108 if matmat is not None:
109
110 self._matmat = matmat
111
112 if rmatmat is not None:
113
114 self._rmatmat = rmatmat
115 else:
116 self._rmatmat = None
117
118 self.dtype = None
119 self.dtypein = None
120 self.dtypeout = None
121 if dtype is not None:
122 self.dtype = dtype
123 self.dtypein = np.dtype(dtype)
124 self.dtypeout = np.dtype(dtype)
125 if dtypein is not None and dtypeout is not None:
126 self.dtype = np.dtype(dtypein)
127 self.dtypein = np.dtype(dtypein)
128 self.dtypeout = np.dtype(dtypeout)
129 elif dtypein is not None:
130 self.dtype = np.dtype(dtypein)
131 self.dtypein = np.dtype(dtypein)
132 elif dtypeout is not None:
133 self.dtype = np.dtype(dtypeout)
134 self.dtypeout = np.dtype(dtypeout)
135
137 """Default matrix-matrix multiplication handler. Falls back on
138 the user-defined matvec() routine, which is always provided.
139 """
140
141 return np.hstack( [ self.matvec(col.reshape(-1,1)) for col in X.T ] )
142
143
145 """Matrix-vector multiplication
146
147 Performs the operation y=A*x where A is an MxN linear
148 operator and x is a column vector or rank-1 array.
149
150 Parameters
151 ----------
152 x : {matrix, ndarray}
153 An array with shape (N,) or (N,1).
154
155 Returns
156 -------
157 y : {matrix, ndarray}
158 A matrix or ndarray with shape (M,) or (M,1) depending
159 on the type and shape of the x argument.
160
161 Notes
162 -----
163 This matvec wraps the user-specified matvec routine to ensure that
164 y has the correct shape and type.
165
166 """
167
168 x = np.asanyarray(x)
169
170 M,N = self.shape
171
172 if x.shape != (N,) and x.shape != (N,1):
173 raise ValueError('dimension mismatch')
174
175 y = self._matvec(x)
176
177 if isinstance(x, np.matrix):
178 y = np.asmatrix(y)
179 else:
180 y = np.asarray(y)
181
182 if x.ndim == 1:
183 y = y.reshape(M)
184 elif x.ndim == 2:
185 y = y.reshape(M,1)
186 else:
187 raise ValueError('invalid shape returned by user-defined matvec()')
188
189
190 return y
191
192
194 """Matrix-matrix multiplication
195
196 Performs the operation y=A*X where A is an MxN linear
197 operator and X dense N*K matrix or ndarray.
198
199 Parameters
200 ----------
201 X : {matrix, ndarray}
202 An array with shape (N,K).
203
204 Returns
205 -------
206 Y : {matrix, ndarray}
207 A matrix or ndarray with shape (M,K) depending on
208 the type of the X argument.
209
210 Notes
211 -----
212 This matmat wraps any user-specified matmat routine to ensure that
213 y has the correct type.
214
215 """
216
217 X = np.asanyarray(X)
218
219 if X.ndim != 2:
220 raise ValueError('expected rank-2 ndarray or matrix')
221
222 M,N = self.shape
223
224 if X.shape[0] != N:
225 raise ValueError('dimension mismatch')
226
227 Y = self._matmat(X)
228
229 if isinstance(Y, np.matrix):
230 Y = np.asmatrix(Y)
231
232 return Y
233
234
236 if np.isscalar(x):
237 matvec = _mat_scalar(self._matvec, x)
238 rmatvec = _mat_scalar(self.rmatvec, x)
239 matmat = _mat_scalar(self._matmat, x)
240 rmatmat = _mat_scalar(self._rmatmat, x)
241 return LinearOperator(self.shape, matvec, rmatvec=rmatvec,
242 matmat=matmat, rmatmat=rmatmat)
243 if isinstance(x, LinearOperator):
244
245 if self.shape[1] != x.shape[0]:
246 raise ValueError('LinearOperator shape mismatch')
247 if self.dtypein != x.dtypeout:
248 raise ValueError('LinearOperator dtypein mismatch')
249 shape = (self.shape[0], x.shape[1])
250 matvec = _mat_mul(self._matvec, x._matvec)
251 if self.rmatvec is not None and x.rmatvec is not None:
252 rmatvec = _mat_mul(x.rmatvec, self.rmatvec)
253 else:
254 rmatvec = None
255 if self._matmat is not None and x._matmat is not None:
256 matmat = _mat_mul(self._matmat, x._matmat)
257 else:
258 matmat = None
259 if self._rmatmat is not None and x._matmat is not None:
260 rmatmat = _mat_mul(x._rmatmat, self._rmatmat)
261 else:
262 rmatmat = None
263 return LinearOperator(shape, matvec, rmatvec=rmatvec,
264 matmat=matmat, rmatmat=rmatmat,
265 dtypein=x.dtypein,
266 dtypeout=self.dtypeout)
267 else:
268 x = np.asarray(x)
269 if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1:
270 warnings.filterwarnings("ignore", category=np.ComplexWarning)
271 return self.matvec(x).astype(self.dtypeout)
272 warnings.resetwarnings()
273 elif x.ndim == 2:
274 warnings.filterwarnings("ignore", category=np.ComplexWarning)
275 return self.matmat(x).astype(self.dtypeout)
276 warnings.resetwarnings()
277 else:
278 raise ValueError('expected rank-1 or rank-2 array or matrix or LinearOperator')
279
281 if isinstance(A, LinearOperator):
282 if self.shape != A.shape:
283 raise ValueError('expected LinearOperator of the same shape')
284 if self.dtype != A.dtype:
285 raise ValueError('LinearOperator dtype mismatch')
286 if self.dtypein != A.dtypein:
287 raise ValueError('LinearOperator dtypein mismatch')
288 if self.dtypeout != A.dtypeout:
289 raise ValueError('LinearOperator dtypeout mismatch')
290 matvec = _mat_add(self._matvec, A._matvec)
291 if self.rmatvec is not None and A.rmatvec is not None:
292 rmatvec = _mat_add(self.rmatvec, A.rmatvec)
293 else:
294 rmatvec = None
295 if self._matmat is not None and A._matmat is not None:
296 matmat = _mat_add(self._matmat, A._matmat)
297 else:
298 matmat = None
299 if self._matmat is not None and A._matmat is not None:
300 rmatmat = _mat_add(self._matmat, A._matmat)
301 else:
302 rmatmat = None
303 return LinearOperator(self.shape, matvec, rmatvec=rmatvec,
304 matmat=matmat, rmatmat=rmatmat,
305 dtype=self.dtype, dtypein=self.dtypein,
306 dtypeout=self.dtypeout)
307 if np.isscalar(A):
308 return self.__add__(A * aslinearoperator(np.eye(self.shape[0],
309 self.shape[1],
310 dtype=self.dtype,
311 )))
312 else:
313 raise ValueError('expected LinearOperator')
314
317
320
322 if np.isscalar(x):
323
324 return self.__mul__(x)
325 if isinstance(x, LinearOperator):
326 return x.__mul__(self)
327 else:
328 x = np.asarray(x)
329 if hasattr(self, 'rmatvec'):
330 if x.ndim == 1 or x.ndim == 2 and x.shape[1] == 1:
331 return self.rmatvec(x)
332 elif x.ndim == 2:
333 return self.rmatmat(x)
334 else:
335 raise ValueError('expected rank-1 or rank-2 array or matrix')
336 else:
337 raise ValueError('LinearOperator does not have rmatvec attribute')
338
341
344
347
350
353
355 from copy import copy
356 if not isinstance(k, int):
357 raise NotImplemented('Only power to an int is implemented')
358 if k < 0:
359 raise NotImplemented('Negative power is not implemented')
360 else:
361 A = copy(self)
362 for i in xrange(k):
363 A *= self
364 return A
365
367 M,N = self.shape
368 if hasattr(self,'dtype'):
369 dt = 'dtype=' + str(self.dtypein)
370 else:
371 dt = 'unspecified dtype'
372 class_name = self.__class__.__name__
373 return '<%dx%d %s with %s>' % (M, N, class_name, dt)
374
375 @property
377 if self.rmatvec is not None:
378 matvec, rmatvec = self.rmatvec, self._matvec
379 else:
380 raise NotImplementedError('rmatvec is not defined')
381 if self._matmat is not None and self._rmatmat is not None:
382 matmat, rmatmat = self._rmatmat, self._matmat
383 else:
384 matmat, rmatmat = None, None
385 dtype = getattr(self, 'dtype', None)
386 dtypein = getattr(self, 'dtypein', None)
387 dtypeout = getattr(self, 'dtypeout', None)
388 return LinearOperator(self.shape[::-1], matvec, rmatvec=rmatvec,
389 matmat=matmat, rmatmat=rmatmat, dtype=dtype,
390 dtypein=dtypeout, dtypeout=dtypein)
391
393 A = np.empty(self.shape, dtype=self.dtype)
394 x = np.zeros(self.shape[1], dtype=self.dtype)
395 for i in xrange(A.shape[1]):
396 x[i] = 1
397 A[:, i] = self * x
398 x[:] = 0
399 return A
400
402 """Functional composition of two matvec functions"""
403 def matvec(x):
404 return matvec1(matvec2(x))
405 return matvec
406
408 """Functional addition of two matvec functions"""
409 def matvec(x):
410 return np.squeeze(matvec1(x)) + np.squeeze(matvec2(x))
411 return matvec
412
414 def matvec(x):
415 return scalar * matvec0(x)
416 return matvec
417
419 """Return A as a LinearOperator.
420
421 'A' may be any of the following types:
422 - ndarray
423 - matrix
424 - sparse matrix (e.g. csr_matrix, lil_matrix, etc.)
425 - LinearOperator
426 - An object with .shape and .matvec attributes
427
428 See the LinearOperator documentation for additonal information.
429
430 Examples
431 --------
432 >>> from scipy import matrix
433 >>> M = matrix( [[1,2,3],[4,5,6]], dtype='int32' )
434 >>> aslinearoperator( M )
435 <2x3 LinearOperator with dtype=int32>
436
437 """
438 if isinstance(A, LinearOperator):
439 return A
440
441 elif isinstance(A, np.ndarray) or isinstance(A, np.matrix):
442 if A.ndim > 2:
443 raise ValueError('array must have rank <= 2')
444
445 A = np.atleast_2d(np.asarray(A))
446
447 def matvec(v):
448 return np.dot(A, v)
449 def rmatvec(v):
450 return np.dot(A.conj().transpose(), v)
451 def matmat(V):
452 return np.dot(A, V)
453 def rmatmat(V):
454 return np.dot(V, A)
455 return LinearOperator(A.shape, matvec, rmatvec=rmatvec,
456 matmat=matmat, rmatmat=rmatmat, dtype=A.dtype)
457
458
459
460
461
462
463
464
465
466
467
468
469
470 else:
471 if hasattr(A, 'shape') and hasattr(A, 'matvec'):
472 rmatvec = getattr(A, 'rmatvec', None)
473 matmat = getattr(A, 'matmat', None)
474 rmatmat = getattr(A, 'rmatmat', None)
475 dtype = getattr(A, 'dtype', None)
476 return LinearOperator(A.shape, A.matvec, rmatvec=rmatvec,
477 matmat=matmat, rmatmat=rmatmat, dtype=dtype)
478 else:
479 raise TypeError('type not understood')
480
482 """
483 Concatenate LinearOperator in rows or in columns.
484
485 Parameters
486 ----------
487 As : list of LinearOperators
488 The list of LinearOperators to concatenate.
489
490 axis : 0 or 1
491 The axis along which to concatenate the LinearOperators.
492
493 Returns
494 -------
495 out: LinearOperator
496 Output a LinearOperator which is the concatenation of a list of
497 LinearOpeartors.
498
499 """
500
501 if axis == 0:
502
503 for A in As:
504
505 if A.shape[1] != As[0].shape[1]:
506 raise ValueError("All LinearOperators must have the same number of row/columns.")
507
508 if A.dtype != As[0].dtype:
509 raise ValueError("All LinearOperators must have the same data-type.")
510
511 sizes = [A.shape[0] for A in As]
512 shape = np.sum(sizes), As[0].shape[1]
513
514 dtype = As[0].dtype
515
516 def matvec(x):
517 return np.concatenate([A.matvec(x) for A in As])
518
519 sizesum = list(np.cumsum(sizes))[:-1]
520 sizes1 = [None,] + sizesum
521 sizes2 = sizesum + [None,]
522 slices = [slice(s1, s2, None) for s1, s2 in zip(sizes1, sizes2)]
523 def rmatvec(x):
524 out = np.zeros(shape[1])
525 for A, s in zip(As, slices):
526 out += A.rmatvec(x[s])
527 return out
528
529 elif axis == 1:
530
531 return concatenate([A.T for A in As], axis=0).T
532 else:
533 raise ValueError("axis must be 0 or 1")
534 return LinearOperator(shape, matvec, rmatvec=rmatvec, dtype=dtype)
535
537 """
538 Defines a block diagonal LinearOperator from a list of LinearOperators.
539 """
540
541 for A in As:
542 if A.dtype != As[0].dtype:
543 raise ValueError("All LinearOperators must have the same data-type.")
544 dtype = As[0].dtype
545
546 shape = np.sum([A.shape[0] for A in As]), np.sum([A.shape[1] for A in As])
547
548 sizes = [A.shape[1] for A in As]
549 sizes = list(np.cumsum(sizes))[:-1]
550 sizes1 = [None,] + sizes
551 sizes2 = sizes + [None,]
552 slices = [slice(s1, s2, None) for s1, s2 in zip(sizes1, sizes2)]
553
554 def matvec(x):
555 return np.concatenate([A.matvec(x[s]) for A, s in zip(As, slices)])
556
557 rsizes = [A.shape[0] for A in As]
558 rsizes = list(np.cumsum(rsizes))[:-1]
559 rsizes1 = [None,] + rsizes
560 rsizes2 = rsizes + [None,]
561 rslices = [slice(s1, s2, None) for s1, s2 in zip(rsizes1, rsizes2)]
562
563 def rmatvec(x):
564 return np.concatenate([A.rmatvec(x[s]) for A, s in zip(As, rslices)])
565 return LinearOperator(shape, matvec, rmatvec=rmatvec, dtype=dtype)
566