SimpleAlgortihm.py 10.2 KB
 jansen31 committed Apr 21, 2020 1 2 ``````# Simple Algorithm # Analog to AlgorithmWrapper in ProxMatlab `````` 3 `````` `````` Matthijs committed Apr 22, 2020 4 ``````from numpy import zeros, angle, trace, exp, sqrt, sum, ndarray `````` 5 6 ``````from numpy.linalg import norm `````` jansen31 committed May 10, 2019 7 `````` `````` Matthijs committed Apr 22, 2020 8 9 10 11 ``````def nd_norm_difference(arr1: ndarray, arr2: ndarray, norm_factor: float = 1.): """ Root-sum-square norm of the difference of 2 numpy arrays, equal to Frobenius norm for arr.ndim == 2 """ `````` Matthijs committed Apr 23, 2020 12 `````` return sqrt(sum(abs(arr1 - arr2) ** 2)) / norm_factor `````` Matthijs committed Apr 22, 2020 13 14 15 16 17 18 19 20 21 22 23 24 25 `````` def phase_offset_compensated_norm(arr1: ndarray, arr2: ndarray, norm_factor: float = 1., norm_type: str = 'fro') -> float: """ Norm of the difference between 2 numpy arrays. :param arr1: numpy array :param arr2: numpy array like arr1 :param norm_factor: normalization factor :param norm_type: for 1d and 2d arrays, pass norm_type to numpy.linalg.norm :return: float """ `````` Matthijs committed Apr 21, 2020 26 27 28 29 30 `````` if arr1.ndim < 3: phase_offset = angle(sum(arr1 * arr2.conj())) return norm(arr1 - arr2 * exp(1j * phase_offset), norm_type) / norm_factor else: phase_offset = angle(sum(arr1 * arr2.conj())) `````` Matthijs committed Apr 22, 2020 31 `````` return nd_norm_difference(arr1, arr2 * exp(1j * phase_offset), norm_factor=norm_factor) `````` jansen31 committed Jun 12, 2019 32 33 `````` `````` 34 35 ``````class SimpleAlgorithm: `````` jansen31 committed May 10, 2019 36 `````` def __init__(self, config): `````` 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 `````` """ Parameters ---------- config : dict Dictionary containing the problem configuration. It must contain the following mappings: proxoperators: 2 ProxOperators Tuple of ProxOperators (the class, no instance) beta_0: number Starting relaxation parmater beta_max: number Maximum relaxation parameter beta_switch: int Iteration at which beta moves from beta_0 -> beta_max `````` 52 `````` norm_data: number `````` 53 54 55 56 57 58 59 60 61 62 63 64 65 `````` ? Nx: int Row-dim of the product space elements Ny: int Column-dim of the product space elements Nz: int Depth-dim of the product space elements dim: int Size of the product space """ self.prox1 = config['proxoperators'][0](config) self.prox2 = config['proxoperators'][1](config) `````` 66 `````` self.norm_data = config['norm_data'] `````` jansen31 committed Jun 12, 2019 67 68 `````` self.Nx = config['Nx'] self.Ny = config['Ny'] `````` jansen31 committed May 10, 2019 69 `````` self.Nz = config['Nz'] `````` 70 71 72 73 74 `````` self.product_space_dimension = config['product_space_dimension'] self.iter = 0 self.config = config if 'truth' in config: `````` markus.meier01 committed Dec 15, 2019 75 76 77 `````` self.truth = config['truth'] self.truth_dim = config['truth_dim'] self.norm_truth = config['norm_truth'] `````` 78 79 `````` if 'diagnostic' in config: `````` Matthijs committed Apr 21, 2020 80 `````` self.diagnostic = config['diagnostic'] `````` 81 82 83 `````` else: self.diagnostic = False `````` jansen31 committed Jun 12, 2019 84 85 86 87 88 89 90 91 `````` def evaluate(self, u): """ Method to override in specific algorithm subclass :param u: old iterate :return: new iterate """ return u `````` 92 93 94 `````` def run(self, u, tol, maxiter): """ Runs the algorithm for the specified input data `````` Matthijs committed Apr 22, 2020 95 96 97 98 99 `````` :param u: iterate :param tol: tolerance for convergence, if change= tol: iter += 1 `````` jansen31 committed May 10, 2019 141 `````` config['iter'] = iter `````` Matthijs committed Apr 23, 2020 142 143 `````` if 'tqdm_progressbar' in self.config and self.config['tqdm_progressbar'] is not None: pbar.update() `````` jansen31 committed May 10, 2019 144 145 `````` # makes call to algorithm: config['u'] = u `````` 146 147 `````` u_new = self.evaluate(u) `````` Matthijs committed Apr 21, 2020 148 `````` if self.diagnostic: `````` Matthijs committed Apr 22, 2020 149 150 `````` # the next prox operation only gets used in the computation of the size of the gap. # This extra step is not required in alternating projections, which makes RAAR a bit slower `````` 151 152 153 `````` u2 = self.prox2.work(u_new) u1 = self.prox1.work(u2) `````` Matthijs committed Apr 22, 2020 154 `````` # compute the normalized change in successive iterates `````` jansen31 committed May 10, 2019 155 156 157 `````` tmp_change = 0 tmp_gap = 0 tmp_shadow = 0 `````` 158 159 `````` # the simple case, where everything can be calculated in one difference if self.product_space_dimension == 1 or (p == 1 and q == 1): `````` jansen31 committed Jun 12, 2019 160 `````` tmp_change = phase_offset_compensated_norm(u, u_new, norm_type='fro', norm_factor=norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 161 `````` if self.diagnostic: `````` 162 163 164 `````` # For Douglas-Rachford,in general it is appropriate to monitor the SHADOWS of the iterates, # since in the convex case these converge even for beta=1. (see Bauschke-Combettes-Luke, # J. Approx. Theory, 2004) `````` jansen31 committed Jun 12, 2019 165 166 `````` tmp_shadow = phase_offset_compensated_norm(u2, shadow, norm_factor=norm_data, norm_type='fro') ** 2 tmp_gap = phase_offset_compensated_norm(u1, u2, norm_factor=norm_data, norm_type='fro') ** 2 `````` 167 `````` if hasattr(self, 'truth'): `````` markus.meier01 committed Dec 15, 2019 168 `````` if self.truth_dim[0] == 1: `````` jansen31 committed May 10, 2019 169 `````` z = u1[0, :] `````` markus.meier01 committed Dec 15, 2019 170 `````` elif self.truth_dim[1] == 1: `````` jansen31 committed May 10, 2019 171 172 173 `````` z = u1[:, 0] else: z = u1 `````` markus.meier01 committed Dec 15, 2019 174 `````` Relerrs[iter] = phase_offset_compensated_norm(self.truth, z, norm_factor=self.norm_truth, `````` jansen31 committed Jun 12, 2019 175 `````` norm_type='fro') `````` jansen31 committed May 10, 2019 176 `````` `````` 177 `````` elif q == 1: # more complex: loop over product space dimension `````` Matthijs committed Apr 22, 2020 178 `````` for j in range(self.product_space_dimension): # this loop might be avoided using the 3d norm function `````` jansen31 committed May 10, 2019 179 `````` tmp_change = tmp_change + (norm(u[:, :, j] - u_new[:, :, j], 'fro') / norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 180 `````` if self.diagnostic: `````` 181 `````` # compute (||P_SP_Mx-P_Mx||/norm_data)^2: `````` jansen31 committed May 10, 2019 182 183 `````` tmp_gap = tmp_gap + (norm(u1[:, :, j] - u2[:, :, j], 'fro') / norm_data) ** 2 tmp_shadow = tmp_shadow + (norm(u2[:, :, j] - shadow[:, :, j], 'fro') / norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 184 `````` if self.diagnostic: `````` 185 `````` if hasattr(self, 'truth'): `````` jansen31 committed May 10, 2019 186 `````` z = u1[:, :, 0] `````` 187 188 `````` Relerrs[iter] = norm((self.truth - exp(-1j * angle(trace(self.truth.T.transpose() * z))) * z), 'fro') / self.norm_truth `````` markus.meier01 committed Dec 15, 2019 189 `````` `````` 190 `````` else: # loop over product space dimension as well as additional z dimension `````` Matthijs committed Apr 21, 2020 191 `````` if self.diagnostic: `````` alexander.dornheim committed Dec 13, 2017 192 193 `````` if hasattr(self, 'truth'): Relerrs[iter] = 0 `````` Matthijs committed Apr 22, 2020 194 195 `````` for j in range(self.product_space_dimension): # this loop might be avoided using a 4d norm function for k in range(self.Nz): # this loop might be avoided using the 3d norm function `````` jansen31 committed May 10, 2019 196 `````` tmp_change = tmp_change + (norm(u[:, :, k, j] - u_new[:, :, k, j], 'fro') / norm_data) ** 2 `````` Matthijs committed Apr 21, 2020 197 `````` if self.diagnostic: `````` jansen31 committed May 10, 2019 198 199 200 `````` # compute (||P_Sx-P_Mx||/norm_data)^2: tmp_gap = tmp_gap + (norm(u1[:, :, k, j] - u2[:, :, k, j], 'fro') / (norm_data)) ** 2 tmp_shadow = tmp_shadow + ( `````` jansen31 committed Jun 12, 2019 201 `````` norm(u2[:, :, k, j] - shadow[:, :, k, j], 'fro') / (norm_data)) ** 2 `````` jansen31 committed May 10, 2019 202 203 `````` if hasattr(self, 'truth') and (j == 0): Relerrs[iter] = Relerrs[iter] + norm( `````` markus.meier01 committed Dec 15, 2019 204 `````` self.truth - exp(-1j * angle(trace(self.truth.T * u1[:, :, k, 1]))) * u1[:, :, k, 1], `````` markus.meier01 committed Dec 10, 2019 205 `````` 'fro') / self.norm_Truth `````` 206 `````` `````` 207 `````` # Add values to the list of change, gap and shadow_change `````` 208 `````` change[iter] = sqrt(tmp_change) `````` Matthijs committed Apr 21, 2020 209 `````` if self.diagnostic: `````` 210 `````` gap[iter] = sqrt(tmp_gap) `````` jansen31 committed May 10, 2019 211 `````` shadow_change[iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to `````` 212 213 214 `````` # the unregularized set. To monitor the Euclidean norm of the gap to the regularized set is # expensive to calculate, so we use this surrogate. Since the stopping criteria is on the change # in the iterates, this does not matter. `````` 215 `````` `````` 216 `````` # update iterate `````` jansen31 committed May 10, 2019 217 `````` u = u_new `````` Matthijs committed Apr 21, 2020 218 `````` if self.diagnostic: `````` Matthijs committed Apr 22, 2020 219 220 `````` # For Douglas-Rachford, in general it is appropriate to monitor the SHADOWS of the iterates, since in # the convex case these converge even for beta=1. (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004) `````` jansen31 committed May 10, 2019 221 `````` shadow = u2 `````` 222 223 `````` ##### POSTPROCESSING `````` jansen31 committed May 10, 2019 224 225 `````` u2 = self.prox2.work(u) u1 = self.prox1.work(u2) `````` 226 `````` if self.Nx == 1: `````` jansen31 committed May 10, 2019 227 228 `````` u1 = u1[:, 0] u2 = u2[:, 0] `````` 229 `````` elif self.Ny == 1: `````` jansen31 committed May 10, 2019 230 231 `````` u1 = u1[0, :] u2 = u2[0, :] `````` 232 `````` elif self.Nz == 1 and u1.ndim > 2: `````` jansen31 committed May 10, 2019 233 234 `````` u1 = u1[:, :, 0] u2 = u2[:, :, 0] `````` 235 `````` `````` jansen31 committed May 10, 2019 236 237 `````` change = change[1:iter + 1] output = {'u': u, 'u1': u1, 'u2': u2, 'iter': iter, 'change': change} `````` Matthijs committed Apr 21, 2020 238 `````` if self.diagnostic: `````` jansen31 committed May 10, 2019 239 240 241 242 `````` gap = gap[1:iter + 1] shadow_change = shadow_change[1:iter + 1] output['gap'] = gap output['shadow_change'] = shadow_change `````` 243 `````` if hasattr(self, 'truth'): `````` jansen31 committed May 10, 2019 244 `````` output['Relerrs'] = Relerrs `````` 245 `` return output``