Mercurial > repos > pointsource_pde / file comparison
comparison: src/quadratic_dataterm.py
src/quadratic_dataterm.py
- changeset 1
- a4137aedcb3a
- child 3
- c3a4f4bb87f7
equal
deleted
inserted
replaced
| 0:7ec1cfe19a24 | 1:a4137aedcb3a |
|---|---|
| 1 import numpy as np | |
| 2 from numpy.linalg import norm | |
| 3 | |
| 4 """ | |
| 5 A quadratic data term | |
| 6 """ | |
| 7 | |
| 8 | |
| 9 class QuadraticDataTerm: | |
| 10 def __init__(self, opA, b, b_norm, xbound=None, λ=1.0): | |
| 11 self.opA = opA | |
| 12 self.b = b | |
| 13 self.b_norm = b_norm | |
| 14 self.xbound = xbound | |
| 15 self.λ = λ | |
| 16 | |
| 17 def apply(self, x): | |
| 18 v = self.opA.apply(x) - self.b | |
| 19 return norm(v) ** 2 * (self.λ / 2) | |
| 20 | |
| 21 def diff(self, x): | |
| 22 v = self.opA.apply(x) - self.b | |
| 23 j = self.λ * v | |
| 24 return self.opA.diff_adjdir(j, x) | |
| 25 | |
| 26 def apply_and_diff(self, x): | |
| 27 v = self.opA.apply(x) - self.b | |
| 28 val = norm(v) ** 2 * (self.λ / 2) | |
| 29 j = self.λ * v | |
| 30 d = self.opA.diff_adjdir(j, x) | |
| 31 return (val, d) | |
| 32 | |
| 33 def diff_lipschitz_factor(self, xbound=None): | |
| 34 if hasattr(self.opA, "opnorm"): | |
| 35 # Linear operator | |
| 36 return self.λ * self.opA.opnorm() ** 2 | |
| 37 else: | |
| 38 raise Exception("diff_lipschitz_factor for non nonlinear operators") | |
| 39 # ‖∇A(x)^*∇F(A(x)) - ∇A(y)^*∇F(A(y))‖ | |
| 40 # = ‖[∇A(x)^*-∇A(y)^*]∇F(A(x)) - ∇A(y)^*[A(y)-A(x)]‖ | |
| 41 # ≤ L_{∇A(x)}(M_A+M_b) + M_{∇A(y)^*} L_A. | |
| 42 la = self.opA.lipschitz_factor() | |
| 43 ma = self.opA.bound(xbound=xbound) | |
| 44 mb = self.b.norm() | |
| 45 mda = self.opA.diff_bound(xbound=xbound) | |
| 46 lda = self.opA.diff_adj_lipschitz_factor() | |
| 47 return self.λ * (lda * (ma + mb) + mda * la) | |
| 48 | |
| 49 # def diff_lipschitz_factor_pair(self, *args): | |
| 50 # return self.diff_lipschitz_factor() | |
| 51 | |
| 52 def diff_bound(self, xbound=None): | |
| 53 return self.λ * (self.opA.codomain_bound(xbound=xbound) + self.b_norm) |
