Mercurial > repos > pointsource_pde / file comparison
comparison: src/compose.py
src/compose.py
- changeset 3
- c3a4f4bb87f7
- parent 1
- a4137aedcb3a
equal
deleted
inserted
replaced
| 1:a4137aedcb3a | 3:c3a4f4bb87f7 |
|---|---|
| 30 res = 0 | 30 res = 0 |
| 31 for f_i in self.fnlist: | 31 for f_i in self.fnlist: |
| 32 res = max(res, f_i.diff_lipschitz_factor()) | 32 res = max(res, f_i.diff_lipschitz_factor()) |
| 33 return res | 33 return res |
| 34 | 34 |
| 35 def diff_bound(self, xbound=None): | 35 def diff_bound(self): |
| 36 res = 0 | 36 res = 0 |
| 37 for f_i in self.fnlist: | 37 for f_i in self.fnlist: |
| 38 res = max(res, f_i.diff_bound(xbound=xbound)) | 38 res = max(res, f_i.diff_bound()) |
| 39 return res | 39 return res |
| 40 | 40 |
| 41 | 41 |
| 42 class ComposeFnWithOperator: | 42 class ComposeFnWithOperator: |
| 43 def __init__(self, f, op, xbound=None, xbound_pair=None): | 43 def __init__(self, f, op): |
| 44 self.f = f | 44 self.f = f |
| 45 self.op = op | 45 self.op = op |
| 46 self.xbound = xbound | |
| 47 self.xbound_pair = xbound_pair | |
| 48 | 46 |
| 49 def apply(self, *args): | 47 def apply(self, *args): |
| 50 return self.f.apply(self.op.apply(*args)) | 48 return self.f.apply(self.op.apply(*args)) |
| 51 | 49 |
| 52 def diff(self, *args): | 50 def diff(self, *args): |
| 59 # TODO: precalculations in apply should be used in diff_adjdir | 57 # TODO: precalculations in apply should be used in diff_adjdir |
| 60 w = self.op.apply(*args) | 58 w = self.op.apply(*args) |
| 61 (a, v) = self.f.apply_and_diff(w) | 59 (a, v) = self.f.apply_and_diff(w) |
| 62 return (a, self.op.diff_adjdir(v, *args, apply_result=w)) | 60 return (a, self.op.diff_adjdir(v, *args, apply_result=w)) |
| 63 | 61 |
| 62 def diff_bound(self): | |
| 63 mf = self.f.diff_bound() | |
| 64 if hasattr(self.op, "opnorm"): | |
| 65 lda = self.op.opnorm() ** 2 | |
| 66 else: | |
| 67 lda = self.op.diff_bound() | |
| 68 return lda * mf | |
| 69 | |
| 70 def diff_bound_pair(self): | |
| 71 mf = self.f.diff_bound() | |
| 72 lda1, lda2 = self.op.diff_bound_pair() | |
| 73 return lda1 * mf, lda2 * mf | |
| 74 | |
| 75 # def lipschitz_factor(self, xbound=None): | |
| 76 # if xbound is None: | |
| 77 # xbound = self.xbound | |
| 78 # lf = self.f.lipschitz_factor(xbound=self.op.codomain_bound(xbound=xbound)) | |
| 79 # la = self.op.lipschitz_factor(xbound=xbound) | |
| 80 # return lf * la | |
| 81 | |
| 82 # def lipschitz_factor_pair(self, xbound=None): | |
| 83 # if xbound is None: | |
| 84 # xbound = self.xbound | |
| 85 # lf = self.f.lipschitz_factor(xbound=self.op.codomain_bound(xbound=xbound)) | |
| 86 # la1, la2 = self.op.lipschitz_factor_pair(xbound=xbound) | |
| 87 # return lf * la1, lf * la2 | |
| 88 | |
| 64 def diff_lipschitz_factor(self): | 89 def diff_lipschitz_factor(self): |
| 65 # ‖∇A(x)^*∇F(A(x)) - ∇A(y)^*∇F(A(y))‖ | 90 """ |
| 66 # = ‖[∇A(x)^*-∇A(y)^*]∇F(A(x)) - ∇A(y)^*[∇F(A(y))-∇F(A(x))]‖ | 91 Calculate the Lipschitz factor of the differential of this composed function. |
| 67 # ≤ L_{∇A(x)} M_{∇F} + M_{∇A(y)^*} L_{∇F}L_A. | 92 We assume that either the operator is linear and implementes `opnorm`, or it is nonlinear, and impliements `diff_chain_lipschitz_factor` to directly calculate a Lipschitz factor of $x ↦ ∇A(x)^*∇F(A(x))$, given a bound $M$ |
| 93 on ∇F. The function `diff_chain_lipschitz_factor` should return the Lipschitz factor divided by $M$: we obtain $M$ through the `diff_bound` on $F$. | |
| 94 """ | |
| 95 | |
| 68 if hasattr(self.op, "opnorm"): | 96 if hasattr(self.op, "opnorm"): |
| 69 # Linear operator | 97 return self.f.diff_lipschitz_factor() * self.op.opnorm() ** 2 |
| 70 lda = 0.0 # This is zero, | |
| 71 mdf = 0.0 # hence this not needed. | |
| 72 else: | 98 else: |
| 73 mdf = self.f.diff_bound(xbound=self.op.codomain_bound(xbound=self.xbound)) | 99 mdf = self.f.diff_bound() |
| 74 lda = self.op.diff_adj_lipschitz_factor() | 100 lda = self.op.diff_chain_lipschitz_factor() |
| 75 | 101 # print( |
| 76 ldf = self.f.diff_lipschitz_factor() | 102 # "LDA %s %f; MDF %f; total %f" |
| 77 la = self.op.lipschitz_factor() | 103 # % (type(self.op).__name__, lda, mdf, mdf * lda), |
| 78 mda = self.op.diff_bound(xbound=self.xbound) | 104 # ) |
| 79 | 105 return mdf * lda |
| 80 return lda * mdf + mda * ldf * la | |
| 81 | 106 |
| 82 def diff_lipschitz_factor_pair(self): | 107 def diff_lipschitz_factor_pair(self): |
| 83 if self.op.hasattr("opnorm"): | 108 """ |
| 84 # Linear operator | 109 This is similar to `diff_lipschitz_factor`, except separates the |
| 85 lda1, lda2 = 0.0, 0.0 # This is zero, | 110 factor for arguments pairs. |
| 86 mdf = 0.0 # hence this not needed. | |
| 87 else: | |
| 88 lda1, lda2 = self.op.diff_adj_lipschitz_factor_pair() | |
| 89 mdf = self.f.diff_bound( | |
| 90 xbound=self.op.codomain_bound_pair(xbound=self.xbound_pair) | |
| 91 ) | |
| 92 | 111 |
| 93 ldf = self.f.diff_lipschitz_factor() | 112 This requires the operator to implement `diff_chain_lipschitz_factor`; |
| 94 la1, la2 = self.op.lipschitz_factor_pair() | 113 there is no special handling of linear opeartors. |
| 95 mda = self.op.diff_bound_pair(xbound=self.xbound_pair) | 114 """ |
| 96 | 115 |
| 97 return lda1 * mdf + mda * ldf * la1, lda2 * mdf + mda * ldf * la2 | 116 mdf = self.f.diff_bound() |
| 117 | |
| 118 lda1, lda2 = self.op.diff_chain_lipschitz_factor_pair() | |
| 119 | |
| 120 # print("LDA %s %f %f; MDF %f; " % (type(self.op).__name__, lda1, lda2, mdf)) | |
| 121 | |
| 122 return mdf * lda1, mdf * lda2 | |
| 98 | 123 |
| 99 | 124 |
| 100 class InjectSecond: | 125 class InjectSecond: |
| 101 def __init__(self, y): | 126 def __init__(self, y): |
| 102 self.y = y | 127 self.y = y |
| 110 # This is not really a linear operator, but for our purposes affine behave essentially | 135 # This is not really a linear operator, but for our purposes affine behave essentially |
| 111 # the same | 136 # the same |
| 112 def opnorm(self, *args): | 137 def opnorm(self, *args): |
| 113 return 1.0 | 138 return 1.0 |
| 114 | 139 |
| 115 def lipschitz_factor(self, *args): | 140 # def lipschitz_factor(self): |
| 141 # return 1.0 | |
| 142 | |
| 143 # def diff_chain_lipschitz_factor(self, *args): | |
| 144 # return 0.0 | |
| 145 | |
| 146 def diff_bound(self): | |
| 116 return 1.0 | 147 return 1.0 |
| 117 | 148 |
| 118 def diff_adj_lipschitz_factor(self, *args): | 149 # def codomain_bound(self, xbound=None): |
| 119 return 0.0 | 150 # if xbound is None: |
| 120 | 151 # raise Exception("Linear operators have unbounded range") |
| 121 def diff_bound(self, xbound=None): | 152 # else: |
| 122 return 1.0 | 153 # return xbound |
| 123 | |
| 124 def codomain_bound(self, xbound=None): | |
| 125 if xbound is None: | |
| 126 raise Exception("Linear operators have unbounded range") | |
| 127 else: | |
| 128 return xbound |
