linear_operators.operators

1 """Definition of specialized LinearOperator subclasses. 2 3 Classes 4 ------- 5 6 - SymmetricOperator : Enforces matvec = rmatvec and square shape. 7 8 - IdentityOperator : Identity operator. 9 10 - HomotheticOperator : Homothetic operators : i.e. of the form a * I. 11 12 - DiagonalOperator : Diagonal operator. 13 14 - MaskOperator : Diagonal operator with ones and zeros on the diagonal. 15 16 - ShiftOperator : Shift the input vector values. 17 18 - PermutationOperator : Perform permutations of the input vector values. 19 20 - FftOperator : Perform one-dimensional Fast-Fourier Transform (np.fft.fft) 21 22 - ReplicationOperator : Replicates the input vector n times. 23 24 - SliceOperator : Performs slicing on an array. 25 26 - TridiagonalOperator : A tridiagonal matrix a 3 1d arrays. 27 28 - SymmetricTridiagonal : A symmetric tridiagonal matrix a 2 1d arrays. 29 30 - BandOperator : A band matrix using lapack storage scheme. 31 32 - SymmetricBandOperator : A symmetric band matrix using lapack storage scheme. 33 34 - LoweTriangularOperator : A special case of band operator with only subdiagonal bands. 35 36 - UpperTriangularOperator : A special case of band operator with only superdiagonal bands. 37 38 - EigendecompositionOperator : A class storing eigenpairs of another operator. 39 40 Functions 41 --------- 42 43 Functions generate instances of the associated classes. The following are available : 44 45 - identity 46 - homothetic 47 - diagonal 48 - mask 49 - shift 50 - permutation 51 - fft 52 - replication 53 - slice_operator 54 - tridiagonal 55 56 - band_approximation : Generate a BandOperator from a dense matrix. 57 - symmetric_band_approximation : Generate a BandOperator from a dense matrix. 58 59 """ 60 import numpy as np 61 from copy import copy 62 from interface import LinearOperator, concatenate

63 64 # subclasses operators 65 66 -class SymmetricOperator(LinearOperator):

67 """ 68 Subclass of LinearOperator for the definition of symmetric 69 operators, i.e. operators such that : M^T == M. It only requires 70 a shape and a matvec function, since the rmatvec should be the 71 same function as matvec. 72 """

73 - def __init__(self, shape, matvec, **kwargs):

74 """Returns a SymmetricOperator of given shape and matvec 75 function. 76 77 Parameters 78 ---------- 79 80 shape : length 2 tuple 81 The shape of the operator. Should be square. 82 83 matvec : function 84 The matrix-vector operation. 85 86 Returns 87 ------- 88 A SymmetricOperator instance. 89 """ 90 if shape[0] != shape[1]: 91 raise ValueError("Symmetric operators are square operators.") 92 kwargs['rmatvec'] = matvec 93 LinearOperator.__init__(self, shape, matvec, **kwargs)

94

95 -class IdentityOperator(SymmetricOperator):

96 """SymmetricOperator with identity matvec (lambda x: x). 97 """

98 - def __init__(self, shape, **kwargs):

99 """Instantiate an IdentityOperator. 100 101 Parameters 102 ---------- 103 104 shape : length 2 tuple 105 The shape of the operator. Should be square. 106 107 Returns 108 ------- 109 An IdentityOperator instance. 110 111 Exemple 112 ------- 113 >>> import linear_operators as lo 114 >>> I = lo.IdentityOperator((3, 3)) 115 >>> I.todense() 116 array([[ 1., 0., 0.], 117 [ 0., 1., 0.], 118 [ 0., 0., 1.]]) 119 120 """ 121 matvec = lambda x: x 122 SymmetricOperator.__init__(self, shape, matvec, **kwargs)

123

124 -class HomotheticOperator(SymmetricOperator):

125 """ 126 Generate a SymmetricOperator performing an homothety, i.e. a 127 multiplication of the input vector by a scalar. 128 129 Arguments 130 ---------- 131 coef : int float 132 Multiplication coefficient 133 """

134 - def __init__(self, shape, coef, **kwargs):

135 """ 136 Parameters 137 ---------- 138 shape : length 2 tuple 139 The shape of the operator. Should be square. 140 141 coef : int or float 142 The multiplication coefficient. 143 144 Returns 145 ------- 146 An HomotheticOperator instance. 147 148 Exemple 149 ------- 150 >>> import linear_operators as lo 151 >>> H = lo.HomotheticOperator((3, 3), 2.) 152 >>> H.todense() 153 array([[ 2., 0., 0.], 154 [ 0., 2., 0.], 155 [ 0., 0., 2.]]) 156 157 """ 158 self.coef = coef 159 matvec = lambda x: coef * x 160 SymmetricOperator.__init__(self, shape, matvec, **kwargs)

161 - def __repr__(self):

162 s = LinearOperator.__repr__(self) 163 return s[:-1] + " and coef=%f >" % self.coef

164

165 -class DiagonalOperator(SymmetricOperator):

166 """ 167 An operator which mimics a diagonal matrix. 168 169 Attributes 170 ---------- 171 diag: ndarray of 1 dim 172 The diagonal of the matrix. 173 """

174 - def __init__(self, diag, **kwargs):

175 """ 176 Parameters 177 ---------- 178 diag : ndarray with ndim == 1. 179 The diagonal of the operator. 180 181 Returns 182 ------- 183 A DiagonalOperator instance. 184 185 Exemple 186 ------- 187 >>> import numpy as np 188 >>> import linear_operators as lo 189 >>> D = lo.DiagonalOperator(np.arange(3)) 190 >>> D.todense() 191 array([[ 0., 0., 0.], 192 [ 0., 1., 0.], 193 [ 0., 0., 2.]]) 194 """ 195 shape = 2 * (diag.size,) 196 self.diag = diag 197 matvec = lambda x: x * diag 198 SymmetricOperator.__init__(self, shape, matvec, **kwargs)

199 - def __repr__(self):

200 s = LinearOperator.__repr__(self) 201 return s[:-1] + ",\n diagonal=" + self.diag.__repr__() + ">"

202

203 -class MaskOperator(DiagonalOperator):

204 """ 205 A subclass of DiagonalOperator with a boolean diagonal. Elements 206 corresponding to zeros of the diagonal are masked out in the 207 output vector. 208 209 Attributes 210 ----------- 211 mask : ndarray of type bool and ndim==1 212 """

213 - def __init__(self, mask, **kwargs):

214 """ 215 Parameters 216 ---------- 217 mask : ndarray of ones and zeros. 218 If mask[i] = 0, the corresponding value will be masked. 219 If mask is not a boolean array, it is converted to boolean. 220 221 Returns 222 ------- 223 A MaskOperator instance. 224 225 Exemple 226 ------- 227 >>> import numpy as np 228 >>> import linear_operators as lo 229 >>> M = lo.MaskOperator(np.arange(4) % 2) 230 >>> M.todense() 231 array([[ 0., 0., 0., 0.], 232 [ 0., 1., 0., 0.], 233 [ 0., 0., 0., 0.], 234 [ 0., 0., 0., 1.]]) 235 236 """ 237 self.mask = mask.astype(np.bool) 238 DiagonalOperator.__init__(self, self.mask, **kwargs)

239

240 -class ShiftOperator(LinearOperator):

241 """ 242 A LinearOperator corresponding to a shift matrix 243 244 Attributes 245 ---------- 246 shift: signed int 247 Output vector elements are shifted according to the shift attribute. 248 If positive, shift toward superdiagonals, and toward subdiagonal if negative. 249 250 References 251 ---------- 252 http://en.wikipedia.org/wiki/Shift_matrix 253 """

254 - def __init__(self, shape, shift, **kwargs):

255 if shape[0] != shape[1]: 256 raise ValueError("Only square operator is implemented.") 257 if np.round(shift) != shift: 258 raise ValueError("The shift argument should be an integer value.") 259 self.shift = shift 260 ashift = np.abs(shift) 261 # square case 262 matvec = lambda x: np.concatenate((x[ashift:], np.zeros(ashift))) 263 rmatvec = lambda x: np.concatenate((np.zeros(ashift), x[:-ashift])) 264 if self.shift < 0: 265 matvec, rmatvec = rmatvec, matvec 266 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs)

267

268 - def __repr__(self):

269 s = LinearOperator.__repr__(self) 270 return s[:-1] + " and shift=%d >" % self.shift

271

272 -class PermutationOperator(LinearOperator):

273 """ 274 Perform permutations to the vector elements. 275 276 Attributes 277 ---------- 278 perm : list or tuple 279 The permutation of coefficients. 280 """

281 - def __init__(self, perm, **kwargs):

282 """ 283 Parameters 284 ---------- 285 perm : list or tuple 286 The permutation of coefficients. 287 288 Returns 289 ------- 290 A PermutationOperator instance. 291 292 Exemple 293 ------- 294 >>> import numpy as np 295 >>> import linear_operators as lo 296 >>> P = lo.PermutationOperator([3, 1, 0, 2]) 297 >>> P.todense() 298 array([[ 0., 0., 0., 1.], 299 [ 0., 1., 0., 0.], 300 [ 1., 0., 0., 0.], 301 [ 0., 0., 1., 0.]]) 302 303 """ 304 305 shape = 2 * (len(perm),) 306 self.perm = perm 307 # reverse permutation 308 self.perm_inv = np.argsort(perm) 309 def matvec(x): 310 return np.asarray([x[pi] for pi in self.perm])

311 def rmatvec(x): 312 return np.asarray([x[pi] for pi in self.perm_inv])

313 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs)

314 - def __repr__(self):

315 s = LinearOperator.__repr__(self) 316 return s[:-1] + ",\n permutation=" + self.perm.__repr__() + ">"

317

318 -class FftOperator(LinearOperator):

319 """ 320 Generate an Operator performing Fast Fourier Transform. 321 Uses rescaled numpy.fft.ftt and numpy.fft.ifft as matvec and rmatvec. 322 """

323 - def __init__(self, shape, **kwargs):

324 matvec = lambda x: np.fft.fft(x, n=shape[0]) / np.sqrt(shape[0]) 325 rmatvec = lambda x: np.fft.ifft(x, n=shape[1]) * np.sqrt(shape[0]) 326 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs)

327

328 -class ReplicationOperator(LinearOperator):

329 """ 330 Generate an operator which replicates the input vector n times. 331 """

332 - def __init__(self, shape, rep_num, **kwargs):

333 """ 334 Parameters 335 ---------- 336 shape : length 2 tuple. 337 The shape of the operator. 338 339 rep_num : int 340 The number of replications. 341 342 Returns 343 ------- 344 A ReplicationOperator instance. 345 346 Exemple 347 ------- 348 >>> import numpy as np 349 >>> import linear_operators as lo 350 >>> R = lo.ReplicationOperator((4, 2), 2) 351 >>> R.todense() 352 array([[ 1., 0.], 353 [ 0., 1.], 354 [ 1., 0.], 355 [ 0., 1.]]) 356 """ 357 self.rep_num = rep_num 358 # ensure coherent input 359 if not shape[0] == shape[1] * rep_num: 360 raise ValueError("Output vector should be n times the size of input vector.") 361 def matvec(x): 362 return np.tile(x, self.rep_num)

363 def rmatvec(x): 364 N = shape[1] 365 return sum([x[i * N:(i + 1) * N] for i in xrange(self.rep_num)])

366 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs)

367 - def __repr__(self):

368 s = LinearOperator.__repr__(self) 369 return s[:-1] + ",\n Replicate %i times" % self.n + ">"

370

371 -class SliceOperator(LinearOperator):

372 """ 373 Perform slicing on the input vector. 374 375 Attributes 376 ---------- 377 slice: slice 378 Slice applied to the input vector. 379 """

380 - def __init__(self, shape, slice, **kwargs):

381 """ 382 383 Exemple 384 ------- 385 >>> S = lo.SliceOperator((2, 4), slice(None, None, 2)) 386 >>> S.todense() 387 array([[ 1., 0., 0., 0.], 388 [ 0., 0., 1., 0.]]) 389 390 >>> S.T.todense() 391 array([[ 1., 0.], 392 [ 0., 0.], 393 [ 0., 1.], 394 [ 0., 0.]]) 395 396 """ 397 def matvec(x): 398 return x[slice]

399 def rmatvec(x): 400 out = np.zeros(shape[1]) 401 out[slice] = x 402 return out

403 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs) 404

405 -class TridiagonalOperator(LinearOperator):

406 """ 407 Store a tridiagonal operator in the form of 3 arrays 408 409 Attributes 410 ----------- 411 shape : length 2 tuple. 412 The shape of the operator. 413 414 diag : ndarray of size shape[0] 415 The diagonal of the matrix. 416 417 subdiag : ndarray of size shape[0] - 1 418 The subdiagonal of the matrix. 419 420 superdiag : ndarray of size shape[0] - 1 421 The superdiagonal of the matrix. 422 """

423 - def __init__(self, shape, diag, subdiag, superdiag, **kwargs):

424 """ 425 Parameters 426 ---------- 427 shape : length 2 tuple. 428 The shape of the operator. 429 430 diag : ndarray of size shape[0] 431 The diagonal of the matrix. 432 433 subdiag : ndarray of size shape[0] - 1 434 The subdiagonal of the matrix. 435 436 superdiag : ndarray of size shape[0] - 1 437 The superdiagonal of the matrix. 438 439 Returns 440 ------- 441 A Tridiagonal matrix operator instance. 442 443 Exemple 444 ------- 445 >>> import numpy as np 446 >>> import linear_operators as lo 447 >>> T = lo.TridiagonalOperator((3, 3), [1, 2, 3], [4, 5], [6, 7]) 448 >>> T.todense() 449 array([[1, 6, 0], 450 [4, 2, 7], 451 [0, 5, 3]]) 452 """ 453 if shape[0] != shape[1]: 454 raise ValueError("Only square operator is implemented.") 455 self.diag = diag 456 self.subdiag = subdiag 457 self.superdiag = superdiag 458 self.kwargs = kwargs 459 460 def matvec(x): 461 out = self.diag * x 462 out[:-1] += self.superdiag * x[1:] 463 out[1:] += self.subdiag * x[:-1] 464 return out

465 466 def rmatvec(x): 467 out = self.diag * x 468 out[:-1] += self.subdiag * x[1:] 469 out[1:] += self.superdiag * x[:-1] 470 return out

471 472 LinearOperator.__init__(self, shape, matvec, rmatvec=rmatvec, **kwargs) 473

474 - def __repr__(self):

475 s = LinearOperator.__repr__(self)[:-1] 476 s += ",\n superdiagonal=" + self.superdiag.__repr__() 477 s += ",\n diagonal=" + self.diag.__repr__() 478 s += ",\n subdiagonal=" + self.subdiag.__repr__() + ">" 479 return s

480

481 - def todense(self):

482 out = np.zeros(self.shape, dtype=self.dtype) 483 out += np.diag(self.diag) 484 out += np.diag(self.subdiag, -1) 485 out += np.diag(self.superdiag, 1) 486 return out

487

488 - def __getitem__(self, y):

489 # if tuple work on two dimensions 490 if isinstance(y, tuple): 491 # test dimension 492 if len(y) > 2: 493 raise IndexError("This is a 2-dimensional array.") 494 yi, yj = y 495 # single element case 496 if isinstance(yi, int) and isinstance(yj, int): 497 n = self.shape[0] 498 i, j = yi % n , yj % n 499 # outside 500 if np.abs(i - j) > 1: 501 return self.dtype(0) 502 # subdiag 503 elif i == j + 1: 504 # border case 505 if i == self.shape[0] - 1: 506 return self.dtype(0) 507 else: 508 return self.subdiag[i] 509 # superdiag 510 elif i == j - 1: 511 # border case 512 if i == self.shape[0]: 513 return self.dtype(0) 514 else: 515 return self.superdiag[i] 516 # diag 517 else: 518 return self.diag[i] 519 # case of tuple of length 1 520 elif len(y) == 1: 521 return self.__getitem__(self, y[0]) 522 # get a column 523 elif yi == slice(None, None) and isinstance(yj, int): 524 x = np.zeros(self.shape[1], dtype=self.dtype) 525 x[yj] = 1. 526 return self * x 527 # general case: no better way than todense 528 else: 529 d = self.todense() 530 return d[y] 531 # Work on lines : same cost as recasting to a dense matrix as 532 # all columns need to be accessed. 533 else: 534 d = self.todense() 535 return d[y]

536 537 @property

538 - def T(self):

539 kwargs = dict() 540 kwargs["dtype"] = getattr(self, 'dtype', None) 541 kwargs["dtypein"] = getattr(self, 'dtypein', None) 542 kwargs["dtypeout"] = getattr(self, 'dtypeout', None) 543 return TridiagonalOperator(self.shape, self.diag, self.superdiag, 544 self.subdiag, **kwargs)

545

546 - def toband(self):

547 """ 548 Convert the TridiagonalOperator into a BandOperator 549 """ 550 kl, ku = 1, 1 551 n = self.shape[1] 552 ab = np.zeros((kl + ku + 1, n)) 553 diags = (self.subdiag, self.diag, self.superdiag) 554 for i, d in zip((-1, 0, 1), diags): 555 ab[_band_diag(ku, i)] = d 556 return BandOperator(self.shape, ab, kl, ku, **self.kwargs)

557

558 # Multiple inheritence 559 -class SymmetricTridiagonal(SymmetricOperator, TridiagonalOperator):

560 - def __init__(self, shape, diag, subdiag, **kwargs):

561 """ 562 Parameters 563 ---------- 564 shape : length 2 tuple. 565 The shape of the operator. 566 567 diag : ndarray of size shape[0] 568 The diagonal of the matrix. 569 570 subdiag : ndarray of size shape[0] - 1 571 The subdiagonal and superdiagonal of the matrix. 572 573 Returns 574 ------- 575 A Tridiagonal matrix operator instance. 576 577 Exemple 578 ------- 579 >>> import numpy as np 580 >>> import linear_operators as lo 581 >>> T = lo.TridiagonalOperator((3, 3), [1, 2, 3], [4, 5]) 582 >>> T.todense() 583 array([[1, 4, 0], 584 [4, 2, 5], 585 [0, 5, 3]]) 586 """ 587 return TridiagonalOperator.__init__(self, shape, diag, subdiag, 588 subdiag, **kwargs)

589

590 - def toband(self):

591 """ 592 Convert into a SymmetricBandOperator. 593 """ 594 u = 2 # tridiagonal 595 n = self.shape[0] 596 # convert to ab format (lower) 597 ab = np.zeros((u, n)) 598 ab[0] = self.diag 599 ab[1, :-1] = self.subdiag 600 # options 601 kwargs = dict() 602 kwargs["dtype"] = getattr(self, 'dtype', None) 603 kwargs["dtypein"] = getattr(self, 'dtypein', None) 604 kwargs["dtypeout"] = getattr(self, 'dtypeout', None) 605 return SymmetricBandOperator(self.shape, ab, lower=True, **kwargs)

606

607 -class BandOperator(LinearOperator):

608 """ 609 Store a band matrix in ab format as defined in LAPACK 610 documentation. 611 612 a[i, j] is stored in ab[ku + 1 + i - j, j] 613 614 for max(1, j -ku) < i < min(m, j + kl) 615 616 Band storage of A (5, 5), kl = 2, ku = 1 : 617 618 * a01 a12 a23 a34 619 a00 a11 a22 a33 a44 620 a10 a21 a32 a43 * 621 a20 a31 a42 * * 622 623 Arguments 624 ---------- 625 shape : 2-tuple 626 Shape of the dense matrix equivalent. 627 kl : int 628 Number of subdiagonals 629 ku : int 630 Number of superdiagonals 631 632 Notes 633 ----- 634 For a description of band matrices see LAPACK doc : 635 636 http://www.netlib.org/lapack/lug/node124.html 637 638 """

639 - def __init__(self, shape, ab, kl, ku, **kwargs):

640 """ 641 Generate a BandOperator instance 642 643 Arguments 644 --------- 645 shape : 2-tuple 646 The shape of the operator 647 ab : ndarray with ndim == 2 648 Store the bands of the matrix using LAPACK storage scheme. 649 kl : int 650 Number of subdiagonals 651 ku : int 652 Number of superdiagonals 653 """ 654 if ab.shape[0] != kl + ku + 1 or ab.shape[1] != shape[1]: 655 raise ValueError("Wrong ab shape.") 656 657 self.ab = ab 658 self.kl = kl 659 self.ku = ku 660 self.kwargs = kwargs 661 662 def matvec(x): 663 # diag 664 out = self.ab[ku] * x 665 # upper part 666 for i in xrange(ku): 667 j = ku - i 668 out[:-j] += self.ab[i, j:] * x[j:] 669 for i in xrange(ku, kl + ku): 670 # lower part 671 out[i:] += self.ab[i + 1, :-i] * x[:-i] 672 return out

673 674 def rmatvec(x): 675 rab = self._rab 676 rkl, rku = ku, kl 677 # diag 678 out = rab[ku] * x 679 # upper part 680 for i in xrange(rku): 681 j = rku - i 682 out[:-j] += rab[i, j:] * x[j:] 683 for i in xrange(rku, rkl + rku): 684 # lower part 685 out[i:] += rab[i + 1, :-i] * x[:-i] 686 return out

687 688 return LinearOperator.__init__(self, shape, matvec, rmatvec, 689 **kwargs) 690 691 @property

692 - def T(self):

693 return BandOperator(self.shape[::-1], self.rab, self.ku, self.kl, 694 **self.kwargs)

695

696 - def diag(self, i=0):

697 """ 698 Returns the i-th diagonal (subdiagonal if i < 0, superdiagonal 699 if i >0). 700 """ 701 return self.ab[_band_diag(self.ku, i)]

702 703 @property

704 - def rab(self):

705 """ 706 Output the ab form of the transpose operator. 707 """ 708 ab = self.ab 709 kl, ku = self.kl, self.ku 710 rku, rkl = kl, ku 711 rab = np.zeros(ab.shape, dtype=ab.dtype) 712 for i in xrange(- kl, ku + 1): 713 rab[_band_diag(rku, -i)] = self.diag(i) 714 return rab

715

716 -def _band_diag(ku, i=0):

717 """ 718 Return a slice to get the i-th line of a band operator 719 """ 720 # diagonal 721 if i == 0: 722 return slice(ku, ku + 1) 723 # superdiagonal 724 if i > 0: 725 return (slice(ku - i, ku - i + 1, None), slice(i, None, None)) 726 # subdiagonal 727 if i < 0: 728 return (slice(ku - i, ku - i + 1, None), slice(None, i, None))

729

730 -class LowerTriangularOperator(BandOperator):

731 - def __init__(self, shape, ab, **kwargs):

732 kl = ab.shape[0] - 1 733 ku = 0 734 BandOperator.__init__(self, shape, ab, kl, ku, **kwargs)

735

736 -class UpperTriangularOperator(BandOperator):

737 - def __init__(self, shape, ab, **kwargs):

738 kl = 0 739 ku = ab.shape[0] - 1 740 BandOperator.__init__(self, shape, ab, kl, ku, **kwargs)

741

742 -class SymmetricBandOperator(SymmetricOperator):

743 - def __init__(self, shape, ab, lower=True, **kwargs):

744 if shape[0] != ab.shape[1]: 745 raise ValueError("ab.shape[1] should be equald to shape[0].") 746 if lower is False: 747 raise NotImplemented 748 749 self.ab = ab 750 self.lower = lower 751 self.kwargs = kwargs 752 753 def matvec(x): 754 out = self.ab[0] * x 755 for i in xrange(1, self.ab.shape[0]): 756 # upper part 757 out[:-i] += self.ab[i, :-i] * x[i:] 758 # lower part 759 out[i:] += self.ab[i, :-i] * x[:-i] 760 return out

761 762 return SymmetricOperator.__init__(self, shape, matvec, **kwargs)

763 764 @property

765 - def rab(self):

766 return self.ab

767

768 - def eigen(self, eigvals_only=False, overwrite_a_band=False, select='a', 769 select_range=None, max_ev=0):

770 """ 771 Solve real symmetric or complex hermitian band matrix 772 eigenvalue problem. 773 774 Uses scipy.linalg.eig_banded function. 775 """ 776 from scipy.linalg import eig_banded 777 778 w, v = eig_banded(self.ab, lower=self.lower, 779 eigvals_only=eigvals_only, 780 overwrite_a_band=overwrite_a_band, 781 select=select, 782 select_range=select_range, 783 max_ev=max_ev) 784 return EigendecompositionOperator(w=w, v=v)

785

786 - def cholesky(self, overwrite_ab=False):

787 """ 788 Chlesky decomposition. 789 Operator needs to be positive-definite. 790 791 Uses scipy.linalg.cholesky_banded. 792 793 Returns a matrix in ab form 794 """ 795 from scipy.linalg import cholesky_banded 796 797 ab_chol = cholesky_banded(self.ab, 798 overwrite_ab=overwrite_ab, 799 lower=self.lower) 800 if self.lower: 801 out = LowerTriangularOperator(self.shape, ab_chol, **self.kwargs) 802 else: 803 out = UpperTriangularOperator(self.shape, ab_chol, **self.kwargs) 804 return out

805

806 -class EigendecompositionOperator(SymmetricOperator):

807 """ 808 Define a SymmetricOperator from the eigendecomposition of another 809 SymmetricOperator. This can be used as an approximation for the 810 operator. 811 812 Inputs 813 ------- 814 815 A: LinearOperator (default: None) 816 The LinearOperator to approximate. 817 v: 2d ndarray (default: None) 818 The eigenvectors as given by arpack.eigsh 819 w: 1d ndarray (default: None) 820 The eigenvalues as given by arpack.eigsh 821 **kwargs: keyword arguments 822 Passed to the arpack.eigsh function. 823 824 You need to specify either A or v and w. 825 826 Returns 827 ------- 828 829 An EigendecompositionOperator instance, which is a subclass of the 830 SymmetricOperator. 831 832 Notes 833 ----- 834 835 This is really a wrapper for 836 scipy.sparse.linalg.eigen.arpack.eigsh 837 """

838 - def __init__(self, A=None, v=None, w=None, **kwargs):

839 from interface import aslinearoperator 840 kwargs['return_eigenvectors'] = True 841 if v is None or w is None: 842 w, v = eigsh(A, **kwargs) 843 shape = A.shape 844 dtype = A.dtype 845 dtypein = A.dtypein 846 dtypeout = A.dtypeout 847 else: 848 shape = 2 * (v.shape[0],) 849 dtype = v.dtype 850 dtypein = dtypeout = dtype 851 W = diagonal(w) 852 V = aslinearoperator(v) 853 M = V * W * V.T 854 # store some information 855 self.eigenvalues = w 856 self.eigenvectors = v 857 self.kwargs = kwargs 858 SymmetricOperator.__init__(self, shape, M.matvec, 859 dtypein=dtypein, 860 dtypeout=dtypeout, 861 dtype=dtype)

862 - def det(self):

863 """ 864 Output an approximation of the determinant from the 865 eigenvalues. 866 """ 867 return np.prod(self.eigenvalues)

868

869 - def logdet(self):

870 """ 871 Output the log of the determinant. Useful as the determinant 872 of large matrices can exceed floating point capabilities. 873 """ 874 return np.sum(np.log(self.eigenvalues))

875

876 - def __pow__(self, n):

877 """ 878 Raising an eigendecomposition to an integer power requires 879 only raising the eigenvalues to this power. 880 """ 881 return EigendecompositionOperator(v=self.eigenvectors, 882 w=self.eigenvalues ** n, 883 kwargs=self.kwargs)

884

885 - def inv(self):

886 """ 887 Returns the pseudo-inverse of the operator. 888 """ 889 return self ** -1

890

891 - def trace(self):

892 return np.sum(self.eigenvalues)

893

894 - def cond(self):

895 """ 896 Output an approximation of the condition number by taking the 897 ratio of the maximum over the minimum eigenvalues, removing 898 the zeros. 899 900 For better approximation of the condition number, one should 901 consider generating the eigendecomposition with the keyword 902 which='BE', in order to have a correct estimate of the small 903 eigenvalues. 904 """ 905 nze = self.eigenvalues[self.eigenvalues != 0] 906 return nze.max() / nze.min()

907

908 # functions 909 -def identity(shape, **kwargs):

910 """ 911 Parameters: 912 shape : the shape of the operator 913 Output: 914 I : identity LinearOperator 915 916 Exemple: 917 >>> I = identity((16, 16), dtype=float32) 918 """ 919 return IdentityOperator(shape, **kwargs)

920

921 -def homothetic(shape, coef, **kwargs):

922 return HomotheticOperator(shape, coef, **kwargs)

923

924 -def diagonal(d, **kwargs):

925 return DiagonalOperator(d, **kwargs)

926 927 # for backward compatibility 928 diag = diagonal

929 930 -def mask(mask, **kwargs):

931 return MaskOperator(mask, **kwargs)

932

933 -def shift(shape, s, **kwargs):

934 return ShiftOperator(shape, s, **kwargs)

935

936 -def permutation(p, **kwargs):

937 return PermutationOperator(p, **kwargs)

938

939 -def fft(shape, **kwargs):

940 return FftOperator(shape, **kwargs)

941

942 -def replication(shape, n, **kwargs):

943 return ReplicationOperator(shape, n, **kwargs)

944

945 -def slice_operator(shape, slice, **kwargs):

946 return SliceOperator(shape, slice, **kwargs)

947

948 -def tridiagonal(shape, diag, subdiag, superdiag, **kwargs):

949 return TridiagonalOperator(shape, diag, subdiag, superdiag, **kwargs)

950

951 -def symmetric_tridiagonal(shape, diag, subdiag, **kwargs):

952 return SymmetricTridiagonal(shape, diag, subdiag, **kwargs)

953

954 -def band_operator(shape, ab, kl, ku, **kwargs):

955 return BandOperator(shape, ab, kl, ku, **kwargs)

956

957 -def symmetric_band_operator(shape, ab, lower=True, **kwargs):

958 return SymmetricBandOperator(shape, ab, lower=lower, **kwargs)

959

960 # utilities 961 962 -def band_approximation(mat, kl=0, ku=0):

963 """ 964 Approximate a dense matrix as a BandOperator. 965 """ 966 ab = np.zeros((kl + ku + 1, mat.shape[1]), dtype=mat.dtype) 967 # diag 968 ab[ku] = np.diag(mat) 969 # upper 970 for i in xrange(ku): 971 j = ku - i 972 ab[i, j:] = np.diag(mat, j) 973 # lower 974 for i in xrange(ku + 1, kl + ku + 1): 975 j = ku - i 976 ab[i, :j] = np.diag(mat, j) 977 return BandOperator(mat.shape, ab, kl, ku, dtype=mat.dtype)

978

979 -def symmetric_band_approximation(mat, k=0):

980 """ 981 Approximate a symmetrix matrix as a SymmetricBandOperator. 982 """ 983 ab = np.zeros((k + 1, mat.shape[1]), dtype=mat.dtype) 984 ab[0] = np.diag(mat) 985 for i in xrange(1, k + 1): 986 ab[i, :-i] = np.diag(mat, i) 987 return SymmetricBandOperator(mat.shape, ab, dtype=mat.dtype)

988

Source Code for Module linear_operators.operators