I want to implement a Gradient Descent Algorithm in FEniCS. I am working in 1D (if this works out I can move onto 2D).
My expression for the energy is ( I choose a simple one for now) E[u]=\int_0^{L} (1-u)^2 dx
I present the MWE below, where i am trying to do the gradient descent step, ie, u^{n+1}=u^n-\gamma \partial_u E[u],
where \gamma is the learning rate.
The MWE is presented below:
import dolfin
from dolfin import *
import ufl
import matplotlib.pyplot as plt
from petsc4py import PETSc
L=10
mesh = IntervalMesh(100,0,L)
element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, "CG", 1)
gamma = float(0.1)# Learning rate.
u = Function(V)
u_up = Function(V)
phi = (1-u)**2
# Total potential energy
Pi = phi*dx
#Defining derivative
F = ufl.derivative(Pi,u)
# Set up the non-linear problem
u_up = u - gamma*F
plot(u)
plt.title(r"$u(x)$",fontsize=26)
plt.show()
plot(u_up)
plt.title(r"$u_up(x)$",fontsize=26)
plt.show()
The error is
Traceback (most recent call last):
File "/Users/DigoSen/Desktop/Books/Codes/Python/FEniCS/MWE/Discourse-Post4.py", line 25, in <module>
u_up = u - gamma*F
~~^~~~~~~~~
File "/Users/DigoSen/anaconda3/envs/fenics2019/lib/python3.12/site-packages/ufl/form.py", line 304, in __rsub__
return other + (-self)
~~~~~~^~~~~~~~~
TypeError: unsupported operand type(s) for +: 'Function' and 'Form'
I understand that F and u are not the same type. But I don’t know how to convert them to the same type.
I am using dolfin version 2019.1.0, which i installed using conda-forge.
Hi, I don’t quite follow the series of commands that you mentioned.
When you say assemble Pi, I’m guessing its something like Pi = assemble(phi*dx)
“Then you should choose an appropriate Riesz representation of the gradient to put it in V.”
Im not sure i understand your suggestion. Are you saying I should supply the variational derivative? I was hoping I could get the code to compute it itself. Please correct me if i misunderstood your suggestion.
Well, I wrote this a bit quickly on my phone. What you should do is first:
Convert F from a symbolic representation into a numerical representation.
This can be done with F_vector = assemble(F)
However, F is in the dual space of u, and you can either use the values of the “vector” directly, or use a Riesz representation. For a Riesz representation choose a coersive bilinear form such that a(u, v) == F and solve that variational problem.
If you do not want to do that, you can simply do
from dolfin import *
L = 10
mesh = IntervalMesh(100, 0, L)
element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, "CG", 1)
gamma = float(0.1) # Learning rate.
u = Function(V)
u_up = Function(V)
phi = (1-u)**2
# Total potential energy
Pi = phi*dx
# Defining derivative
F = derivative(Pi, u)
F_vector = assemble(F)
# Set up the non-linear problem
u_up.vector()[:] = u.vector()[:] - gamma*F_vector[:]
This is a continuation of the question. Should I add this as a separate post?
I want to implement the Nesterov momentum method. Where instead of computing the gradient at the current point, the gradient at a different point is used.
To talk briefly about the method u_{i+1}=u_i+v_{i+1}, v_{i+1}=\tau v_i -\eta\nabla f(u_i+\tau v_i), \text{with}~v_0=0,~~\tau\in(0,1).
I was wondering what changes i would have to do to the line below.
# Defining derivative
F = derivative(Pi, u)
F_vector = assemble(F)
In the second line above, I end up with F_vector being evaluated as the gradient at the current iteration step. How do i enforce that it be evaluated at u_i+\tau v_i?
You should be replacing the values currently stored in u with those of the updated solution u + tau v. For instance, using assign(u, function_with_u_plus_tau_u) or, as you do in your previous code u.vector()[:] = u.vector()[:] + tau*v.vector()[:]