Commit ef7e0b6d authored by Russell Luke's avatar Russell Luke
Browse files

quadratic prox in progress

parent 6897e43e
......@@ -19,23 +19,7 @@ class ADM_Context:
prox operators and ADM related classes.
"""
def __init__(self, experiment):
# this is a bit heavy handed, but seems to work:
# self=experiment
""" self.product_space_dimension = experiment.product_space_dimension
# primal_scaling may be set as well by the caller of eval method
self.primal_scaling = experiment.primal_scaling
# image_scaling may be set as well by the caller of eval method
self.image_scaling = experiment.image_scaling
# augmentation_scaling may be set as well by the caller of eval method
self.augmentation_scaling = experiment.augmentation_scaling
# lambda is a relaxation paramter and may be set as well by the caller of eval method
self.lmbda = experiment.lambda_0
# for Bregman prox:
self.dBregman_potential = experiment.dBregman_potential
# why wouldn't the following work:
self.coupling_function = experiment.coupling_function
"""
# the rest below is phase retreival specific
# the rest below is phase retreival specific
if hasattr(experiment, 'illumination'):
self.illumination = experiment.illumination
else:
......@@ -72,7 +56,6 @@ class ADM_Context:
else:
# put the image_dim into the Ny ``slot"
self.Ny = experiment.image_dim
if hasattr(experiment, 'domain_objective'):
self.domain_objective = experiment.domain_objective
else:
......@@ -181,11 +164,11 @@ class Prox_NoLips_primal_l0_plus_Bregman(ProxOperator, ADM_Context):
"""
# This first part completes the action of the operator T_2 in
# the fixed point formulation of NoLips ADM:
dQ_xk_mat = dQuadratic.eval(self, u[0])
Q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = Q_xk_vec - u[1]
dq_xk_mat = dQuadratic.eval(self, u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = q_xk_vec - u[1]
# augmentation scaling is the parameter rho in the paper
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dQ_xk_mat), image_resid_vec)
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dq_xk_mat), image_resid_vec)
if self.dBregman_potential==None:
dphi_xk_vec = 0
else:
......@@ -288,11 +271,11 @@ class Prox_NoLips_primal_l1_plus_Bregman(ProxOperator, ADM_Context):
"""
# This first part completes the action of the operator T_2 in
# the fixed point formulation of NoLips ADM:
dQ_xk_mat = dQuadratic.eval(self, u[0])
Q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = Q_xk_vec - u[1]
dq_xk_mat = dQuadratic.eval(self, u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = q_xk_vec - u[1]
# augmentation scaling is the parameter rho in the paper
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dQ_xk_mat), image_resid_vec)
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dq_xk_mat), image_resid_vec)
if self.dBregman_potential==None:
dphi_xk_vec = 0
else:
......@@ -394,11 +377,11 @@ class Prox_NoLips_primal_l2_plus_Bregman(ProxOperator, ADM_Context):
"""
# This first part completes the action of the operator T_2 in
# the fixed point formulation of NoLips ADM:
dQ_xk_mat = dQuadratic.eval(self, u[0])
Q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = Q_xk_vec - u[1]
dq_xk_mat = dQuadratic.eval(self, u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = q_xk_vec - u[1]
# augmentation scaling is the parameter rho in the paper
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dQ_xk_mat), image_resid_vec)
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dq_xk_mat), image_resid_vec)
if self.dBregman_potential==None:
dphi_xk_vec = 0
else:
......@@ -513,11 +496,11 @@ class Prox_NoLips_primal_0_plus_Bregman(ProxOperator, ADM_Context):
"""
# This first part completes the action of the operator T_2 in
# the fixed point formulation of NoLips ADM:
dQ_xk_mat = dQuadratic.eval(self, u[0])
Q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = Q_xk_vec - u[1]
dq_xk_mat = dQuadratic.eval(self, u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = q_xk_vec - u[1]
# augmentation scaling is the parameter rho in the paper
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dQ_xk_mat), image_resid_vec)
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dq_xk_mat), image_resid_vec)
if self.dBregman_potential==None:
dphi_xk_vec = 0
else:
......@@ -569,6 +552,232 @@ class Prox_NoLips_primal_0_plus_Bregman(ProxOperator, ADM_Context):
u_vec = tmp * u_prime
return u_vec
class Prox_NoLips_primal_0_plus_3prox(ProxOperator, ADM_Context):
"""
Prox operator for the NoLips experiment with zero primal objective
and a proximal term to the power 3.
In this experiment we are trying quadratic
approximation of the coupling function. The primal prox is approximated
numerically with a call to the samsara nonlinear optimization routine.
"""
def __init__(self, experiment):
# because of multiple inheritance we call explicitely
# the __init__ method of the two parent classes
ProxOperator.__init__(self, experiment)
ADM_Context.__init__(self, experiment)
if hasattr(experiment, 'A_tens'):
self.A_tens = experiment.A_tens
else:
self.A_tens = 0
if hasattr(experiment, 'b_mat'):
self.b_mat = experiment.b_mat
else:
self.b_mat = 0
if hasattr(experiment, 'c_vec'):
self.c_vec = experiment.c_vec
else:
self.c_vec = 0
if hasattr(experiment, 'Lip_grad_max'):
self.Lip_grad_max = experiment.Lip_grad_max
else:
self.Lip_grad_max = 1.0
if hasattr(experiment, 'Lip_grad_k'):
self.Lip_grad_k = experiment.Lip_grad_k
else:
self.Lip_grad_k = 1.0
if hasattr(experiment, 'augmentation_scaling'):
self.augmentation_scaling = experiment.augmentation_scaling
else:
self.augmentation_scaling = 1.0
if hasattr(experiment, 'Gprox'):
self.Gprox = experiment.Gprox
else:
self.Gprox = None
if hasattr(experiment, 'dGprox'):
self.dGprox = experiment.dGprox
else:
self.dGprox = None
if hasattr(experiment, 'dBregman_potential'):
self.dBregman_potential = experiment.dBregman_potential
else:
self.dBregman_potential = None
if hasattr(experiment, 'coupling_approximation'):
self.coupling_approximation = experiment.coupling_approximation
else:
self.coupling_approximation = 'linear'
if hasattr(experiment, 'primal_scaling'):
self.primal_scaling = experiment.primal_scaling
else:
self.primal_scaling = 1.0
if hasattr(experiment, 'image_scaling'):
self.image_scaling = experiment.image_scaling
else:
self.image_scaling = 1.0
def eval(self, u):
"""
Computes the generalized prox of the 0-function via numerical approximation
Had to write this in an odd way. Formally, the ADM method
for the NoLips problem structure
works on 2 blocks of variables, (x,y) with different
dimensions. In the fixed point interpretation of this
the ADM algorithm is written as the composition of two operators
T_1T_2, each self-mappings on the domain of x. For monitoring
purposes, I keep the (x,y) block structure, but then have to split
the T_2 operator between this primal operator and the operator in the
image domain. If that makes any sense.
Input
----------
u : cell of 2 vectors, the primal domain variable of dim n
and the auxilliary image variable of dim m
prox_idx : int, optional
Index of this prox operator
scaling : scaling in the l0-prox
Returns
-------
u_vec : vector of dim n
"""
# This first part computes to quadratic approximation of the
# primal part of the coupling:
dq_xk_mat = dQuadratic.eval(self, u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
image_resid_vec = q_xk_vec - u[1]
# augmentation scaling is the parameter rho in the paper
dQxkyk_vec = self.augmentation_scaling * matmul(transpose(dq_xk_mat), image_resid_vec)
zk_vec = dQxkyk_vec
# Initial guess: the linear prox point
# self.lmbda is alpha/Lbar in the paper
uold_vec = u[0] - zk_vec
unew_vec = None
stepsize = 999.
prox_objective_new = None
dprox_objective_new_vec = None
# This second part is just the action of the operator T_1 in
# the fixed point formulation of NoLips ADM:
# self.lambda is alpha/Lbar in the paper
# end initial guess
if self.coupling_approximation=='quadratic':
HQ_xk_tens = HQuadratic.eval(self, u[0])
HQ_xk_mat = zeros(shape(HQ_xk_tens[0,:,:]))
for j in range(len(image_resid_vec)):
HQ_xk_mat += image_resid_vec[j]*HQ_xk_tens[j,:,:]
HQ_xk_mat += matmul(transpose(dQ_xk_mat,[1,0]), dQ_xk_mat)
HQ_xk_mat *= self.augmentation_scaling
else:
HQ_xk_mat = zeros(shape(u[0]))
# the objective for the generalized prox is:
prox_objective0_vec = dQxkyk_vec + matmul(HQ_xk_mat, u[0])
prox_objective1_vec = matmul(HQ_xk_mat, uold_vec)
prox_objective_old = matmul(prox_objective0_vec+prox_objective1_vec, uold_vec)
primal_resid_vec = uold_vec - u[0]
norm_primal_resid = matmul(primal_resid_vec, primal_resid_vec)
prox_objective_old += ((1/3)/self.lambd)*norm_primal_resid**3
dprox_objective_old_vec = prox_objective0_vec + prox_objective1_vec
dprox_objective_old_vec += 1/self.lambd*norm_primal_resid*primal_resid_vec
norm_dprox_objective_new_vec = 999.
# Samsara options
options = {}
options['verbose'] = True
options['alpha_0'] = 5e-12
options['gamma'] = .5
options['c'] = .01
options['beta'] = .9999999
options['eta1'] = .995
options['eta2'] = .8
options['eta3'] = .25
# options['eta4'] = .025
options['maxmem'] = 8
# options['lookback'] = 0
options['tr'] = 1e+15
options['ngrad_TOL'] = 1e-6
options['step_TOL'] = 1e-12
options['Maxit'] = 1000
# options['QN_method'] = Broyden1_product
# options['Step_finder'] = Explicit_TR
# options['History'] = cdSY_mat
# options['update_conditions'] = 'Exact'
options['initial_step'] = 1e+5*options['step_TOL']
options['QN_method'] = BFGS1_product
# options['Step_finder'] = Explicit_TR
options['History'] = cdSY_mat
options['update_conditions'] = 'Trust Region'
options['initial_step'] = 1e-7
s = Samsara(**options)
# Tolerances
ngrad_TOL = 2e-14
step_TOL = 0.0001*ngrad_TOL
Maxit = 1500
it = 0
# enter Samsara main loop
while it < Maxit and norm_dprox_objective_new_vec > ngrad_TOL and stepsize > step_TOL:
unew_vec, uold_vec, prox_objective_old, dprox_objective_old_vec, stepsize =\
s.run(uold_vec, unew_vec, prox_objective_old, prox_objective_new,
dprox_objective_old_vec, dprox_objective_new_vec)
it += 1
# evaluate new objective value and new gradient:
primal_resid_vec = unew_vec - u[0]
norm_primal_resid = matmul(primal_resid_vec, primal_resid_vec)
prox_objective_new += ((1/3)/self.lambd)*norm_primal_resid**3
dprox_objective_new_vec = prox_objective0_vec + prox_objective1_vec
dprox_objective_new_vec += 1/self.lambd*norm_primal_resid*primal_resid_vec
norm_dprox_objective_new_vec = matmul(dprox_objective_new_vec, dprox_objective_new_vec)
if stepsize <= step_TOL:
print('Quasi_Newton Algorithm stagnated: stepsize tolerance violated.')
if it >= Maxit:
print('Quasi-Newton Maxit exceeded: maximum step count violated.')
u_vec = unew_vec
return u_vec
# The following only holds if (F+phi)=0. If
# we are in this class, F=0, so the only thing to check is
# that phi=0
if (dphi_xk_vec == 0):
u_vec = call to samsara
else:
print('Message from ADM_prox.py, class Prox_NoLips_primal_0_plus_Bregman:')
print('The only implementation so far is for no Bregman prox.')
u_vec = None
# This should throw an error
else:
# the linear approximation is the default
# self.lmbda is alpha/Lbar in the paper
u_prime = u[0] - self.lmbda * (zk_vec)
# This second part is just the action of the operator T_1 in
# the fixed point formulation of NoLips ADM:
# self.lambda is alpha/Lbar in the paper
prox_parameter = self.lmbda / self.primal_scaling
if self.dBregman_potential==None:
tmp=1
else:
# using map avoids for... loops at each of the inner functions!
a = sum(u_prime**2) * prox_parameter * self.Lip_grad_max
# tmp = real root of a*t^3 + t - 1 via Cardano's formula
if abs(a)>=1e-5:
descr = sqrt(0.25/(a**2) + (1/27)*(1/(a**3)))
tmp = cbrt(0.5/a + descr) + cbrt(0.5/a - descr)
else:
tmp=1
u_vec = tmp * u_prime
return u_vec
def prox_objective(self, u):
class Prox_NoLip_image_l1_at_q(ProxOperator, ADM_Context):
"""
Prox operator for the NoLips experiment with l1-penalized quadratic constriants
......@@ -630,9 +839,9 @@ class Prox_NoLip_image_l1_at_q(ProxOperator, ADM_Context):
-------
y_new : vector of dim u[1]
"""
Q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
q_xk_vec = Quadratic(self.A_tens, self.b_mat, self.c_vec,u[0])
shrinkage_parameter = self.augmentation_scaling / self.image_scaling
y_new = proxoperators.Prox_l1.eval(Q_xk_vec, shrinkage_parameter)
y_new = proxoperators.Prox_l1.eval(q_xk_vec, shrinkage_parameter)
return y_new
class Prox_NoLip_image_linfty_at_q(ProxOperator, ADM_Context):
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment