# Computational Mathematics and Applications Seminar

Since the legendary 1972 encounter of H. Montgomery and F. Dyson at tea time in Princeton, a statistical correspondence of the non-trivial zeros of the Riemann Zeta function with eigenvalues of high-dimensional random matrices has emerged. Surrounded by many deep conjectures, there is a striking analogyto the energy levels of a quantum billiard system with chaotic dynamics. Thanks

to extensive calculation of Riemann zeros by A. Odlyzko, overwhelming numerical evidence has been found for the quantum analogy. The statistical accuracy provided by an enormous dataset of more than one billion zeros reveals distinctive finite size effects. Using the physical analogy, a precise prediction of these effects was recently accomplished through the numerical evaluation of operator determinants and their perturbation series (joint work with P. Forrester and A. Mays, Melbourne).

$$

\def\curl#1{\left\{#1\right\}}

\def\vek#1{\mathbf{#1}}

$$

lll-posed problems arise often in the context of scientific applications in which one cannot directly observe the object or quantity of interest. However, indirect observations or measurements can be made, and the observable data $y$ can be represented as the wanted observation $x$ being acted upon by an operator $\mathcal{A}$. Thus we want to solve the operator equation \begin{equation}\label{eqn.Txy} \mathcal{A} x = y, \end{equation} (1) often formulated in some Hilbert space $H$ with $\mathcal{A}:H\rightarrow H$ and $x,y\in H$. The difficulty then is that these problems are generally ill-posed, and thus $x$ does not depend continuously on the on the right-hand side. As $y$ is often derived from measurements, one has instead a perturbed $y^{\delta}$ such that ${y - y^{\delta}}_{H}<\delta$. Thus due to the ill-posedness, solving (1) with $y^{\delta}$ is not guaranteed to produce a meaningful solution. One such class of techniques to treat such problems are the Tikhonov-regularization methods. One seeks in reconstructing the solution to balance fidelity to the data against size of some functional evaluation of the reconstructed image (e.g., the norm of the reconstruction) to mitigate the effects of the ill-posedness. For some $\lambda>0$, we solve \begin{equation}\label{eqn.tikh} x_{\lambda} = \textrm{argmin}_{\widetilde{x}\in H}\left\lbrace{\left\|{b - A\widetilde{x}} \right\|_{H}^{2} + \lambda \left\|{\widetilde{x}}\right\|_{H}^{2}} \right\rbrace. \end{equation} In this talk, we discuss some new strategies for treating *discretized *versions of this problem. Here, we consider a discreditized, finite dimensional version of (1), \begin{equation}\label{eqn.Axb} Ax = b \mbox{ with } A\in \mathbb{R}^{n\times n}\mbox{ and } b\in\mathbb{R}^{n}, \end{equation} which inherits a discrete version of ill conditioning from [1]. We propose methods built on top of the Arnoldi-Tikhonov method of Lewis and Reichel, whereby one builds the Krylov subspace \begin{equation}

\mathcal{K}_{j}(\vek A,\vek w) = {\rm span\,}\curl{\vek w,\vek A\vek w,\vek A^{2}\vek w,\ldots,\vek A^{j-1}\vek w}\mbox{ where } \vek w\in\curl{\vek b,\vek A\vek b}

\end{equation}

and solves the discretized Tikhonov minimization problem projected onto that subspace. We propose to extend this strategy to setting of augmented Krylov subspace methods. Thus, we project onto a sum of subspaces of the form $\mathcal{U} + \mathcal{K}_{j}$ where $\mathcal{U}$ is a fixed subspace and $\mathcal{K}_{j}$ is a Krylov subspace. It turns out there are multiple ways to do this leading to different algorithms. We will explain how these different methods arise mathematically and demonstrate their effectiveness on a few example problems. Along the way, some new mathematical properties of the Arnoldi-Tikhonov method are also proven.

The FEniCS system [1] allows the description of finite element discretisations of partial differential equations using a high-level syntax, and the automated conversion of these representations to working code via automated code generation. In previous work described in [2] the high-level representation is processed automatically to derive discrete tangent-linear and adjoint models. The processing of the model code at a high level eases the technical difficulty associated with management of data in adjoint calculations, allowing the use of optimal data management strategies [3].

This previous methodology is extended to enable the calculation of higher order partial differential equation constrained derivative information. The key additional step is to treat tangent-linear

equations on an equal footing with originating forward equations, and in particular to treat these in a manner which can themselves be further processed to enable the derivation of associated adjoint information, and the derivation of higher order tangent-linear equations, to arbitrary order. This enables the calculation of higher order derivative information -- specifically the contraction of a Kth order derivative against (K - 1) directions -- while still making use of optimal data management strategies. Specific applications making use of Hessian information associated with models written using the FEniCS system are presented.

[1] "Automated solution of differential equations by the finite element method: The FEniCS book", A. Logg, K.-A. Mardal, and G. N. Wells (editors), Springer, 2012

[2] P. E. Farrell, D. A. Ham, S. W. Funke, and M. E. Rognes, "Automated derivation of the adjoint of high-level transient finite element programs", SIAM Journal on Scientific Computing 35(4), C369--C393, 2013

[3] A. Griewank, and A. Walther, "Algorithm 799: Revolve: An implementation of checkpointing for the reverse or adjoint mode of computational differentiation", ACM Transactions on Mathematical Software 26(1), 19--45, 2000

A posteriori error estimators are a key tool for the quality assessment of given finite element approximations to an unknown PDE solution as well as for the application of adaptive techniques. Typically, the estimators are equivalent to the error up to an additive term, the so called oscillation. It is a common believe that this is the price for the `computability' of the estimator and that the oscillation is of higher order than the error. Cohen, DeVore, and Nochetto [CoDeNo:2012], however, presented an example, where the error vanishes with the generic optimal rate, but the oscillation does not. Interestingly, in this example, the local $H^{-1}$-norms are assumed to be computed exactly and thus the computability of the estimator cannot be the reason for the asymptotic overestimation. In particular, this proves both believes wrong in general. In this talk, we present a new approach to posteriori error analysis, where the oscillation is dominated by the error. The crucial step is a new splitting of the data into oscillation and oscillation free data. Moreover, the estimator is computable if the discrete linear system can essentially be assembled exactly.

In many applications requiring the solution of a linear system Ax=b, the matrix A has been shown to have a low-rank property: its off-diagonal blocks have low numerical rank, i.e., they can be well approximated by matrices of small rank. Several matrix formats have been proposed to exploit this property depending on how the block partitioning of the matrix is computed.

In this talk, I will discuss the block low-rank (BLR) format, which partitions the matrix with a simple, flat 2D blocking. I will present the main characteristics of BLR matrices, in particular in terms of asymptotic complexity and parallel performance. I will then discuss some recent advances and ongoing research on BLR matrices: their multilevel extension, their use as preconditioners for iterative solvers, the error analysis of their factorization, and finally the use of fast matrix arithmetic to accelerate BLR matrix operations.

In this talk I will present two recent findings concerning the preconditioning of elliptic problems. The first result concerns preconditioning of elliptic problems with variable coefficient K by an inverse Laplacian. Here we show that there is a close relationship between the eigenvalues of the preconditioned system and K.

The second results concern the problem on mixed form where K approaches zero. Here, we show a uniform inf-sup condition and corresponding robust preconditioning.

Standard quadratic programs have numerous applications and play an important role in copositivity detection. We consider reformulating a standard quadratic program as a mixed integer linear programming (MILP) problem. We propose alternative MILP reformulations that exploit the specific structure of standard quadratic programs. We report extensive computational results on various classes of instances. Our experiments reveal that our MILP reformulations significantly outperform other global solution approaches.

This is joint work with Jacek Gondzio.