4
$\begingroup$

As a model for a nonlinear least-squares problem with a large linear part problem, consider $$ \Delta u = 0 \quad\text{in } \Omega,\\ n\cdot\nabla u = 0 \quad \text{on } \Gamma,\\ (u(x_i) - u(y_i))^2 = 1 \quad \text{ for } i \in\{1,\dots,10\}, $$ i.e., on a domain $\Omega$ we're looking for a function $u$ that's somewhat smooth and for a bunch of given coordinate pairs, we want the function values to differ by 1. This system is overdetermined and will generally have no exact solution, but one could treat it with your favorite nonlinear least-squares method.

Now, when discretizing the problem with, e.g., FEM and quite a fine mesh, the linear operator $\Delta_h$ will be much larger than the 10 nonlinear equations at the end, and also one can solve linear systems with $\Delta_h$ really fast.

Is it possible to exploit this fact?

$\endgroup$
6
$\begingroup$

If you change variables to optimize for the residual of the linear part, then the Hessian will be a low-rank update to the identity. Then L-BFGS would work very well. Specifically, your problem takes the form $$ \min_x \frac{1}{2}\|Ax-b\|^2 + \frac{\mu}{2}\|g(x)\|^2 $$ where $Ax=b$ is the linear PDE and $g$ is the nonlinear part, and $\mu$ is a tradeoff parameter. Let $r:=Ax-b$. Then the problem becomes $$ \min_r \frac{1}{2}\|r\|^2 + \frac{\mu}{2}\|g(A^{-1}(r+b))\|^2. $$ The gradient of the problem is $$ \nabla f(r) = r + \mu A^{-T} G^T g(A^{-1}(r+b)); $$ the (Gauss-Newton) Hessian for this problem is $$ H = I + \mu A^{-T}G^TG A^{-1}, $$ where $G$ is the Jacobian of $g$, which has rank $10$. Hence the Hessian is a rank-10 update to the identity. So, on this modified problem L-BFGS will probably converge in around 10 iterations.

Each iteration requires solving a system of the form $Ax=z$ to evaluate the objective function and $A^Tx=w$ to compute the gradient.

Edit: Heres a simple python mockup. Seems to work.

from numpy import ones, zeros, arange, dot, sqrt
from numpy.random import choice, randn
from numpy.linalg import norm
from scipy.sparse import eye, diags, kron, coo_matrix, bmat
from scipy.sparse.linalg import factorized
from scipy.optimize import minimize
from matplotlib.pyplot import matshow, show

# See: https://scicomp.stackexchange.com/a/29205/1502

n=40
N=n*n
dx = 1./n

m=10
inds = choice(N,2*m)
xx_inds = inds[:m]
yy_inds = inds[m:]

mu = 1./(m*(dx**2)) # how strongly to force the differences being 1 vs forcing the PDE being solved

# A: 2D Neumann Laplacian
A1D = diags([-ones(n - 1), 2*ones(n), -ones(n - 1)], [-1, 0, 1], shape=(n, n)).tocsr()
A1D[0,0] = 1
A1D[-1,-1] = 1
I=eye(n)
A0 = (kron(A1D, I) + kron(I, A1D)).tocsr()
A = bmat([[A0, ones([N,1])],[ones([1,N]), zeros(1)]])
solve_A = factorized(A)
solve_At = solve_A # symmetric

B = (coo_matrix((ones(m),(arange(m),xx_inds)), shape=(m,N+1))
     - coo_matrix((ones(m),(arange(m),yy_inds)), shape=(m,N+1))).tocsr()
def g(u):
    return (B*u)**2 - 1.

def G(u):
    return 2*diags(B*u,0)*B

def objective(r):
    return 0.5*norm(r)**2 + mu*0.5*norm(g(solve_A(r)))**2

def gradient(r):
    u = solve_A(r)
    return r + mu*solve_At(G(u).T * g(u))

# Check gradient with finite differences
r0 = randn(N+1)
s = 1e-7
dr = randn(N+1)
r1 = r0 + s*dr

j0 = objective(r0)
j1 = objective(r1)
dj_diff = (j1-j0)/s
dj = dot(gradient(r0),dr)
gradient_err = norm(dj - dj_diff)/norm(dj_diff)
print 's=', s, ', gradient_err=', gradient_err

# Solve with L-BFGS
result = minimize(objective, ones(N+1), jac=gradient, method='L-BFGS-B', options={'maxiter':50})
r_optimal = result['x']
u_optimal = solve_A(r_optimal)
g_optimal = g(u_optimal)
num_iter = result['nit']
print 'num_iter=', num_iter
print 'norm(r_optimal*(dx**2))=', norm(r_optimal*(dx**2))
print 'norm(g_optimal)/m=', norm(g_optimal)/m
print 'u_optimal[xx_inds] - u_optimal[yy_inds]=', u_optimal[xx_inds] - u_optimal[yy_inds]
matshow(u_optimal[:N].reshape([n,n]))
show()

Edit2: Enforcing Neumann condition

Let $\{\phi\}_{i=1}^N$ be a set of finite element basis functions, let $$A^0_{ij}:= \int_\Omega \nabla\phi_i \cdot \nabla \phi_j ~dx,$$ and let $$\mathbb{1}_j:= \int_\Omega \phi_j ~dx.$$ Then at the discrete level solving the pure Neumann Laplacian problem with average value zero constraint is equivalent to solving the following saddle point system: $$\underbrace{\begin{bmatrix}A^0 & \mathbb{1} \\ \mathbb{1}^T & 0\end{bmatrix}}_{A}\underbrace{\begin{bmatrix}u \\ \lambda\end{bmatrix}}_{x} = 0.$$ The matrix $A$ is invertible since $A^0$ is invertible on $\text{ker}(\mathbb{1}^T)$. For some background on invertibility of saddle point systems I recommend this: https://arxiv.org/abs/1202.3330 .

If you have some specialized or black box solver that solves the problem $$A^0 \widehat{u} = f, \quad\quad \text{and} \quad\quad \mathbb{1}^T \widehat{u}=c^0,$$ you can use it to solve problems of the form $$\begin{bmatrix}A^0 & \mathbb{1} \\ \mathbb{1}^T & 0\end{bmatrix}\begin{bmatrix}u \\ \lambda\end{bmatrix} = \begin{bmatrix}f \\ c\end{bmatrix}$$ by setting $u = \widehat{u}+(c-c^0)$ and $\lambda = (A^0 u - f)_j/\mathbb{1}_j$, any component $j$ will do; Since the first row block equation is satisfied, the quantity $(A^0 u - f)_j/\mathbb{1}_j$ will be the same for all $j$.

$\endgroup$
  • $\begingroup$ Thanks for this nice answer! One question: In the example, A as a Neumann-Laplacian is rank-deficient, so $A^{-1}$ isn't well-defined. Not sure what factorized(A) does in the code, but I'd assume it'll be flawed by round-off. Perhaps this can be saved when avoiding the (admittedly nice) reformulation in $r$ by using the Hessian of the original problem, $H=A^TA + G^TG$. This is also a low-rank update on a known matrix. Not sure if $A^TA$ can be given to scipy's BFGS though. $\endgroup$ – Nico Schlömer Apr 1 '18 at 11:23
  • $\begingroup$ Aha, the code uses a Dirichlet-Laplacian which is why you don't see a break-down. $\endgroup$ – Nico Schlömer Apr 1 '18 at 12:13
  • $\begingroup$ The current version uses the Neumann laplacian with average value zero constraint enforced via Lagrange multipliers: $A=\begin{bmatrix}A_0 & 1 \\ 1^T & 0\end{bmatrix}$. The first version I uploaded did use dirichlet bcs, which was on here for about an hour yesterday before I changed it to Neumann in edit 3. $\endgroup$ – Nick Alger Apr 1 '18 at 13:42
  • $\begingroup$ Well, adding another constraint will only work if the system is consistent, i.e., if the right hand side has no components in the nullspace of $A$. Is this guaranteed here? $\endgroup$ – Nico Schlömer Apr 1 '18 at 14:24
  • 1
    $\begingroup$ I'm fairly confident that it is correct as written; you can run the code yourself and see that it works. But it will be hard to explain everything in the comments. Let me edit the post a bit to try to clarify. $\endgroup$ – Nick Alger Apr 1 '18 at 18:03
4
$\begingroup$

If you want to solve such problems with least squares, you need to first address two things:

  • You are interested in a continuous problem, not one particular discretization of it. Your statement "quite a fine mesh, the linear operator $\Delta$_h will be much larger than the 10 nonlinear equations at the end" shows that you are somehow aware of the problem, but you are not drawing the right conclusion. The correct conclusion is that any formulation for the least squares minimization must be posed in a way that is independent of the discretization. In particular, this will mean that you need to pose it in a functional analytic, rather than a finite dimensional setting. An example would be to describe the residual $r=\Delta u$ in term of the $H^{-1}$ norm; you'd then have to find a way to compute this norm for your finite dimensional approximation.

  • In this context, it is important to understand that the solution of the Laplace equation are in $H^1$, and consequently evaluating the solution at individual points is not a well-defined operation. This means that your least-squares problem is not well posed: it is conceivable that you can find a sequence of functions that all satisfy the Laplace equation and for which the constraint becomes smaller and smaller, but the limit is not in $H^1$.

$\endgroup$
  • $\begingroup$ Thanks Wolfgang for the comment. You're right: The problem should better be posed as an optimization problem of least-squares type, i.e., $\|\Delta u\|^2 + \|g(u)\|^2 \to \min$. About actions on points not being well-defined: Fair enough, let's get the trace spaces out. Since this is usually not mentioned (when talking about Dirichlet conditions, for example), I thought I just leave it here, too. Honestly I think it'd overcomplicate the question. $\endgroup$ – Nico Schlömer Apr 9 '18 at 13:21

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.