
In this paper we combine the theory of reproducing kernel Hilbert spaces with
the field of collocation methods to solve boundary value problems with special
emphasis on reproducing property of kernels. From the reproducing property of
kernels we proposed a new efficient algorithm to obtain the cardinal functions
of a reproducing kernel Hilbert space which can be apply conveniently for
multidimensional domains. The differentiation matrices are constructed and also
we drive pointwise error estimate of applying them. In addition we prove the
nonsingularity of collocation matrix. The proposed method is truly meshless and
can be applied conveniently and accurately for high order and also
multidimensional problems. Numerical results are presented for the several
problems such as second and fifth order two point boundary value problems, one
and two dimensional unsteady Burgers equations and a parabolic partial
differential equation in three dimensions. Also we compare the numerical
results with those reported in the literature to show the high accuracy and
efficiency of the proposed method

Dense kernel matrices $\Theta \in \mathbb{R}^{N \times N}$ obtained from
point evaluations of a covariance function $G$ at locations $\{ x_{i} \}_{1
\leq i \leq N} \subset \mathbb{R}^{d}$ arise in statistics, machine learning,
and numerical analysis. For covariance functions that are Green's functions of
elliptic boundary value problems and homogeneouslydistributed sampling points,
we show how to identify a subset $S \subset \{ 1 , \dots , N \}^2$, with $\# S
= O ( N \log (N) \log^{d} ( N /\epsilon ) )$, such that the zero fillin
incomplete Cholesky factorisation of the sparse matrix $\Theta_{ij} 1_{( i, j )
\in S}$ is an $\epsilon$approximation of $\Theta$. This factorisation can
provably be obtained in complexity $O ( N \log( N ) \log^{d}( N /\epsilon) )$
in space and $O ( N \log^{2}( N ) \log^{2d}( N /\epsilon) )$ in time, improving
upon the state of the art for general elliptic operators; we further present
numerical evidence that $d$ can be taken to be the intrinsic dimension of the
data set rather than that of the ambient space. The algorithm only needs to
know the spatial configuration of the $x_{i}$ and does not require an analytic
representation of $G$. Furthermore, this factorization straightforwardly
provides an approximate sparse PCA with optimal rate of convergence in the
operator norm. Hence, by using only subsampling and the incomplete Cholesky
factorization, we obtain, at nearly linear complexity, the compression,
inversion and approximate PCA of a large class of covariance matrices. By
inverting the order of the Cholesky factorization we also obtain a solver for
elliptic PDE with complexity $O ( N \log^{d}( N /\epsilon) )$ in space and $O (
N \log^{2d}( N /\epsilon) )$ in time, improving upon the state of the art for
general elliptic operators.

Lane changing is one of the most common maneuvers on motorways. Although,
macroscopic traffic models are well known for their suitability to describe
fast moving crowded traffic, most of these models are generally developed in
one dimensional framework, henceforth lane changing behavior is somehow
neglected. In this paper, we propose a macroscopic model, which accounts for
lanechanging behavior on motorway, based on a twodimensional extension of the
Aw and Rascle [Aw and Rascle, SIAM J.Appl.Math., 2000] and Zhang [Zhang,
Transport.Res.BMeth., 2002] macroscopic model for traffic flow. Under
conditions, when lane changing maneuvers are no longer possible, the model
"relaxes" to the onedimensional AwRascleZhang model. Following the same
approach as in [Aw, Klar, Materne and Rascle, SIAM J.Appl.Math., 2002], we
derive the twodimensional macroscopic model through scaling of time
discretization of a microscopic followtheleader model with driving direction.
We provide a detailed analysis of the spacetime discretization of the proposed
macroscopic as well as an approximation of the solution to the associated
Riemann problem. Furthermore, we illustrate some features of the proposed model
through some numerical experiments.

Wavelet (Besov) priors are a promising way of reconstructing indirectly
measured fields in a regularized manner. We demonstrate how wavelets can be
used as a localized basis for reconstructing permeability fields with sharp
interfaces from noisy pointwise pressure field measurements in the context of
the elliptic inverse problem. For this we derive the adjoint method of
minimizing the Besovnormregularized misfit functional (this corresponds to
determining the maximum a posteriori point in the Bayesian point of view) in
the Haar wavelet setting. As it turns out, choosing a waveletbased prior
allows for accelerated optimization compared to established
trigonometricallybased priors.

Optimization over low rank matrices has broad applications in machine
learning. For large scale problems, an attractive heuristic is to factorize the
low rank matrix to a product of two much smaller matrices. In this paper, we
study the nonconvex problem $\min_{U\in\mathcal{R}^{n\times r}} g(U)=f(UU^T)$
under the assumptions that $f(X)$ is restricted $\mu$strongly convex and
$L$smooth on the set $\{X:X\succeq 0,rank(X)\leq r\}$. We propose an
accelerated gradient method with alternating constraint that operates directly
on the $U$ factors and show that the method has local linear convergence rate
with the optimal dependence on the condition number of $\sqrt{L/\mu}$.
Globally, our method converges to the critical point with zero gradient from
any initializer. Our method also applies to the problem with the asymmetric
factorization of $X=\widetilde U\widetilde V^T$ and the same convergence result
can be obtained. Extensive experimental results verify the advantage of our
method.

This paper investigates gradient recovery schemes for data defined on
discretized manifolds. The proposed method, parametric polynomial preserving
recovery (PPPR), does not require the tangent spaces of the exact manifolds,
and they have been assumed for some significant gradient recovery methods in
the literature. Another advantage of PPPR is that superconvergence is
guaranteed without the symmetric condition which has been asked in the existing
techniques. There is also numerical evidence that the superconvergence by PPPR
is high curvature stable, which distinguishes itself from the others. As an
application, we show its capability of constructing an asymptotically exact
\textit{a posteriori} error estimator. Several numerical examples on
twodimensional surfaces are presented to support the theoretical results and
comparisons with existing methods are documented.

We study the exponent of the exponential rate of convergence in terms of the
number of degrees of freedom for various nonstandard {$p$version} finite
element spaces employing reduced cardinality basis. More specifically, we show
that serendipity finite element methods and discontinuous Galerkin finite
element methods with total degree $\mathcal{P}_p$ basis have a faster
exponential convergence with respect to the number of degrees of freedom than
their counterparts employing the tensor product $\mathcal{Q}_p$ basis for
quadrilateral/hexahedral elements, for piecewise analytic problems under
$p$refinement. The above results are proven by using a new $p$optimal error
bound for the $L^2$orthogonal projection onto the total degree $\mathcal{P}_p$
basis, and for the $H^1$projection onto the serendipity finite element space
over tensor product elements with dimension $d\geq2$. These new $p$optimal
error bounds lead to a larger exponent of the exponential rate of convergence
with respect to the number of degrees of freedom. Moreover, these results show
that part of the basis functions in $\mathcal{Q}_p$ basis {plays} no roles in
achieving the $hp$optimal error bound in the Sobolev space. The sharpness of
theoretical results is also verified by a series of numerical examples.

Algebraic multigrid (AMG) is often an effective solver for symmetric positive
definite (SPD) linear systems resulting from the discretization of general
elliptic PDEs, or the spatial discretization of parabolic PDEs. However,
convergence theory and most variations of AMG rely on $A$ being SPD. Hyperbolic
PDEs, which arise often in largescale scientific simulations, remain a
challenge for AMG, as well as other fast linear solvers, in part because the
resulting linear systems are often highly nonsymmetric. Here, a novel
convergence framework is developed for nonsymmetric, reductionbased AMG, and
sufficient conditions derived for $\ell^2$convergence of error and residual.
In particular, classical multigrid approximation properties are connected with
reductionbased measures to develop a robust framework for nonsymmetric,
reductionbased AMG.
Matrices with blocktriangular structure are then recognized as being
amenable to reductiontype algorithms, and a reductionbased AMG method is
developed for upwind discretizations of hyperbolic PDEs, based on the concept
of a Neumann approximation to ideal restriction ($n$AIR). $n$AIR can be seen as
a variation of local AIR ($\ell$AIR) introduced in previous work, specifically
targeting matrices with triangular structure. Although less versatile than
$\ell$AIR, setup times for $n$AIR can be substantially faster for problems with
high connectivity. $n$AIR is shown to be an effective and scalable solver of
steady state transport for discontinuous, upwind discretizations, with
unstructured meshes, and up to 6thorder finite elements, offering a
significant improvement over existing AMG methods. $n$AIR is also shown to be
effective on several classes of `nearly triangular' matrices, resulting from
curvilinear finite elements and artificial diffusion.

This paper presents a weakly intrusive strategy for computing a lowrank
approximation of the solution of a system of nonlinear parameterdependent
equations. The proposed strategy relies on a Newtonlike iterative solver which
only requires evaluations of the residual of the parameterdependent equation
and of a preconditioner (such as the differential of the residual) for
instances of the parameters independently. The algorithm provides an
approximation of the set of solutions associated with a possibly large number
of instances of the parameters, with a computational complexity which can be
orders of magnitude lower than when using the same Newtonlike solver for all
instances of the parameters. The reduction of complexity requires efficient
strategies for obtaining lowrank approximations of the residual, of the
preconditioner, and of the increment at each iteration of the algorithm. For
the approximation of the residual and the preconditioner, weakly intrusive
variants of the empirical interpolation method are introduced, which require
evaluations of entries of the residual and the preconditioner. Then, an
approximation of the increment is obtained by using a greedy algorithm for
lowrank approximation, and a lowrank approximation of the iterate is finally
obtained by using a truncated singular value decomposition. When the
preconditioner is the differential of the residual, the proposed algorithm is
interpreted as an inexact Newton solver for which a detailed convergence
analysis is provided. Numerical examples illustrate the efficiency of the
method.

As a counterpoint to classical stochastic particle methods for diffusion, we
develop a deterministic particle method for linear and nonlinear diffusion. At
first glance, deterministic particle methods are incompatible with diffusive
partial differential equations since initial data given by sums of Dirac masses
would be smoothed instantaneously: particles do not remain particles. Inspired
by classical vortex blob methods, we introduce a nonlocal regularization of our
velocity field that ensures particles do remain particles, and we apply this to
develop a numerical blob method for a range of diffusive partial differential
equations of Wasserstein gradient flow type, including the heat equation, the
porous medium equation, the FokkerPlanck equation, the KellerSegel equation,
and its variants. Our choice of regularization is guided by the Wasserstein
gradient flow structure, and the corresponding energy has a novel form,
combining aspects of the wellknown interaction and potential energies. In the
presence of a confining drift or interaction potential, we prove that
minimizers of the regularized energy exist and, as the regularization is
removed, converge to the minimizers of the unregularized energy. We then
restrict our attention to nonlinear diffusion of porous medium type with at
least quadratic exponent. Under sufficient regularity assumptions, we prove
that gradient flows of the regularized energies converge to solutions of the
porous medium equation. As a corollary, we obtain convergence of our numerical
blob method, again under sufficient regularity assumptions. We conclude by
considering a range of numerical examples to demonstrate our method's rate of
convergence to exact solutions and to illustrate key qualitative properties
preserved by the method, including asymptotic behavior of the FokkerPlanck
equation and critical mass of the twodimensional KellerSegel equation.

Algebraic multigrid (AMG) solvers and preconditioners are some of the fastest
numerical methods to solve linear systems, particularly in a parallel
environment, scaling to hundreds of thousands of cores. Most AMG methods and
theory assume a symmetric positive definite operator. This paper presents a new
variation on classical AMG for nonsymmetric matrices (denoted lAIR), based on a
local approximation to the ideal restriction operator, coupled with
Frelaxation. A new block decomposition of the AMG errorpropagation operator
is used for a spectral analysis of convergence, and the efficacy of the
algorithm is demonstrated on systems arising from the discrete form of the
advectiondiffusionreaction equation. lAIR is shown to be a robust solver for
various discretizations of the advectiondiffusionreaction equation, including
timedependent and steadystate, from purely advective to purely diffusive.
Convergence is robust for discretizations on unstructured meshes and using
higherorder finite elements, and is particularly effective on upwind
discontinuous Galerkin discretizations. Although the implementation used here
is not parallel, each part of the algorithm is highly parallelizable, avoiding
common multigrid adjustments for strong advection such as linerelaxation and
K or Wcycles that can be effective in serial, but suffer from high
communication costs in parallel, limiting their scalability.

QuasiMonte Carlo (QMC) method is a useful numerical tool for pricing and
hedging of complex financial derivatives. These problems are usually of high
dimensionality and discontinuities. The two factors may significantly
deteriorate the performance of the QMC method. This paper develops an
integrated method that overcomes the challenges of the high dimensionality and
discontinuities concurrently. For this purpose, a smoothing method is proposed
to remove the discontinuities for some typical functions arising from financial
engineering. To make the smoothing method applicable for more general
functions, a new path generation method is designed for simulating the paths of
the underlying assets such that the resulting function has the required form.
The new path generation method has an additional power to reduce the effective
dimension of the target function. Our proposed method caters for a large
variety of model specifications, including the BlackScholes, exponential
normal inverse Gaussian L\'evy, and Heston models. Numerical experiments
dealing with these models show that in the QMC setting the proposed smoothing
method in combination with the new path generation method can lead to a
dramatic variance reduction for pricing exotic options with discontinuous
payoffs and for calculating options' Greeks. The investigation on the effective
dimension and the related characteristics explains the significant enhancement
of the combined procedure.

One way of getting insight into nonGaussian measures, posed on infinite
dimensional Hilbert spaces, is to first obtain best fit Gaussian
approximations, which are more amenable to numerical approximation. These
Gaussians can then be used to accelerate sampling algorithms. This begs the
questions of how one should measure optimality and how the optimizers can be
obtained. Here, we consider the problem of minimizing the distance with respect
to relative entropy. We examine this minimization problem by seeking roots of
the first variation of relative entropy, taken with respect to the mean of the
Gaussian, leaving the covariance fixed. Adapting a convergence analysis of
RobbinsMonro to the infinite dimensional setting, we can justify the
application of this algorithm and highlight necessary assumptions to ensure
convergence, not only in the context of relative entropy minimization, but
other infinite dimensional problems as well. Numerical examples in path space,
showing the robustness of this method with respect to dimension, are provided.

The conditioning of implicit RungeKutta (RK) integration for linear finite
element approximation of diffusion equations on general anisotropic meshes is
investigated. Bounds are established for the condition number of the resulting
linear system with and without diagonal preconditioning for the implicit Euler
and general implicit RK methods. Two solution strategies are considered for the
linear system resulting from general implicit RK integration: the simultaneous
solution (the system is solved as a whole) and a successive solution which
follows the commonly used implementation of implicit RK methods to first
transform the system into smaller systems using the Jordan normal form of the
RK matrix and then solve them successively.
For the simultaneous solution in case of a positive semidefinite symmetric
part of the RK coefficient matrix and for the successive solution it is shown
that . If the smallest eigenvalue of the symmetric part of the RK coefficient
matrix is negative and the simultaneous solution strategy is used, an upper
bound on the time step is given so that the system matrix is positive definite.
The obtained bounds for the condition number have explicit geometric
interpretations and take the interplay between the diffusion matrix and the
mesh geometry into full consideration. They show that there are three
meshdependent factors that can affect the conditioning: the number of
elements, the mesh nonuniformity measured in the Euclidean metric, and the mesh
nonuniformity with respect to the inverse of the diffusion matrix. They also
reveal that the preconditioning using the diagonal of the system matrix, the
mass matrix, or the lumped mass matrix can effectively eliminate the effects of
the mesh nonuniformity measured in the Euclidean metric. Numerical examples are
given.

We consider solving the exterior Dirichlet problem for the Helmholtz equation
with the $h$version of the boundary element method (BEM) using the standard
secondkind combinedfield integral equations. We prove a new, sharp bound on
how the number of GMRES iterations must grow with the wavenumber $k$ to have
the error in the iterative solution bounded independently of $k$ as
$k\rightarrow \infty$ when the boundary of the obstacle is analytic and has
strictly positive curvature. To our knowledge, this result is the firstever
sharp bound on how the number of GMRES iterations depends on the wavenumber for
an integral equation used to solve a scattering problem. We also prove new
bounds on how $h$ must decrease with $k$ to maintain $k$independent
quasioptimality of the Galerkin solutions as $k \rightarrow \infty$ when the
obstacle is nontrapping.

We introduce a new family of numerical algorithms for approximating solutions
of general highdimensional semilinear parabolic partial differential equations
at single spacetime points. The algorithm is obtained through a delicate
combination of the FeynmanKac and the BismutElworthyLi formulas, and an
approximate decomposition of the Picard fixedpoint iteration with multilevel
accuracy. The algorithm has been tested on a variety of semilinear partial
differential equations that arise in physics and finance, with very
satisfactory results. Analytical tools needed for the analysis of such
algorithms, including a semilinear FeynmanKac formula, a new class of
seminorms and their recursive inequalities, are also introduced. They allow us
to prove for semilinear heat equations with gradientindependent nonlinearity
that the computational complexity of the proposed algorithm is bounded by
$O(d\,\varepsilon^{(4+\delta)})$ for any $\delta \in (0,\infty)$ under
suitable assumptions, where $d\in \mathbb{N}$ is the dimensionality of the
problem and $\varepsilon\in(0,\infty)$ is the prescribed accuracy.

Temporal or spatial structures are readily extracted from complex data by
modal decompositions like Proper Orthogonal Decomposition (POD) or Dynamic Mode
Decomposition (DMD). Subspaces of such decompositions serve as reduced order
models and define either spatial structures in time or temporal structures in
space. On the contrary, convecting phenomena pose a major problem to those
decompositions. A structure traveling with a certain group velocity will be
perceived as a plethora of modes in time or space respectively. This manifests
itself for example in poorly decaying singular values when using a POD. The
poor decay is counterintuitive, since a single structure is expected to be
represented by a few modes. The intuition proves to be correct and we show that
in a properly chosen reference frame along the characteristics defined by the
group velocity, a POD or DMD reduces moving structures to a few modes, as
expected. Beyond serving as a reduced model, the resulting entity can be used
to define a constant or minimally changing structure in turbulent flows. This
can be interpreted as an empirical counterpart to exact coherent structures. We
present the method and its application to a head vortex of a compressible
starting jet.

In this work, we study the numerical approximation of local fluctuations of
certain classes of parabolic stochastic partial differential equations (SPDEs).
Our focus is on effects for small spatiallycorrelated noise on a time scale
before large deviation effects have occurred. In particular, we are interested
in the local directions of the noise described by a covariance operator. We
introduce a new strategy and prove a Combined ERror EStimate (CERES) for the
five main errors: the spatial discretization error, the local linearization
error, the noise truncation error, the local relaxation error to steady state,
and the approximation error via an iterative lowrank matrix algorithm. In
summary, we obtain one CERES describing, apart from modelling of the original
equations and standard roundoff, all sources of error for a local fluctuation
analysis of an SPDE in one estimate. To prove our results, we rely on a
combination of methods from optimal Galerkin approximation of SPDEs, covariance
moment estimates, analytical techniques for Lyapunov equations, iterative
numerical schemes for lowrank solution of Lyapunov equations, and working with
related spectral norms for different classes of operators.

In this paper we present algorithms for an efficient implementation of the
Localized Orthogonal Decomposition method (LOD). The LOD is a multiscale method
for the numerical simulation of partial differential equations with a continuum
of inseparable scales. We show how the method can be implemented in a fairly
standard Finite Element framework and discuss its realization for different
types of problems, such as linear elliptic problems with rough coefficients and
linear eigenvalue problems.

We develop a framework that allows the use of the multilevel Monte Carlo
(MLMC) methodology (Giles2015) to calculate expectations with respect to the
invariant measure of an ergodic SDE. In that context, we study the
(overdamped) Langevin equations with a strongly concave potential. We show
that, when appropriate contracting couplings for the numerical integrators are
available, one can obtain a uniform in time estimate of the MLMC variance in
contrast to the majority of the results in the MLMC literature. As a
consequence, a root mean square error of $\mathcal{O}(\varepsilon)$ is achieved
with $\mathcal{O}(\varepsilon^{2})$ complexity on par with Markov Chain Monte
Carlo (MCMC) methods, which however can be computationally intensive when
applied to large data sets. Finally, we present a multilevel version of the
recently introduced Stochastic Gradient Langevin Dynamics (SGLD) method
(Welling and Teh, 2011) built for large datasets applications. We show that
this is the first stochastic gradient MCMC method with complexity
$\mathcal{O}(\varepsilon^{2}\log {\varepsilon}^{3})$, in contrast to the
complexity $\mathcal{O}(\varepsilon^{3})$ of currently available methods.
Numerical experiments confirm our theoretical findings.

Deep neural networks have become invaluable tools for supervised machine
learning, e.g., classification of text or images. While often offering superior
results over traditional techniques and successfully expressing complicated
patterns in data, deep architectures are known to be challenging to design and
train such that they generalize well to new data. Important issues with deep
architectures are numerical instabilities in derivativebased learning
algorithms commonly called exploding or vanishing gradients. In this paper we
propose new forward propagation techniques inspired by systems of Ordinary
Differential Equations (ODE) that overcome this challenge and lead to
wellposed learning problems for arbitrarily deep networks.
The backbone of our approach is our interpretation of deep learning as a
parameter estimation problem of nonlinear dynamical systems. Given this
formulation, we analyze stability and wellposedness of deep learning and use
this new understanding to develop new network architectures. We relate the
exploding and vanishing gradient phenomenon to the stability of the discrete
ODE and present several strategies for stabilizing deep learning for very deep
networks. While our new architectures restrict the solution space, several
numerical experiments show their competitiveness with stateoftheart
networks.

The generalized Debye source representation of timeharmonic electromagnetic
fields yields wellconditioned secondkind integral equations for a variety of
boundary value problems, including the problems of scattering from perfect
electric conductors and dielectric bodies. Furthermore, these representations,
and resulting integral equations, are fully stable in the static limit as
$\omega \to 0$ in multiply connected geometries. In this paper, we present the
first highorder accurate solver based on this representation for bodies of
revolution. The resulting solver uses a Nystr\"om discretization of a
onedimensional generating curve and highorder integral equation methods for
applying and inverting surface differentials. The accuracy and speed of the
solvers are demonstrated in several numerical examples.

In this paper we present criteria for the choice of the shape parameter c
contained in the famous radial function multiquadric. It may be of interest to
RBF people and all people using radial basis functions to do approximation.

The topic of these notes could be easily expanded into a full onesemester
course. Nevertheless, we shall try to give some flavour along with theoretical
bases of spectral and pseudospectral methods. The main focus is made on
Fouriertype discretizations, even if some indications on how to handle
nonperiodic problems via Tchebyshev and Legendre approaches are made as well.
The applications presented here are diffusiontype problems in accordance with
the topics of the PhD school.

In this paper, we prove necessary and sufficient conditions for a
hybridizable discontinuous Galerkin (HDG) method to satisfy a multisymplectic
conservation law, when applied to a canonical Hamiltonian system of partial
differential equations. We show that these conditions are satisfied by the
"hybridized" versions of several of the most commonlyused finite element
methods, including mixed, nonconforming, and discontinuous Galerkin methods.
(Interestingly, for the continuous Galerkin method in dimension greater than
one, we show that multisymplecticity only holds in a weaker sense.)
Consequently, these generalpurpose finite element methods may be used for
structurepreserving discretization (or semidiscretization) of canonical
Hamiltonian systems of ODEs or PDEs. This establishes multisymplecticity for a
large class of arbitrarilyhighorder methods on unstructured meshes.