0
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ your notation is a bit confusing. i gather $f^k$ is your current guess at an optimum? then why is the if test written in terms of $f^k$? or is it supposed to indicate the cost function that you minimize? oh, and i think you should not include links to not-so-legal copies of textbooks... $\endgroup$ – GoHokies Mar 7 '17 at 6:56
  • $\begingroup$ You are doing a separate line search for each point (i,j) -- this is of course wrong. The (i,j) loop should only be around the lines where you are updating fk1. $\endgroup$ – Christian Clason Mar 7 '17 at 8:28
  • 2
    $\begingroup$ perhaps it'd be best to look at a line search code that works and try to understand that before you code your own implementation. $\endgroup$ – GoHokies Mar 7 '17 at 18:49
  • 2
    $\begingroup$ @Integrals write down the mathematical formulation of your optimization problem in your original post. we can't give much useful advice about the code if we don't know what problem you're trying to solve. and if this is homework, you should label it as such. $\endgroup$ – GoHokies Mar 7 '17 at 19:02
  • $\begingroup$ @GoHokies okay, I will re-formulate my question so it is very clear. I really appreciate your help and discussion. Thanks $\endgroup$ – Jeff Faraci Mar 7 '17 at 19:03
1
$\begingroup$

[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.

$\endgroup$
  • $\begingroup$ What is $J(f^h)$, is $J(f^h)=(f^h)^\top A^h f^h -(b^h)^\top f^h$? (It seems it must be based on the result of $\nabla_f J=0$. I just want to be very clear how you define $J$. Thanks for your time! I also found that matlab code but that is more difficult to understand than just the backtracking algorithm so I figured I'd start with the simpler case. Based on understanding exactly the form of J, you've showed me now the severe lack of understanding I've had. I will work on the code and update the post if necessary. Thanks a lot. $\endgroup$ – Jeff Faraci Mar 8 '17 at 17:13
  • 1
    $\begingroup$ yes, $J$ is defined exactly as you've written. Is my answer clear enough for you to get the ball rolling? $\endgroup$ – GoHokies Mar 8 '17 at 18:23
  • $\begingroup$ We'll see :) I have a few questions, I will also gladly check it as the answer I just want to be a clear on a bit more first. Please correct me where I am wrong; In order to use this backtracking algorithm, I need to store the following in an array at time iteration k depending on space $i,j$, $$ J(i,j)=f^\top(i,j) A(i,j) f(i,j) -b^\top f(i,j) $$ So I need to FIRST have the numerical solution $f(i,j)$ in order to do this procedure. But this requires having used a time step. So do I first guess an initial step and get a solution and iterate till it finds the best $\alpha_k$? $\endgroup$ – Jeff Faraci Mar 9 '17 at 0:38
  • 1
    $\begingroup$ @Integrals i'll try to post some pseudocode later today. $J$ is a scalar (a single real value), so there is no such thing as $J(i,j)$. you can store and update the solution $f$ as a matrix (2D array), but for the purposes of the line search algorithm, you work with $f$ reshaped into a column vector. that is why you end up with a scalar $J$ $\endgroup$ – GoHokies Mar 9 '17 at 7:01
  • 1
    $\begingroup$ yes a vector = a 1D array. if you store your discrete differential operator $A^h$ as a sparse matrix, and apply it to $r^h$ through a matrix-vector operation, you can avoid using 2D arrays altogether - just keep all variables in vector (1D array shape) to avoid needing to re-shape things at every time step. $\endgroup$ – GoHokies Mar 9 '17 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.