# Compression Challenges in Large Scale Partial Differential Equation Solvers
## 5 Mass Storage
### Adjoint Solutions
#### 5.1.1 PDE-Constrained Optimization  
For PDE-constrained optimization, typically \(X=Y\times U,x=(y,u)\) in the abstract problem (7), where the influence of the control \(u\) on the state \(y\) is given by the PDE. Here, \(J\) is the objective to be minimized, e.g., penalizing the deviation of \(y\) from some desired state. Especially in time-dependent problems, often the reduced form is considered: There, the PDE (7) is used to compute for a given control \(u\) the associated (locally) unique solution \(y=y(u)\). With only the control remaining as the optimization variable, the reduced problem reads \(\min_{u}j(u)\), with \(j(u):=J(y(u),u)\). Computation of the reduced gradient then leads to the adjoint equation for \(p\in Z^{\star}\)  
\[c_{y}^{\star}(p;(y,u),\varphi)=-J_{y}((y,u),\varphi), \tag{8}\]  
where \(\star\) denotes the dual operator/dual function spaces, and \(c_{y},J_{y}\) are the derivatives of \(c(y,u;\varphi),J(y,u)\) with respect to the \(y\)-component. Exemplarily, we consider optimal control of the monodomain equations on a simple 2D unit square domain \(\Omega\). This system describes the electrical activity of the heart (see, e.g., [78]), and consists of a parabolic PDE for the transmembrane voltage \(v\), coupled to pointwise ODEs for the gating variable \(w\) describing the state of ion channels:  
\[ \begin{split} v_{t}&=\operatorname{div}(\sigma \
abla v)-I_{\text{ion}}(v,w)+I_{\text{e}}&\text{in }\Omega\times(0,T)\\ w_{t}&=G(v,w)&\text{in }\Omega\times(0,T). \end{split}  \tag{9}\]  
Homogeneous Neumann boundary conditions are prescribed. The functions \(I_{\text{ion}}(v,w)\) and \(G(v,w)\) are specified by choosing a membrane model. For the optimal control problem an initial excitation in some subdomain \(\Omega_{\text{exi}}\) is prescribed. The external current stimulus is \(I_{\text{e}}(x,t)=\chi_{\Omega_{\text{e}}}(x)u(t)\), where the control \(u\) is spatially constant on a control domain \(\Omega_{\text{c}}\). Defining some observation domain \(\Omega_{\text{obs}}\), the objective is given by  
\[J(y,u)=\frac{1}{2}\left\|v\right\|_{L^{2}(\Omega_{\text{obs}}\times(0,T))}^{2} +\frac{\alpha}{2}\left\|u\right\|_{L^{2}(0,T)}^{2}, \tag{10}\]  
i.e., we aim at damping out the excitation wave. For details, see [79]. Solution of the optimization problem with inexact Newton-CG methods and lossy compression is investigated in [46; 71]; here we use the quasi-Newton method due to [PERSON], [PERSON], [PERSON], and [PERSON] (BFGS, see, e.g., [80; 81]). For time discretization, we use a linearly implicit Euler method with fixed timestep size \(dt=0.04\). Using linear finite elements, spatial adaptivity is performed individually for state and adjoint using a hierarchical error estimator [82], with a restriction to at most \(25,000\) vertices in space. The adaptively refined grids were stored using the methods from [83], which reduced the storage space for the mesh to less than 1 bit/vertex (see [34]). Lossy compression of state values, i.e., the finite element solutions \(v\) and \(w\), at all time steps, affects the accuracy of the reduced gradient computed by adjoint methods, and results in inexact quasi-Newton updates. Error analysis [22] shows that BFGS with inexact gradients converges linearly, if the gradient error \(e_{g}\) in each step fulfills  
\[\left\|e_{g}\right\|\leq\frac{\varepsilon}{\kappa(B)^{1/2}}\left\|\tilde{g}\right\| \tag{11}\]  
for \(\varepsilon<1/2\). Here, \(\kappa(B)\) is the condition number of the approximate Hessian \(B\), and \(\tilde{g}\) denotes the inexactly computed gradient. The error bound (11) allows computing adaptive compression error tolerances from pre-computed worst-case gradient error estimates analogously to [46], see [22] for details. Figure 7 shows the progress of the optimization method. For trajectory compression, different fixed as well as the adaptively chosen quantization tolerances were used. We estimate the spatial discretization error in the reduced gradient by using a solution on a finer mesh as a reference.