import numpy as np
from numpy import zeros, isreal, reshape, real, imag
from proxtoolbox.utils.size import size_matlab
from proxtoolbox.utils.cell import Cell, isCell
import importlib
[docs]class QNAccelerator:
"""
Quasi-Newton acceleration of fixed point iterations
based on Samsara, a reverse communication nonlinear
optimization package. Samsara must be installed in
the folder proxtoolbox/algorithms/samsara.
Based on Matlab code written by Russell Luke (Inst. Fuer
Numerische und Angewandte Mathematik, Universitaet
Gottingen) on Feb 23, 2018.
"""
def __init__(self, experiment):
self.samsara_data = []
samsara_module_name = 'proxtoolbox.algorithms.samsara'
try:
module = importlib.import_module(samsara_module_name)
samsara_attr = getattr(module, 'Samsara')
self.samsara = samsara_attr() # instantiate Samsara
except Exception as e:
errMsg = e.__str__()
print('*********************************************************************')
print('* SAMSARA package missing. Please download from *')
print('* http://vaopt.math.uni-goettingen.de/en/software.php *')
print('* Save and unpack the datafile then copy the python/samsara folder *')
print('* in the proxtoolbox/algorithms folder. *')
print("* You may need to execute setup.py install again *")
print('******************************************************************')
raise
self.unew_vec = 0
self.gradfnew_vec = 999
[docs] def evaluate(self, u, alg):
"""
Quasi-Newton acceleration evaluation
Parameters
----------
u : ndarray or a list of ndarray objects
The point to be accelerated.
alg : Algorithm object
The algorithm that uses this accelerator
Returns
-------
u_new : ndarray or a list of ndarray objects
(same type as input parameter u)
The step generated by the quasi-Newton routine (same type as input parameter `u`).
"""
# Warning: only the complex case when u is not a cell was tested
# The original Matlab code, when u is a list, may not have been used
# in demos and may not work.
# Reshape u into a vector. Use order 'F' to match Matlab column-major order
# Warning: if the user sends in u which is
# already a point on the diagonal of the product space, this vector will be
# larger than necessary.
if isCell(u):
# determine data type and store shapes and sizes of each component in u
# assume data is either real or complex
dim_u = len(u)
u_dtype = np.dtype('float64') # default data type
u_shapes = Cell(dim_u)
u_sizes = Cell(dim_u)
isReal = True
for j in range(dim_u):
if np.iscomplexobj(u[j]):
u_dtype = u[j].dtype
isReal = False
u_shapes[j] = u[j].shape
u_sizes[j] = u[j].size
# flatten data - follows Matlab column-major order
unew_vec = None
for j in range(dim_u):
u_j_flat = u[j].flatten('F')
if isReal:
if unew_vec is None:
unew_vec = u_j_flat
else:
unew_vec = np.concatenate((unew_vec, u_j_flat))
else: # complex case
u_j_flat_real = np.concatenate((real(u_j_flat), imag(u_j_flat)), axis=None)
if unew_vec is None:
unew_vec = u_j_flat_real
else:
unew_vec = np.concatenate((unew_vec, u_j_flat_real))
else:
u_shape = u.shape
u_size = u.size
u_flat = u.flatten('F')
m, n, p, q = size_matlab(u)
isReal = True
if np.isrealobj(u):
unew_vec = u_flat
else: # complex case
isReal = False
unew_vec = np.concatenate((real(u_flat), imag(u_flat)), axis=None)
uold_vec = self.unew_vec # we do not need to do a copy here
gradfold_vec = self.gradfnew_vec # we do need to do a copy here
gradfnew_vec = uold_vec - unew_vec
# make call to SAMSARA for acceleration step. method_input.stats.gap
# is the objective function (surrogate). Requires running in diagnostic
# mode since this is not always generated by the algorithm. If it
# is absent, then samsara makes the norm of the residual the surrogate
# objective function
if alg.iter > 3:
if hasattr(alg.iterateMonitor, 'gaps'):
if isCell(u):
dim_prod = u[0].size
else:
dim_prod = m*n*p*q # could also use u.size
gaps = alg.iterateMonitor.gaps
unew_vec, dummy1, dummy2, dummy3, dummy4 = \
self.samsara.run(uold_vec, unew_vec,
gaps[alg.iter-2]*dim_prod,
gaps[alg.iter-1]*dim_prod,
gradfold_vec, gradfnew_vec)
else:
unew_vec, dummy1, dummy2, dummy3, dummy4 = \
self.samsara.run(uold_vec, unew_vec,
[],
[],
gradfold_vec, gradfnew_vec)
# reshape and store information
if isCell(u):
u_new = Cell(len(u))
start = 0
end = 0
for j in range(len(u)):
size_j = u_sizes[j]
if isReal:
end += size_j
u_new[j] = reshape(unew_vec[start:end], u_shapes[j], order='F')
else: # complex case
end += 2*size_j
u_new_real = unew_vec[start:start+size_j]
u_new_im = unew_vec[start+size_j:end]
u_new[j] = reshape(u_new_real + 1j*u_new_im, u_shapes[j], order='F')
start = end
else:
if isReal:
u_new = reshape(unew_vec, u_shape, order='F')
else: # complex case
u_new_real = unew_vec[0:u_size]
u_new_im = unew_vec[u_size:2*u_size]
u_new = reshape(u_new_real + 1j*u_new_im, u_shape, order='F')
self.unew_vec = unew_vec
self.gradfnew_vec = gradfnew_vec
return u_new