nonlinear minimization

[gdmcbain]

Some finite element problems can be cast as minimizations rather than boundary value problems; can they be cast as such in scikit-fem and be solved with scipy.optimize.minimize?

For example, the solution x of the Young–Laplace equation in ex10 is a minimal surface, constrained by Dirichlet boundary conditions; can it be solved by minimizing a functional for the area?

[gdmcbain]

Yes. See #643.

scikit-fem/docs/examples/ex39.py

Lines 16 to 43 in 52c01bf

	@Functional
	def area(w):
	return np.sqrt(1 + dot(w["w"].grad, w["w"].grad))
	def fun(z: np.ndarray) -> float:
	return area.assemble(basis, w=basis.interpolate(z))
	def jac(z: np.ndarray) -> np.ndarray:
	return rhs.assemble(basis, w=basis.interpolate(z))
	def hess(z: np.ndarray) -> np.ndarray:
	return jacobian.assemble(basis, w=basis.interpolate(z))
	trace = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis.N))
	constraint = LinearConstraint(trace, *[np.sin(np.pi * basis.mesh.p[0, D])]*2)
	x = minimize(
	fun,
	basis.zeros(),
	jac=jac,
	hess=hess,
	method="trust-constr",
	constraints=constraint,
	).x

[gdmcbain]

This is all a little bit like #119, which attempted to use scipy.optimize.root, except that minimize is much more accepting of the jacobian if provided by the user. (It calls the jacobian the hessian and the residual the jacobian because it's thinking about the functional rather than the solution.)

[gdmcbain]

An awkwardness here is the ndarray that has to be constructed for the scipy.optimize.LinearConstraint.

scikit-fem/docs/examples/ex39.py

Lines 33 to 34 in 52c01bf

	trace = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis.N))
	constraint = LinearConstraint(trace, *[np.sin(np.pi * basis.mesh.p[0, D])]*2)

Initially I had tried something like

trace = LinearOperator(lambda z: z[D], (D.size, basis.N))

but this fails as internally SciPy applies it with np.dot(trace, z) instead of trace.dot(z) or trace @ z. This is disastrous and returns something of completely wrong shape with dtype=object. It might be worth trying to get that fixed upstream if much minimization is going to be done with scikit-fem.

https://github.com/scipy/scipy/blob/04e9644ed52b7df431656d5b3d1358e87d77ec9d/scipy/optimize/_constraints.py#L372-L374

——

Edit:— This is incorrect. This code isn't on the path. See below.

[gdmcbain]

According to whpowell96 Feb 19 '20 at 17:27 (Computational Science Stack Exchange, ‘Which SciPy nonlinear solver when Jacobian is analytically known and sparse?’) :

Nonlinear solvers are almost always a better choice than applying optimization algorithms to the residual

Is that right?

(Other answers and comments on the question recover many of the disappointing findings of #119 about scipy.optimize.root.)

[gdmcbain]

Actually, correcting the error above, the linear constraint passes through:

• class PreparedConstraint
• class LinearVectorFunction

This doesn't work for a LinearOperator either though.

[gdmcbain]

The attempt to use scipy.optimize.minimize(method="trust-constr") #643 was abandoned after it was noticed and realized that the inner linear solves were done using Krylov iteration without effectively preconditioning which therefore blew out in count as the mesh is refined. It works but doesn't scale beyond small problems.

[gdmcbain]

As a further demonstration though, I did extend the example to include a constraint on the volume, thereby generalizing from a minimal surface to a constant mean-curvature surface. Also flattening the ring of ex10 to a planar unit square and taking the required volume as one half, putting

trace = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis.N))
volume = csr_matrix(unit_load.assemble(basis)[None, :])

constraint = [
    LinearConstraint(trace, *[0] * 2),
    LinearConstraint(
        volume,
        *[0.5] * 2,
    ),
]

in the code of

https://github.com/gdmcbain/scikit-fem/blob/182907de4f800a07abe276eb7b8c92da3394caf5/docs/examples/ex39.py#L33

produces

[gdmcbain]

Although scipy.optimize.minimize is a bit disappointing, I'm continuing to look into this as the minimization formulation provides an interesting equivalent alternative for some problems and might be required for others (like obstacle problems? Glowinski 1984, §II.2). Also this is more about demonstrating how to set up the functional, jacobian, and hessian in scikit-fem than the performance of the actual minimizer which presumably might be fixed upstream or swapped out, or not be so relevant for smaller problems.

• Glowinski, R. (1984). Numerical Methods for Nonlinear Variational Problems. New York: Springer

As a next example, I thought to redo ex01, the ‘Poisson equation with unit load’. Since it has the form a (u, v) = b (v) for all v, it is equivalent (Glowinski 1984, App. I, §2.4) to the minimization of J (v) = a(v, v) / 2 − b (v), which has jacobian J' (v) = a (v, ⋅) − b and hessian J'' = a . In scikit-fem, the Dirichlet conditions can again be imposed as a scipy.optimize.LinearConstraint.

https://github.com/gdmcbain/scikit-fem/blob/8a900e21cbcb7d8ea6cc21c63bb6525a9d014b46/docs/examples/ex01_minimize.py#L18-L42

def functional(x: np.ndarray) -> float:
    return x @ A @ x / 2 - b @ x

def jacobian(x: np.ndarray) -> np.ndarray:
    return x @ A - b

def hessian(x: np.ndarray) -> np.ndarray:
    return A

D = m.boundary_nodes()
trace = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis.N))
constraint = LinearConstraint(trace, *[np.zeros(D.size)]*2)

x = minimize(
    functional,
    basis.zeros(),
    jac=jacobian,
    hess=hessian,
    method="trust-constr",
    constraints=constraint,
    options={"verbose": 3},
).x

The solution is indistinguishable from that of ex01, and passes the same test.

[gdmcbain]

This technique also leads to a very succinct expression of the Stokes problem; e.g. redoing ex30 ‘Krylov-Uzawa method for the Stokes equation’ using what Ern & Guermond (2004, §4.1.3) call ‘the constraint formulation’. Basically it's the same functional, jacobian, and hessian, as for ex01, but with the bilinear form being the viscous dissipation of the velocity field. As E. & G. say (p. 179)

Note that the pressure is no longer involved.

(Here a pressure basis is constructed for the assembling of the divergence operator which appears in the constraints.)

https://github.com/gdmcbain/scikit-fem/blob/bb14737515063c6d578eaf6b24b68769e1e5ee2c/docs/examples/ex30_minimize.py#L39-L63

def functional(u: np.ndarray) -> float:
    return u @ A @ u / 2 - f @ u

def jacobian(u: np.ndarray) -> np.ndarray:
    return u @ A - f

def hessian(u: np.ndarray) -> np.ndarray:
    return A

incompressibility = LinearConstraint(B, 0, 0)
slip = csr_matrix((np.ones(D.size), (np.arange(D.size), D)), (D.size, basis["u"].N))
noslip = LinearConstraint(slip, 0, 0)

velocity = minimize(
    functional,
    basis["u"].zeros(),
    jac=jacobian,
    hess=hessian,
    method="trust-constr",
    constraints=[incompressibility, noslip],
    options={"verbose": 3},
).x

This passes the same test as ex30: 40c94aa.

[gdmcbain]

I realize now, as an aside, that ex30 could switch from the unit square

scikit-fem/docs/examples/ex30.py

Line 64 in fb649d8

mesh = MeshQuad.init_tensor(*(np.linspace(-.5, .5, 2**6),)*2)

to use the same unit circle domain as ex18

scikit-fem/docs/examples/ex18.py

Line 54 in fb649d8

mesh = MeshTri.init_circle(4)

but on a finer mesh. The original motivation for using the square rather than the disk here was that before #476, ex18 loaded a static mesh

scikit-fem/docs/examples/ex18.py

Line 54 in 7eb3016

mesh = from_file(Path(__file__).with_name("disk.json")).to_meshtri()

whereas a finer mesh was wanted for ex30 to demonstrate that the Krylov–Uzawa scheme did scale reasonably to bigger problems.

[kinnala]

Oops. Interesting findings you have here.