from proxtoolbox.algorithms.iterateMonitor import IterateMonitor
from proxtoolbox.utils.cell import Cell, isCell
import numpy as np
from numpy import zeros, angle, trace, exp, sqrt, sum, matmul, array, reshape
from numpy.linalg import norm
[docs]class FeasibilityIterateMonitor(IterateMonitor):
"""
Algorithm analyzer for monitoring iterates of
projection algorithms for feasibility problems.
Specialization of the IterateMonitor class.
"""
def __init__(self, experiment):
super(FeasibilityIterateMonitor, self).__init__(experiment)
self.gaps = None
self.shadow_changes = None
self.rel_errors = None
# instantiate prox operators
self.proxOperators = []
for proxClass in experiment.proxOperators:
if proxClass != '':
prox = proxClass(experiment)
self.proxOperators.append(prox)
[docs] def preprocess(self, alg):
"""
Allocate data structures needed to collect
statistics. Called before running the algorithm.
Parameters
----------
alg : instance of Algorithm class
Algorithm to be monitored.
"""
super(FeasibilityIterateMonitor, self).preprocess(alg)
nProx = len(self.proxOperators)
# Even if more detailed diagnostics are not run, the mathematical
# exit criteria will be on step-sizes, for which one needs to monitor the
# norm of the difference between successive iterates of the blocks, in
# the case that u is a cell of blocks.
# To do this, we need to convert the relevant block elements into ARRAYS
if self.diagnostic:
# if the algorithm is cyclic and one still wants to do a full
# diagnostic, computing the generalized difference vector
# (see Luke, Nguyen&Tam, Math of Oper. Res., 2017/2018) is
# expensive. Need to expand u to the product space. Use cells for this
self.u_monitor = Cell(nProx)
# assumes that u0 is an array
self.u_monitor[0] = self.u0
for j in range(1, nProx):
self.u_monitor[j] = self.proxOperators[j].eval(self.u_monitor[j-1], j)
else: # always creates a u_monitor, issue is whether to monitor the intermediate
# steps for the gap, as the above prepares for in the case that
# diagnostic ist true.
self.u_monitor = self.u0
# set up diagnostic arrays
if self.diagnostic:
self.gaps = self.changes.copy()
self.gaps[0] = sqrt(999)
self.shadow_changes = self.changes.copy()
self.shadow_changes[0] = sqrt(999)
if self.truth is not None:
self.rel_errors = self.changes.copy()
self.rel_errors[0] = sqrt(999)
[docs] def updateStatistics(self, alg):
"""
Update statistics. Called at each iteration
while the algorithm is running.
Parameters
----------
alg : instance of Algorithm class
Algorithm being monitored.
Notes
-----
The call to the Numpy norm function should not
specify 'fro' since this is the default. This
would not be an issue except for the fact that
this norm function throws an exception when we apply
it to a vector, that is 'fro' is not considered
a valid option for a vector.
"""
normM = self.norm_data
nProx = len(self.proxOperators)
u = alg.u_new
prev_u_mon = self.u_monitor.copy()
self.u_monitor[0] = u
if self.rotate:
# the iterates of algorithms like cyclic projections for phase retrieval do not (in general) converge to
# an absolute fixed point, but rather a fixed point relative to
# a global phase shift. Our observation in the case of cyclic projections is that
# the iterates converge to the correct relative phases, but on each iteration the
# global phase shifts by a constant amount, making it look like the algorithm
# is not converging to something reasonable. We don't know how to characterize this
# yet (as of 2018) but a rotation of the iterates to a fixed global reference
# reveals the RELEVANT convergence of the (shadow sequence of the) algorithm:
if isCell(u):
for j in range(len(u)):
self.u_monitor[0][j] = exp(-1j*angle(trace(matmul(prev_u_mon[0][j].T.conj(), u[j]))))*u[j]
else:
# assumes that u is a ndarray
self.u_monitor[0] = exp(-1j*angle(trace(matmul(prev_u_mon[0].T.conj(), u))))*u
u1 = self.proxOperators[nProx-1].eval(self.u_monitor[0])
# for Douglas Rachford and DRlambda the iterates themselves do not
# converge to the intersection, even if this is nonempty. The Shadows, however,
# do go to the intersection, if this is nonempty.
if self.diagnostic:
self.u_monitor[1] = self.proxOperators[0].eval(u1, 0)
for j in range(1, nProx-1):
self.u_monitor[j+1] = self.proxOperators[j].eval(self.u_monitor[j], j)
# convert u1 to an array if necessary
if isCell(u1):
p, q = self.getIterateSize(u1[len(u1)-1])
#[~,~,p,q]=size(u1{length(u1)});
else:
p, q = self.getIterateSize(u1)
#[~,~,p,q]=size(u1);
# compute the normalized change in successive iterates:
# change(iter) = sum(sum((feval('P_M',M,u)-tmp).^2))/normM;
tmp_change = 0
tmp_gap = 0
tmp_shadow = 0
if p == 1 and q == 1:
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
tmp_change = tmp_change + (norm(self.u_monitor[0][j] - prev_u_mon[0][j])/normM)**2
else:
tmp_change = (norm(self.u_monitor[0] - prev_u_mon[0])/normM)**2
if self.diagnostic:
tmp_shadow = tmp_change
# For Douglas-Rachford,in general it is appropriate to monitor the
# SHADOWS of the iterates, since in the convex case these converge
# even for lambda=1.
# (see Bauschke-Combettes-Luke, J. Approx. Theory, 2004)
if isCell(u1):
for jj in range(len(u1)):
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1][jj] - prev_u_mon[1][jj])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[1][jj] - u1[jj])/normM)**2
else:
tmp_shadow = tmp_shadow + (norm(self.u_monitor[1] - prev_u_mon[1])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[1] - u1)/normM)**2
for j in range(2, nProx):
# 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)
if isCell(self.u_monitor[j]):
for jj in range(len(self.u_monitor[j])):
tmp_shadow = tmp_shadow + (norm(self.u_monitor[j][jj] - prev_u_mon[j][jj])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[j][jj] - self.u_monitor[j-1][jj])/normM)**2
else:
tmp_shadow = tmp_shadow + (norm(self.u_monitor[j] - prev_u_mon[j])/normM)**2
tmp_gap = tmp_gap + (norm(self.u_monitor[j] - self.u_monitor[j-1])/normM)**2
if self.truth is not None:
# convert u1 to an array if necessary
if isCell(u1):
z = u1[len(u1)-1]
else:
z = u1
if self.truth_dim[0] == 1:
z = z[0,:]
z = reshape(z, (1,len(z))) # we want a true matrix not just a vector
elif self.truth_dim[1] == 1:
z = z[:,0]
z = reshape(z, (len(z),1)) # we want a true matrix not just a vector
"""
# TODO deal properly with the following code which
# is experiment-dependent
if(any(strcmp('experiment', fieldnames(method_input))))
if(strcmp(method_input.experiment,'Sparse1'))
% only looks to see if the support of the signal has been found
% up to linear shifts and reflections.
[~,Iz]=sort(z, 'descend');
Iz=Iz(1:length(method_input.truth_supp_diff));
Iz=sort(Iz,'ascend');
z_supp_diff=[Iz(2:end);Iz(1)+method_input.Nx*method_input.Ny]-Iz;
Z=circulant(z_supp_diff);
Z2=circulant(z_supp_diff');
Relerrs=min(sqrt(min(sum(abs(Z-method_input.truth_supp_diff).^2,1))),sqrt(min(sum(abs(Z2-method_input.truth_supp_diff).^2,1))));
elseif(strcmp(method_input.experiment,'Sparse2'))
% only looks to see if the support of the signal has been found
% up to linear shifts and reflections.
tmp=method_input.sparsity_parameter;
method_input.sparsity_parameter=method_input.true_sparsity;
z=feval(method_input.Prox{1}, method_input,z);
method_input.sparsity_parameter=tmp;
Iz=z~=0;
z(Iz)=z(Iz)./z(Iz);
zcov=fftshift(real(ifft2(fft2(z).*ifft2(z)))./method_input.Nx);
Relerrs=norm(zcov-method_input.cov_true_support);
"""
rel_error = norm(self.truth - exp(-1j*angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
elif q == 1:
if isCell(self.u_monitor[0]):
for j in range(len(self.u_monitor[0])):
for jj in range(p):
tmp_change = tmp_change + (norm(self.u_monitor[0][j][:,:,jj] - prev_u_mon[0][j][:,:,jj])/normM)**2
if self.diagnostic:
tmp_shadow = tmp_change
for k in range(1, nProx):
# compute (||P_SP_Mx-P_Mx||/normM)^2:
for j in range(len(self.u_monitor[0])):
for jj in range(p):
tmp_gap = tmp_gap + (norm(self.u_monitor[k][j][:,:,jj] - self.u_monitor[k-1][:,:,jj])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[k][j][:,:,jj] - prev_u_mon[k][j][:,:,jj])/normM)**2
else:
for j in range(p):
tmp_change = tmp_change + (norm(self.u_monitor[0][:,:,j] - prev_u_mon[0][:,:,j])/normM)**2
if self.diagnostic:
tmp_shadow = tmp_change
for k in range(1, nProx):
# compute (||P_SP_Mx-P_Mx||/normM)^2:
tmp_gap = tmp_gap + (norm(self.u_monitor[k][:,:,j] - self.u_monitor[k-1][:,:,j])/normM)**2
tmp_shadow = tmp_shadow + (norm(self.u_monitor[k][:,:,j] - prev_u_mon[k][:,:,j])/normM)**2
if self.diagnostic and self.truth is not None:
z = u1[:,:,0]
rel_error = norm(self.truth - exp(-1j*angle(trace(matmul(self.truth.T.conj(), z)))) * z) / self.norm_truth
else: # cells of 4D arrays???
pass
self.changes[alg.iter] = sqrt(tmp_change)
if self.diagnostic:
self.gaps[alg.iter] = sqrt(tmp_gap)
self.shadow_changes[alg.iter] = sqrt(tmp_shadow) # this is the Euclidean norm of the gap to
# 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.
if self.truth is not None:
self.rel_errors[alg.iter] = rel_error
[docs] def postprocess(self, alg, output):
"""
Called after the algorithm stops. Store statistics in
the given 'output' dictionary
Parameters
----------
alg : instance of Algorithm class
Algorithm that was monitored.
output : dictionary
Contains the last iterate and various statistics that
were collected while the algorithm was running.
Returns
-------
output : dictionary into which the following entries are added
(when diagnostics are required)
gaps : ndarray
Squared gap distance normalized by the magnitude
constraint
shadow_changes : ndarray
Euclidean norm of the gap to the unregularized set.
To monitor the Euclidean norm of the gap to the
regularized set is expensive to calculate, so we use this
surrogate.
"""
output = super(FeasibilityIterateMonitor, self).postprocess(alg, output)
if self.diagnostic:
stats = output['stats']
stats['gaps'] = self.gaps[1:alg.iter+1]
stats['shadow_changes'] = self.shadow_changes[1:alg.iter+1]
if self.truth is not None:
stats['rel_errors'] = self.rel_errors[1:alg.iter+1]
return output