Add preconditioning in GD #1359

Currently we have the following for the GD algorithm:

$$ x_{k+1} = x_{k} - \gamma_{k} \nabla f (x_{k}) $$

However, we need to have a preconditioning before the f.gradient step as discussed in #1345. This is both useful in the deterministic and stochastic step.

$$ x_{k+1} = x_{k} - \gamma_{k} D(x_{k})\nabla f (x_{k}) $$

Notice that we already have this functionality in the SIRT algorithm but it is hardcoded. In fact in SIRT, we have two fixed preconds, which are

$$ \frac{1}{A^{T}\mathbb{1}} $$

and

$$ \frac{1}{A \mathbb{1}} $$

Of course SIRT is designed for a specific optimisation problem where we know that the objective function is

$$ ||Ax - d||^{2} $$

and our users will only need to pass the operator A and the data d. However, if GD , ISTA and FISTA allow preconditioning we can have an algorithm similar to SIRT

• without non-negativity constraint = GD
• with non-negativity constraint = ISTA/FISTA

And also, we can have their stochastic versions, i.e., OS-SIRT etc. Actually, since SIRT does not accept an objective function, we cannot setup OS-SIRT with the SIRT class via SubsetSumFunction.

As described in #1345, we can have the following that can be used in both deterministic and stochastic cases:

a) preconditioner is fixed $\frac{1}{A^{T}\mathbb{1}}$
b) depends on the iterates $\frac{x_{k}+\delta}{A^{T}\mathbb{1}}$
c) depends on the subsets_num, i.e., activate/deactivate it after some epochs or iterations.
d) None

To cover all of the above cases, I suggest to have sth like

GD-update

    def update(self):
        '''Single iteration'''
        
        self.objective_function.gradient(self.x, out=self.x_update)
        if self.precond is True:
           self.x_update *= self.precond( x iterate, num_subset, iteration)  
        if self.update_step_size:
            # the next update and solution are calculated within the armijo_rule
            self.step_size = self.armijo_rule()
        else:
            self.x.sapyb(1.0, self.x_update, -self.step_size, out=self.x)

And the user can predefine a callable function for the precond step:

def fixed_precond(x, num_subset, iteration):
         return 1./A.adjoint(ag.allocate(1.0)) # sometimes we need to check if there are no 0's in the denominator

def fixed_precond(x, num_subset, iteration):
         return (x + delta)/(A.adjoint(ag.allocate(1.0))

or in the stochastic case

def fixed_precond(x, num_subset, iteration):
         if epochs<5:
               return (x + delta)/(A.adjoint(ag.allocate(1.0))
        else:
              self.precond = False # deactivate precond after 5 epochs  

What do you think @paskino @gfardell @jakobsj