unconstrainded optimization aproach to singular poisson problem with tao petsc solvers

Question

I am trying to define a funcional obj = gamma * norm(u, 2) + norm(a - L, 2)
in order to use some PETSc TAO solver like tron or cg to an anternative
approach to Singular Poisson problem (singular_poisson.py).

How can I define this functional in a correc way?

from dolfin import *

parameters["linear_algebra_backend"] = "PETSc"

# Create mesh
k = 16
mesh = UnitSquareMesh(k, k)

# Create function space
V = FunctionSpace(mesh, "Lagrange", 1)

# Create solution, trial and test functions
u, du, v = Function(V), TrialFunction(V), TestFunction(V)

# Variational formulation
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")

a = inner(grad(u), grad(u))*dx
L = f * u * dx + g * u * ds

gamma = 0.1

obj = gamma * norm(u, 2) + norm(a - L, 2)

grad_obj = derivative(obj, u, v)
hess_obj = derivative(grad_obj, u, du)

# Box constrains
constraint_l = Expression("x_min", x_min=-1.2)
constraint_u = Expression("x_max", x_max=1.2)

u_min = interpolate(constraint_l, V)
u_max = interpolate(constraint_u, V)

# Define the minimisation problem by using OptimisationProblem class
class Problem(OptimisationProblem):

    def __init__(self):
        OptimisationProblem.__init__(self)

    # Objective function
    def f(self, x):
        u.vector()[:] = x
        return assemble(obj)

    # Gradient of the objective function
    def F(self, b, x):
        u.vector()[:] = x
        assemble(grad_obj, tensor=b)

    # Hessian of the objective function
    def J(self, A, x):
        u.vector()[:] = x
        assemble(hess_obj, tensor=A)

# Create the PETScTAOSolver
solver = PETScTAOSolver()

# Set some parameters
solver.parameters["method"] = "tron"
solver.parameters["monitor_convergence"] = True
solver.parameters["report"] = True
solver.parameters["maximum_iterations"] = 1000
#info(parameters, True)    # Uncomment this line to see the available parameters

# Parse (PETSc) parameters
parameters.parse()

# Solve the problem
set_log_level(DEBUG)
solver.solve(Problem(), u.vector(), u_min.vector(), u_max.vector())

# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u

# Plot solution
plot(u, interactive=True)

I get the error message:

---

RuntimeError Traceback (most recent call last)
/Users/andremachado/Desktop/demo_singular-poisson_v4_TAO.py in ()
23 gamma = 0.1
24
---> 25 obj = gamma * norm(u, 2) + norm(a - L, 2)
26
27 grad_obj = derivative(obj, u, v)

/Users/andremachado/.hashdist/bld/profile/q3rqj2knkdy7/lib/python2.7/site-packages/dolfin/fem/norms.py in norm(v, norm_type, mesh)
100 cpp.dolfin_error("norms.py",
101 "compute norm",
--> 102 "Norm type must be a string, not " + str(type(norm_type)))
103 if mesh is not None and not isinstance(mesh, Mesh):
104 cpp.dolfin_error("norms.py",

/Users/andremachado/.hashdist/bld/profile/q3rqj2knkdy7/lib/python2.7/site-packages/dolfin/cpp/common.py in dolfin_error(location, task, reason)
2483
2484 """
-> 2485 return _common.dolfin_error(location, task, reason)
2486
2487 def deprecation(feature, version_deprecated, version_remove, message):

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** fenics@fenicsproject.org

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to compute norm.
*** Reason: Norm type must be a string, not .
*** Where: This error was encountered inside norms.py.
*** Process: unknown

---

*** DOLFIN version: 1.5.0
*** Git changeset:
*** -------------------------------------------------------------------------

What is the correct way to set a funcional
like obj = gamma * norm(u, 2) + norm(a - L, 2),
gamma = 0.01 with ufl in order to use the
powerfull (box constrain) PETSc TAO optimization
solver?

Answer 1

The solution to the Poisson problem is a minimizer of the functional

0.5 * grad(u)**2 * dx - L,

so I think your objective should be defined as

obj = gamma * u**2 * dx + .5 * grad(u)**2 * dx - L

You should also be aware the "singular Poisson" problem your are solving only has solutions if the right hand side L annihilates the constants. This seems not to be the case here.

Comments

Magne, take a look at it.

I am using a optimization solver less powerfull than TAO
one - CVXPY and I get the same solution minimizing
gamma * norm(u, 2) + norm(Ax -b, 2) with bound constrains
x >= -0.5, A, b = assembler_system(a, L) via

gamma = 1.0e-4

obj = cvxpy.Minimize(cvxpy.norm2(gamma * x) + cvxpy.norm2(A.array()*x-b.array()))

cvxpy.constrains = [ x >= -0.5 ]


import cvxpy

from dolfin import *

k = 16
mesh = UnitSquareMesh(k, k)

V = FunctionSpace(mesh, "Lagrange", 1)

u = TrialFunction(V)
v = TestFunction(V)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")

a = inner(grad(u), grad(v))*dx
L = f * v * dx + g * v * ds

A, b = assemble_system(a, L)

u = Function(V)

#CVXPY solver
n_dic = (A.array()[0, :]).size

x = cvxpy.Variable(n_dic)

gamma = 1.0e-4
obj = cvxpy.Minimize(cvxpy.norm2(gamma * x) + cvxpy.norm2(A.array()*x-b.array()))
cvxpy.constrains = [ x >= -0.5 ]

prob = cvxpy.Problem(obj, cvxpy.constrains)
prob.solve(solver='ECOS', verbose=True, max_iters=250) #SCS, ECOS_BB, CVXOPT

x_value = np.array(x.value)
u.vector()[:] = x_value[:, 0]

plot(u, interactive=True)

I have tried the functional obj = gamma * u2 * dx + 0.5 * grad(u)2 * dx - L with box constrains -0,5 <= x <= 10.0 and all TAO solvers tron, blmvm, nls, nm, ntr, cg but the solution is different.

( This 0.5 * grad(u)**2 * dx - L is norm(a - L, 2) for poisson equation? )

It doesn't make sense to me...

---

First of all you need to fix your problem formulation since you currently try to solve a problem that doesn't have solutions. You are using the same expressions for f and g as in singular_poisson.py, but in the demo the assembled right hand side is changed to be consistent before the linear solve.

A minimizer of the functional J defined

L = f * u * dx + g * u * ds
J = 0.5 * grad(u)**2 * dx - L

is a solution to the Poisson problem. Conversely, a solution of the Poisson equation is a minimizer of J. Solutions exists if and only if the linear functional L vanishes on the constant functions. When solutions exists, they are not unique.

Now, assume that L permits solutions. If you define your objective functional as

obj = .5 * gamma * u**2 * dx + J

then it can be cosidered a regularization of the functional J. For small gamma, a minimizer of obj may be close to a minimizer of J, hence close to the solution of the Poisson problem.
Note that this functional has unique minimizer, regardless of L.

Finally, consider that you can define an iterate procedure to eliminate the dependence on gamma, for example by setting

obj = .5 * gamma * (u- u_prev)**2 * dx + J

This approach is equivalent solving the heat equation with implicit Euler and time step 1/gamma, until it converges to a stationary solution.