Backtracking-Armijo Line Search Algorithm

Question

EDIT based on comments below: I add the mathematical formulation of my problem below. I am trying to solve an equation of the form $$ \partial_t f(x,y,t)= (\partial^2_x +\partial^2_y) f(x,y,t) \equiv G(x,y,t), $$ discretizing this equation we have $$ f^{k+1}_{i,j}= f^k_{i,j} + \Delta t G^k_{i,j} $$ where $i,j$ refer to discretized spatial coordinates $x,y$ and $k$ corresponds to the iteration step.
However here $\Delta t $ is a fixed step size, I want to use a line search to find an optimal step size. Defining $\Delta t \equiv \alpha_k$, I want to find $\alpha_k $ such that $f^{k+1}(i,j) < f^{k}(i,j)-c\alpha_k G^\top G$ which is a backtracking Armijo line search. So the equation I am trying to solve is : $$ f^{k+1}_{i,j}= f^k_{i,j}+\alpha^k G^k (i,j) $$ Below is a back tracking line search algorithm to find $\alpha_k$ but it is not being computed correctly I realize.

I updated my algorithm based on the comment below however it still seems my stepsize at each iteration, $\alpha_k$ is not being updated properly. When I print it out it just prints the initial value I inputted for it. Is this algorithm not updating $\alpha_k$ correctly? I thought the point of the backtracking line search was to find me an optimal value $\alpha_k$ such that I get to the minimum. How can I fix this? thanks!

I am trying to code the backtracking-Armijo line search algorithm on page 10 here https://people.maths.ox.ac.uk/hauser/hauser_lecture2.pdf.

Below is a sample code for a back tracking line search algorithm . I can check that the algorithm is not correct but I am not sure where I am going wrong. A few errors i've realized are possibly the if statement condition and updating alphak at the very end ; lastly I'm not sure what else is wrong (I am unsure how to fix these problems). I tried to follow the algorithm in the book but it is not too detailed. It is very clear there is a problem with the algorithm but I am not sure what.

Note, the search direction I choose given by Pk is in the direction of the negative gradient, which I call -g. Assume below that g(i,j) and fk(i,j) are given at the first iteration, and are 2D arrays since they depend on spatial positions i,j.

integer, parameter :: nx=10,ny=10, k=10
real, dimension(-nx:nx,-ny:ny) :: fk,fk1,g,gt,Pk
integer :: i,j,m
real :: alphak,c,rho !step size at iteration k

c=0.0001
rho=0.5

do m=1,k


        alphak = 1.0
        Pk(i,j) = -g(i,j) !search direction = -gradient
        fk1(i,j) = fk(i,j) + alphak*Pk(i,j)
        gt(i,j) = g(j,i) !transpose of g

           if (fk1(i,j) > fk(i,j)-c*alphak*gt*g) then

                alphak = rho*alphak

                   do j=-ny+1,ny-1
                   do i=-nx+1,nx-1
                      fk1(i,j) = fk1(i,j)-alphak*g(i,j)
                   end do
                   end do

           end if

    print*, "print out alphak for m=', m
    print "(//(5(5x,e22.14)))", alphak

end do

Answer 1

[The following references a paper by Ascher et al. ...]

You want to solve

$$ \partial_t f(x,y,t)= (\partial^2_x +\partial^2_y) f(x,y,t), $$

subject to some initial and boundary conditions. However, as I understood from your comments, you are in fact interested in the steady state solution (and not in the transient behavior of $f$). This steady state is of course the $f$ that satisfies

$$ \frac{\partial f}{\partial t} = 0. $$

You start by choosing a semi-discrete form of your PDE, namely: $$ \partial_t f^h(t) = - \left( A^h f^h - b^h \right), $$

where, if your discrete computational grid is $N \times N$, then $f^h \in R^{N^2}$ is your state vector (one DOF per grid point but reshaped into a column vector), $A^h$ is the discrete Laplacian operator (symmetric and positive definite), and $b^h$ includes forcing terms and boundary conditions (if any). See equations (2a) - (2b) in the paper I linked to for an example finite volume discretization.

Now, an obvious second step would be to use Euler's (explicit) method to march in time towards the steady state, i.e. go from $t^n$ to $t^{n+1}$ using

$$ f^{n+1,h} = f^{n,h} + \alpha \left( b^h - A^h f^{n,h} \right). $$

Essentially what you're doing here is to take a step $\alpha$ along the direction of the residual

$$r^{n,h} := b^h - A^h f^{n,h}. $$

Trouble is, a "default" Euler time step $\alpha$ will be small (it is constrained by numerical stability) so converging to the steady state will take a while with this approach. That's why you want to use a higher time step if possible, and this is where the line search algorithm comes into play. You will essentially try to take the best step along the residual direction at each iteration.

To do that, observe that the Euler discretization above can be interpreted as a gradient descent method for a convex numerical optimization problem:

$$ \min_{f^h} J(f^h) := \frac{1}{2} \left( (f^h)^T A^h (f^h) - (b^h)^T f^h \right). $$

Here the cost function $J$ is a scalar. The necessary condition for a minimum is $\nabla_f{^h} J = 0$ i.e.

$$A^h f^{h} = b^h \leftrightarrow r^h = 0$$

Hence the necessary (and sufficient, since the cost function is convex) conditions give you the steady state of the discrete PDE.

Returning to the line search, a natural choice for the $n$-th time step is that who minimizes the cost function $J$ along the residual line, i.e., the solution to

$$ \min_\alpha J(f^{n,h} + \alpha r^{n,h}) $$

This gives you the steepest descent step

$$ \alpha = \frac{(r^{n,h})^T(r^{n,h})}{(r^{n,h})^T A^h (r^{n,h})}. $$

This is what's called an exact line search. However, minimizing $J$ may not be cost effective for more complicated cost functions. Instead, people have come up with Armijo-type backtracking searches that do not look for the exact minimizer of $J$ along the search direction, but only require sufficient decrease in $J$: you iterate over $\alpha$ until

$$ J(f^{n,h} + \alpha r^{n,h}) < J(f^{n,h}) + \alpha \beta (r^{n,h})^T r^{n,h} $$

for some fixed $\beta \in (0,1)$.

Other step choices are of course possible (see Ascher's paper). If you're looking for an implementation of steepest descent with Armijo, check out this MATLAB code.

Pseudocode:

1. $n \leftarrow 0$, $t^n \leftarrow t^0, f^h = f^{0,h} $ the initial conditions, $f$ is a 1D array
2. $r^h = b^h - A^h f^h$ (the residual, also a 1D array)

3. while $r^h$ not small enough

   (a) run the Armijo iteration using $J$ (as shown in the slides you've referenced); this gives you the time step $\alpha$

   (b) update the PDE state using the Euler formula $f^{n+1,h} \leftarrow f^{n,h} + \alpha r^{n,h}$

   (c) update the residual $r^{n+1,h} = b^h - A^h f^{n+1,h}$

   (d) check if the new residual is small enough i.e. $\| r^{n+1,h} \| < \epsilon$ (for small $\epsilon > 0$, the exact choice of norm and $\epsilon$ is up to you); if so break

   (e) $t^{n+1} \leftarrow t^n + \alpha$; $t^n \leftarrow t^{n+1}$; $n \leftarrow n+1$

   end while

I'd suggest you try implementing this in the programming language you're most familiar with, then port that working implementation to Fortran.

An alternative technique is pseudo-transient continuation.