First discretize, then optimize

4.3. First discretize, then optimize#

We want to discretize the functional

\[E(u) = \frac{1}{2}||f-Au||^2_2+\frac{\alpha}{2}||\nabla u||^2_2\]

first before calculating and solving the optimality condition of the discretized functional. We only need to discretize the gradient. In 1D, following a standard numerics course, the derivative can be approximated on a discrete grid. Let \(h>0\) be the mesh size, then the unit intervall can for example be discertized as \(\Omega_h=\{0, h,...,(N-1)h,1\}\), where \(N = \frac{1}{h}\). We can approximate the gradient in the following ways

  1. forward differences: \(u'(x_i)=u'(0+i\cdot h)\approx \left(u(x_{i+1})-u(x_i)\right)/h, \quad \text{error:}\mathcal{O}(h)\)

  2. backward differences: \(u'(x_i)=u'(0+i\cdot h)\approx \left(u(x_{i})-u(x_{i-1})\right)/h, \quad \text{error:}\mathcal{O}(h)\)

  3. central differences: \(u'(x_i)=u'(0+i\cdot h)\approx \left(u(x_{i+1})-u(x_{i-1})\right)/2h, \quad \text{error:}\mathcal{O}(h^2)\)

Boundary conditions can be incorporated as well:

  1. Periodic: \(u(x_{N+1})=u(x_1),~u(x_0)=u(x_N)\)

  2. Neumann-zero: \(u'(x_N)=u'(x_1)=0\implies u(x_0)=u(x_1)~ \& ~u(x_N)=u(x_{N+1})\)

  3. Dirichlet: \(u(x_{N+1})=u(x_0)=0\)

In higher dimensions one can easily discretize the partial derivatives as well.

Example 4.7 (Two-dimensional Dirichlet boundary condition)

For \(D\approx\nabla\) we choose forward differences:

\[\begin{aligned} (Du)_{ij}&=\left((D^1u)_{ij},(D^2u)_{ij}\right) \\ (D^1u)_{ij}&=\begin{cases} (u_{i+1,j}-0)/h \quad \text{for } ~i=0 \\ (u_{i+1,j}-u_{i,j})/h \quad \text{for } ~0

We see that the gradient as a linear operator remains a linear operator in the discrete setting and we can formulate our optimality condition. Let

\[E_h(u) = \frac{1}{2}||f-Au||^2_2+\frac{\alpha}{2}||Du||^2_2\,,\]

where \(D\) denotes the approximated gradient. Then the optimality condition yields

\[\delta E_h(u) = 0 \Leftrightarrow A^T(Au-f) + \alpha D^*Du =0\,.\]

We can either solve this linear system of equations directly or we can employ a gradient flow discretization (gradient descent) scheme by introducing an artificial time variable. This method works as well for nonlinear systems. The artificial time dependency of \(u\) writes as \(u_t\) and we have the approximated time derivative (velocity) of \(u_t\) given as \(\partial_t u_t \approx (u_{t+1}-u_t)/{\tau}\) and the gradient flow is given as

\[\begin{aligned} \partial_t u_t &= -\delta E_h(u_t)\\ \implies u_{t+1} &=u_t-\tau (A^T(Au_t -f) +\alpha D^*D u_t), \end{aligned}\]

which can be computed as an iterative scheme with \(u_0=f\).