## Timetable

This is the preliminary timetable. If you have any remarks, please contact the organizers (email address at the bottom of this page).

You can also take a look at the preliminary book of abstracts (1.7 MB) which contains some more technical information about the lecture theatre setup.

You can hover over the names of the speakers to see the title of the talk; or you can click on the name of the speaker to read the abstract of that talk. You can use the back button to arrive back at the timetable.

We are using jsMath to typeset the mathematics on these pages, so you might see a message from jsMath that it is going to use images for the fonts if it can not find the TeX fonts on your own machine.

**Note: Internet Explorer is particularly slow rendering the mathematics
on this page, just be patient...**

Monday | Tuesday | Wednesday | Thursday | Friday | |

09:00-09:30 | Registration/Opening | Corless | Niu | Dick | Bailey |

09:30-10:00 | Hackbusch | Garcke | Ziveyi | Baldeaux | Borwein |

10:00-10:30 | Ullrich | Fasshauer | Kang | Gnewuch | Chernov |

Coffee/tea | |||||

11:00-11:30 | Leopardi | Pflüger | Pommier | Joe | Sørevik |

11:30-12:00 | Le Gia | Lamichhane | Gerstner | Nuyens | Prizes and closing |

Lunch | |||||

14:00-14:30 | Heinrich | Mhaskar | Sloan | Yserentant | |

14:30-15:00 | Schwab | Iske | Griebel | Schneider | |

15:00-15:30 | Surana | Sainudiin | Kucherenko | Heber | |

Coffee/tea | |||||

16:00-16:30 | Jakeman | Daum | Papageorgiou | Hegland | |

16:30-17:00 | Wozniakowski | Powell | Kuo | Jahnke | |

- | Reception+BBQ | Conference dinner |

We have a reception + barbecue on Monday evening at around 6 pm on the university campus at JG's.

The conference dinner is on Wednesday evening at 7 pm at the Real Thai restaurant at the Spot in Randwick. (This is walking distance from the university and so hopefully also from where you are staying.)

## Abstracts in order of appearance

## A new scheme for the tensor representation

A representation scheme for tensors is described which is more flexible than the usual ones. It is more powerful in the sense that r-term representations, subspace representations (Tucker) as well as sparse grid representations can be mimiced by this scheme. One advantage is stability. Furthermore, the important truncation procedure can be performed by standard algebraic tools.

## Tensor products of Sobolev-Besov spaces and approximation from the hyperbolic cross

Besov as well as Sobolev spaces of dominating mixed smoothness on \mathbb{R}^d are shown to be tensor products of d spaces defined on \mathbb{R}. Using such a representation we can transfer wavelet isomorphisms from the one-dimensional to the d-dimensional situation. Based on that we point out consequences for hyperbolic cross and best m-term approximation, in particular with respect to tensor product splines.

This work is based on a previous work by the authors, see [1].

*East Journal on Approximations*13(4):387–425, 2007.

## Quadrature using sparse grids on products of spheres

This talk will examine the use of sparse grids for quadrature on products of unit spheres in \mathbb{R}^3 and compare it with the component-by-component method of Hesse, Kuo and Sloan [1].

*Journal of Complexity*, 23(1):25–51, 2007.

## Preconditioners for interpolation on the unit sphere using spherical basis functions

The interpolation problem on the unit sphere using scattered data has many applications in geodesy and earth science in which the sphere is taken as a model for the earth. When spherical basis functions are used to construct the interpolant, the underlying linear systems are ill-conditioned. In this talk, we present an additive Schwarz preconditioner to accelerate the solution process. An estimate for the condition number of the preconditioned system will be discussed. Numerical experiments using MAGSAT satellite data will be presented.

## Stochastic approximation of smooth functions

We study the approximation of functions from a Sobolev space W_p^r(Q) in the norm of another Sobolev space W_q^s(Q) (and, in particular, in the norm of L_q(Q)) by randomized linear algorithms of the form

with n a natural number, (\Omega,\Sigma, {\mathbb P}) a probability space, x_{i,\omega}\in Q and \psi_{i,\omega}\in W_q^s(Q). Here Q\subset{\mathbb R}^d is a bounded Lipschitz domain, 1\le p,q \le \infty, r is a non-negative integer, and s is an arbitrary integer. The parameters are assumed to be such that the embedding J:W_p^r(Q)\to W_q^s(Q) exists. The error is defined as

where \mathcal{B}_{W_p^r(Q)} denotes the unit ball of W_p^r(Q). For fixed n the minimal error among all approximations of the form (1) is

In this talk we survey known results about the minimal error of randomized and also of deterministic approximation. Then we present new results establishing the order of the optimal error, as a function of n, for all p,q, r, s as above, this way solving two problems posed in a recent monograph by Novak and Wozniakowski [1].

It turns out that there is a number of situations where randomized approximation yields better rates than deterministic one. Of particular interest is the case s<0 in view of applications to the complexity of elliptic equations in weak form.

## Sparse tensor FEM for elliptic SPDEs

We consider the Finite Element Solution of second order elliptic problems in a physical domain D\subset \mathbb{R}^d with spatially inhomogeneous random coefficients.

We present convergence rates and complexity estimates for sparse Galerkin semidiscretization in the probability domain of the random solution. It is parametric in the first M Karhúnen-Loève (KL) variables of the input data [1].

Two cases are distinguished:

(i) Exponential decay of the input's KL expansion based on [2]

and

(ii) algebraic decay of the input's KL expansion.

In (i), a "polynomial chaos" type Galerkin discretization
is shown to yield spectral convergence rates * in terms
of N_\Omega, the number of deterministic elliptic problems
to be solved*. In (ii), first
approximation rates in terms of N_\Omega are available.

Finally, in

(iii) ongoing work [3,4,5] on the total complexity vs. accuracy of (adaptive) tensor Galerkin discretizations in both, stochastic as well as in the deterministic domain D will be addressed.

Sufficient conditions on the joint pdf's of the random field input to ensure better complexity than with (Quasi) Monte Carlo in the probability domain and with Galerkin discretization in D will be identified and implementational issues will be addressed in each case.

**217**(2006), 100-122.

## Scalable uncertainty quantification in complex dynamic networks

The issue of management of uncertainty for robust system operation is of interest in a large family of complex networked systems. Such systems typically involve a large number of heterogeneous, connected components, whose dynamics is affected by possibly an equally large number of parameters. Uncertainty quantification (UQ) methods like probabilistic collocation suffer from the curse of dimensionality and become practically infeasible when applied to networks as a whole. However, many networks of interest (e.g. power systems, biological networks) are often composed of weakly interacting subsystems. We propose an iterative UQ approach that exploits these weak connections to overcome the dimensionality curse and radically accelerate uncertainty propagation in large systems. This approach relies on integrating network decomposition techniques and waveform relaxation scheme with UQ methods like probabilistic collocation and polynomial chaos.

Network decomposition can be realized by spectral graph theoretic techniques to identify weakly interacting subsystems. Waveform relaxation, a parallelizable iterative method, on the other hand, exploits this decomposition and evolves each subsystem forward in time independently but coupled with the other subsystems through their solutions from the previous iteration. At each waveform relaxation iteration we propose to apply probabilistic collocation at subsystem level and then use polynomial chaos to propagate the uncertainty among the subsystems. Since UQ methods are applied to relatively simpler subsystems which typically involve a few parameters, this renders a scalable iterative approach to UQ in networks. We illustrate this iterative scheme on a power system network and numerically analyze its convergence properties.

This work was in part supported by DARPA DSO (Dr. Cindy Daniell PM) under AFOSR contract FA9550-07-C-0024 (Dr. Fariba Fahroo PM).

## Investigating the performance of high-order stochastic collocation methods

Differential equations are used to model a wide range of systems and processes that are subject to a wide range of uncertainty in initial and boundary conditions, model coefficients, forcing terms and geometry. The effects of such uncertainty should be traced through the system thoroughly enough to allow one to evaluate their effects on prediction of model outputs. The large amount of computational time needed to solve the underlying deterministic equations limits the effectiveness of traditional Monte-Carlo techniques. Stochastic collocation (SC) has recently emerged as more efficient alternative. SC utilise sparse grid interpolation, to obtain a functional approximation of the solution, and and sparse grid quadrature to easily evaluate important statistics such as mean and variance. SC has the ability to deal with steep non-linear dependence of the solution on random model data and achieves a much faster rate of convergence than traditional Monte Carlo methods when the solutions possess sufficient smoothness in random space. However, when the regularity of the solution decreases, the rate of SC deteriorates. Multi-Element Probabilistic Collocation and dimension-adaptive SC attempt to addressed this limitation. Here compare the performance of these two methods with traditional global SC on a number of benchmark problems.

## On intractability of approximation of infinitely differentiable multivariate functions

Department of Applied Mathematics, University of Warsaw, Poland

We prove that the L_\infty approximation of C^\infty functions defined on [0,1]^d is intractable and suffers from the curse of dimensionality. This is done by showing that the minimal number of linear functionals needed to obtain for an algorithm with worst case error at most \varepsilon\in(0,1) is exponential in d. This holds despite the fact that the rate of convergence is infinite.

## Pseudospectra of matrix polynomials expressed in alternative bases

Spectra and pseudospectra of matrix polynomials are of interest in geometric intersection problems, vibration problems, and analysis of dynamical systems. In this talk we consider the effect of the choice of polynomial basis on the pseudospectrum and on the conditioning of the spectrum of regular matrix polynomials. In particular, we consider the direct use of the Lagrange basis on distinct interpolation nodes, and give a geometric characterization of "good" nodes. We also give some tools for computation of roots at infinity via a new, natural, reversal. The principal achievement of the paper is to connect pseudospectra to the well-established theory of Lebesgue functions and Lebesgue constants, by separating the influence of the scalar basis from the natural scale of the matrix polynomial, which allows many results from interpolation theory to be applied.

This talk is based on previous (one-dimensional) work by the authors, see [1]. There is some reason to hope that it would be valuable to extend this work to higher dimensions (other than by tensor products, which seem hopelessly expensive), because the conditioning results of [2] seem as extendable as those of [3].

*Mathematics and Computer Science*,

**1**2

*pp.*353–374 December 2007.

*Proceedings of SYNASC, Timisoara*,

*pp.*141–153. September 2004.

*Advances in Computational Mathematics*,

**20**

*pp.*149–159, January 2004.

## Machine learning with sums of separable functions

Institut für Mathematik, Technische Universität Berlin, Germany

We present an algorithm for learning (or estimating) a function of many variables from scattered data. The function is approximated by a sum of separable functions, following the paradigm of separated representations. The central fitting algorithm is linear in both the number of data points and the number of variables, and thus is suitable for large data sets in high dimensions. This is valid for the least squares error as well as loss functions more suited for classification like the negative log likelihood or the huberised hinge loss.

We give numerical evidence for the utility of these representations for classification and regression.

## "Optimal" scaling and stable computation of meshfree kernel methods

Meshfree reproducing kernel methods are frequently used in the recovery of an unknown function given by a set of values sampled at a (scattered) set of points in \mathbb{R}^s. Even if we decide on a particular reproducing kernel Hilbert space as our approximation space we usually need to determine one or more *shape parameters* of the kernel –- the variance of a Gaussian kernel may serve as a typical example. We will report on recent studies performed in our group at IIT on determining such shape parameters in an "optimal" way.

Closely connected to the determination of "optimal" shape parameters is the need for stable computation of the kernel approximation since the best results are often obtained when the corresponding system of linear equations is ill-conditioned. One of the best-known approaches to deal with such problems involves the singular value decomposition which was shown to be computationally feasible by Gene Golub, William Kahan and Christian Reinsch in the late 1960s and early 1970s. Another technique –- perhaps much less known –- is due to James Riley [1]. We will explain Riley's algorithm (which we just recently became aware of) and compare it to SVD-based implementations.

## Approximation of non-smooth functions with spatially adaptive sparse grids in classification tasks

The classification problem arising in Data Mining – learning an unknown function by a set of scattered training data – requires the representation of functions on a typically high dimensional feature space. Recent developments show that spatially adaptive methods – implementing a successive grid refinement – allow tackling currently up to 160 dimensions. To understand the underlying mechanisms, suitable basis functions, adequate refinement criterias and error measures have to be developed and examined.

We show techniques which enable us to reach from about 40 to more than 160 dimensions, with a main focus on considerations on the choice of basis functions and the regularization operator. Results on typical real-world datasets as well as artificial ones are demonstrated. To gain further insight, we started the examination of the approximation of non-smooth indicator functions. First results will be shown.

This work is based on previous work by the authors, see [1] and [2].

*Modelling, Simulation and Optimization of Complex Processes, Proc. Int. Conf. HPSC, Hanoi, Vietnam*, pages 121–130. Springer, March 2006.

*ICCS 2007*, volume 4487 of

*LNCS*, pages 708–715. Springer, May 2007.

## A new mixed finite element method for approximating thin plate spline

We consider a mixed finite element method for the discretization of thin plate splines. The mixed formulation is obtained by introducing the gradient of the corresponding function as an additional unknown. An efficient numerical scheme is obtained using a biorthogonal system to discretize the gradient of the corresponding function. The curse of dimensality can be overcome by using sparse grids and combination techniques.

## Function approximation on data defined manifolds

Many applications require approximation of functions based on data that may be assumed to lie on a small dimensional manifold embedded in high dimensional ambient Euclidean space, but the manifold is not otherwise known. We will discuss some recent work in this direction.

Part of this work is joint with M. Maggioni and F. Filbir.

## Analysis of high-dimensional signal data by manifold learning and convolution transforms

Recent advances in nonlinear dimensionality reduction and manifold learning have provided novel methods in the analysis of high-dimensional signals. In this problem, a very large data set U \subset {\mathbf R}^n of scattered points is given, where the data points are assumed to lie on a compact submanifold {\mathcal M} of {\mathbf R}^n, i.e. U \subset {\mathcal M} \subset {\mathbf R}^n. Moreover, the dimension k={\rm dim}({\mathcal M}) of {\mathcal M} is assumed to be much smaller than the dimension of the ambient space {\mathbf R}^n, k \ll n. Now, the primary goal in the data analysis through dimensionality reduction is to contruct a low-dimensional representation of U. To this end, we combine suitable techniques from manifold learning with signal transformations to construct a projection map P : {\mathbf R}^n \to {\mathbf R}^k which then outputs the desired low-dimensional representation P(U) \subset {\mathbf R}^k. But the projection map P is required to preserve instrinsic geometrical and topological properties of the manifold {\mathcal M} in order to obtain a sufficiently accurate (low-dimensional) approximation to U. Therefore, the construction of P needs particular care. In our construction of P, customized convolution filters and suitable wavelet transformations are utilized to analyze the geometric distortion of the manifold {\mathcal M}. The good performance of the resulting nonlinear dimensionality reduction method is illustrated by numerical examples concerning low-dimensional parameterizations of scale modulated signals and solutions to the wave equation at varying initial conditions.

## Auto-validating trans-dimensional rejection sampler

In Bayesian statistical inference and computationally intensive frequentist inference, one is interested in obtaining independent samples from a high dimensional, and possibly multi-modal target density. The challenge is to obtain samples from this target without any knowledge of the normalizing constant. Several approaches to this problem rely on Monte Carlo methods. One of the simplest such methods is the rejection sampler due to von Neumann. Here we introduce an auto-validating version of a trans-dimensional extension of the rejection sampler via interval analysis. We show that our rejection sampler does provide us with independent samples from a large class of target densities in a guaranteed manner. These samples along with their importance weights can be used in rigorous estimates of challenging integrals. We illustrate the efficiency of the sampler by theory and by examples in up to 10 dimensions. Our sampler is immune to the `pathologies' of some infamous densities including the witch's hat and can rigorously draw samples from piece-wise Euclidean spaces of small phylogenetic trees with different dimensions.

This work is based on previous work by the authors, see [1].

## Nonlinear filters with particle flow induced by log-homotopy

We derive and test a new nonlinear filter that uses particle flow to implement Bayes' rule with an ODE rather than a pointwise multiplication of two functions. We compare the computational complexity of our new algorithm with a carefully designed particle filter for an interesting class of smooth fully coupled (i.e., not sparse) nonlinear filter problems of increasing dimension (d = 1 to 20) for optimal estimation accuracy. The computational complexity of the new log-homotopy filter is many orders of magnitude less than the classic particle filte (for optimal estimation accuracy) for d greater than 2. The dimension of the state vector of the Markov process to be estimated is denoted by d. The flow of particles induced by the log-homotopy allows us to migrate the particles smoothly using an ODE, thereby avoiding one of the fundamental problems with particle filters, namely `degeneracy' or `particle collapse' as a result of Bayes' rule. This problem is especially severe for applications with accurate measurements and/or high dimensional state vectors. The initial condition of the ODE for our log-homotopy is the log of the unnormalized conditional density prior to the kth measurement, and the final solution of the ODE is the log of the unnormalized conditional density after the kth measurement. It turns out that a homotopy of the density itself does not work at all, owing to singularity of the initial condition, whereas a homotopy of the log of the density removes this singularity and works well. We derive the particle flow induced by the log-homotopy using the chain rule, the generalized inverse and Liouville's criterion for particle flow in physics. We use adaptive Metropolis-Hastings with optimal scaling to resample particles after computing Bayes' rule with log-homotopy.

We show new performance results for our log-homotopy particle filter, including: (1) problems with linear, quadratic & cubic nonlinearities in the measurement equations; (2) various signal-to-noise ratios (SNR = 0 dB, 10 dB & 20 dB); (3) various dimensions of measurement vector (m = 3 & m = d); (4) various dynamics (stable & unstable, parameterized by the eigenvalues of the state transition matrix); (5) movies of particle flow, including highly non-Gaussian probability densities which are not unimodal; (6) comparison with EKF estimation accuracy; (7) estimation accuracy vs. number of particles. In particular, our log-homotopy filter is an order of magnitude more accurate than the EKF for problems with quadratic measurement nonlinearities.

## How many data are needed to estimate least values of convex quadratic functions in high dimensions

Let the least value of a convex quadratic function of n variables be required, the function being specified by a subroutine that calculates its value for any vector of variables. It is elementary to solve the problem using (n+1)(n+2)/2 values, because then there are enough data to define the function completely. Numerical experiments with the NEWUOA software of the author, however, suggest that the position of the minimum in the space of the variables can be estimated to high accuracy using only of magnitude n function values when n is large. Some results of computations that give this conclusion will be presented. The possibility of establishing the discovery theoretically will also be considered.

## Evaluating expectations of functionals of Brownian motions: A multilevel idea

Pricing a path-dependent financial derivative, such as an Asian option, requires the computation of \mathbb{E}[g(B(\cdot))], the expectation of a payoff functional, g, of a Brownian motion, (B(t))_{t=0}^T. One natural finite dimensional approximation of the Brownian motion is, B^{(d)}, a Karhunen-Loève expansion of the Brownian motion truncated at d terms.

A multilevel idea involves Monte Carlo simulation with different truncated dimensions d_{l}, l = 1,2,\cdots, L. At each level l, the number of sample paths is n_l. The expectation may be decomposed as follows:

The estimator for each piece above is a Quasi-Monte Carlo estimator. At each level l, the cost of evaluating g(B^{(d_{l+1})}(\cdot))-g(B^{(d_{l})}(\cdot)) for a single sample path is assumed to be proportional to d_{l+1}. The cost of approximating the expected payoff by an average over n_l paths is proportional to N_l=n_{l}d_{l+1}, and the total computational cost is \mathcal{O}(N), where N = \sum_{l=1}^{L-1}n_{l}d_{l+1}. This idea was used by [1].

This talk investigates the worst case error as a function of both n_l and d_l for payoff functionals that arise from weighted reproducing kernel Hilbert spaces with moderate smoothness, similar to [2]. The definition of the Hilbert space allows us to write the square worst-case error in evaluating \mathbb{E}[g(B(\cdot))] as the sum of square worst-case error of the terms on the right hand side of (1). They correspond to difference of square discrepancies. The optimal choice of n_l and d_l for a fixed computation budget N depends on how quickly the importance of the higher numbered variables decays.

*Math. Comp.*,

**71**, 1641-1661, 2002.

## Fourier transform approach to pricing multiple asset American Options

This paper extends the Integral Transform Approach of [1] and [2] to the pricing of American Options written on two underlying assets. We derive the bivariate density function of the two underlying assets using Fourier Transform Techniques and solve the resulting integral equation using numerical integration techniques.

*Industrial Management Review - 1965*, pages 32–39, 1965.

*Quantitative Finance Research Centre*, paper number 118, 2004.

## The evaluation of American compound option prices under stochastic volatility using the sparse grid approach

We consider the problem of evaluating numerically American compound option prices when the dynamics of the underlying are driven by stochastic volatility, in particular following the square root process of Heston (1993) in [4]. A compound option (called the mother option) gives the holder the right, but not obligation to buy (long) or sell (short) the underlying option (called the daughter option). Geske (1979) in [2] developed the first closed-form solution for the price of a vanilla European call on a European call.

The incorporation of stochastic volatility into the price of compound options was first attempted by Han (2003) in his thesis [3]. Fouque and Han (2004) in [1] introduce a fast, efficient and robust approximation to compute the prices of compound options such as call-on-call options within the context of multi-scale stochastic volatility models. However, they only consider the case of European option on European option. Also their method relies on certain expansions so the range of validity of their approach is not entirely clear.

In this paper, we set up the partial differential equation (PDE) approach to pricing European and American-type compound options. We assume that both the mother and the daughter options may be American-type. The compound option prices are modelled as a solution of a two-pass free boundary PDE problem.

It seems computationally demanding to solve those two nested PDEs, however, we have applied a modified sparse grid combination approach which was developed by Reisinger in his PhD thesis [6] to solve high dimensional PDEs in a fast and accurate manner. Since the underlying share price has different scale characteristics compared with the levels of the volatility, it is difficult to implement the techniques in Reisinger's thesis directly. Instead, we have found that by modifying their approach slightly, namely by adding some fixed number of points to the volatility direction in each of the subspaces results in a relative "balance" in both direction and this modification produces accurate and efficient results.

We implemented our modified sparse grid combination technique to solve the free boundary PDE followed by the price of the daughter option so that we obtain the desired prices and free boundaries by further interpolation and extrapolation. In this way we obtain the initial and boundary conditions for the mother option. Next, we apply this technique again to solve the free boundary PDE followed by the prices of the mother option to obtain the prices of the compound option. In fact, we solve those PDEs in each of the subspaces on a parallel cluster, which makes the process very efficient.

We implement Monte Carlo Simulation together with the Method of Lines to evaluate the prices of American compound option as well. Comparing the two approaches, we find that the modified sparse grid combination technique works well in producing both efficient and accurate prices for the compound option under stochastic volatility dynamics.

*Journal of Computational Finance*,

**9(1)**, Fall 2005.

*Journal of Financial Economics*,

**7**, 63-81, 1979.

*Review of Financial Studies*,

**6(2)**, 327-343, 1993.

*Journal of Financial and Quantitative Analysis*,

**39(2)**, (2004), 253-275.

*SIAM Journal on Scientific Computing*,

**29(1)**, 440-458, 2007.

## High dimensional PDE's methods applied to option pricing

Laboratoire Jacques-Louis Lions, Université Paris 6, France

Recent developments have shown that it may be possible to use deterministic numerical methods for elliptic or parabolic problems in dimension d, for 4\le d\le 20: these methods are based either on sparse grids [1] or on sparse tensor product approximation spaces [2,3].

In this talk, we discuss such methods for option pricing in finance. We will focus on the difference between the two above mentioned approximation techniques and consider two practical problems for which only one these techniques works properly.

In the first part of the talk, we consider a European vanilla contract with a multivariate generalization of the one dimensional Ornstein-Ulenbeck-based stochastic volatility model. This type of model is used in practice to capture information on new products like options on variance. A relevant generalization is to assume that the underlying asset is driven by a jump process, leading to partial integro-differential equation (PIDE). The jump size depend of the volatility, so the classical Merton jump operator is replaced by

Due to the curse of dimensionality, classical deterministic methods are not competitive to Monte Carlo methods. We discuss sparse grids finite difference methods for solving the PIDE arising in this model with d=4 In particular, the numerical approximation of the integral term involves a wavelet collocation method. We discuss the consistency of the method and we present numerical experiments, both for the PDE and the PIDE problems.

In the second part of the talk, we consider a Basket option on several assets (five in our example) in the Black & Scholes model. We discuss Galerkin methods in a sparse tensor product space constructed with wavelets. We insist on the boundary conditions and the projection of the initial condition. We also give some remarks on the practical implementation problems.

*Computing*, 61(2):151–179, 1998.

*Advances of Computational Mathematics*, 4:171–206, 1995.

*Mathematical Modelling and Numerical Analysis*, 38(1):93–128, 2004.

*Parallel Algorithms for Partial Differential Equations*, volume 31 of

*Notes on Numerical Fluid Mechanics*, Braunschweig/Wiesbaden, 1991. Vieweg.

## Sparse binomial trees in computational finance

Tree or lattice methods are widely used by practitioners for the pricing of Bermudan, American and many other complex derivative securities as they are among the most intuitive, simple and flexible numerical approaches.

In this talk we focus on the pricing of options which depend on several underlying assets such as basket,rainbow, product or performance-dependent options. For multi-asset options which are only based on a small number of underlying assets, d-dimensional extensions of the binomial method have been proposed, e.g., by [1] and [2]. The main drawback of these methods is that their computational cost increases exponentially with the number of assets which prevents their application for problems with more than three or four assets. Furthermore, they often cannot guarantee that the probabilitiesof all up- and down-steps in the tree are positive.

Here, we propose a new multidimensional extension of the binomial method based on the sparse grid approach. This method uses specialtensor-products of one-dimensional trees which by construction only have positive probabilities. A sparsebinomial tree is then constructed using a certain combination of these tensor-product trees. This way, the exponential cost in the number of assets is avoided.

Numerical experiments demonstrate the usefulness of this method for the pricing of American-style multi-asset options with ten or more underlying assets. In particular, the performance of the method with respect to the smoothness of the payoff is investigated.

*The review of the financial studies*, 2(2):241-250, 1989.

*Management science*, 3:1640-1652, 1991.

## The smoothing effect of the ANOVA decomposition

In this talk we establish a smoothing property of the ANOVA decomposition, first remarked upon by Owen and Liu, namely that the early terms of the ANOVA decomposition of a continuous function f defined on [0,1]^d can be smooth even when f itself is not smooth. This may explain why quasi-Monte Carlo and sparse grid methods often exhibit a better-than-Monte Carlo rate of convergence for functions f with a gradient discontinuity, as in option pricing problems. The application to the pricing of Asian options is considered in detail.

## Dimension-wise integration of high-dimensional functions with applications to finance

We present a new general class of methods for the computation of high-dimensional integrals. The quadrature schemes result by truncation and discretization of the anchored-ANOVA decomposition. They are designed to exploit low effective dimensions and include sparse grid methods as a special case. To derive bounds for the resulting modelling and discretization errors, we introduce effective dimensions for the anchored-ANOVA decomposition. We show that the new methods can be applied in a locally-adaptive and dimension-adaptive way and demonstrate their efficiency by numerical experiments with high-dimensional integrals from finance.

## The identification of model effective dimensions using global sensitivity analysis

Modern mathematical models can exhibit very high complexities and include hundreds or even thousands of variables. Quite often such models have effective dimensions in the truncation d_T and/or superposition d_S senses much lower than their nominal dimensions. We suggested classification of models based on their effective dimensions [1]: 1) type A models with not equally important variables: d_T < n; 2) type B models with equally or almost equally important variables and with dominant low order interaction terms in ANOVA decomposition: d_S < n and d_T = n; 3) type C functions with equally or almost equally important variables and with dominant high order interaction terms in ANOVA decomposition: d_S = d_T = n. We show that global sensitivity analysis allows the estimation of the effective dimensions at reasonable computational costs. Namely, d_T can be found by calculation of the Sobol' sensitivity indices for subsets of variables [2]. d_S can be estimated by using the random sampling high dimensional model representation (RS-HDMR). We further developed and improved the RS-HDMR method by applying quasi random sampling based on Sobol' sequences.

Evaluation of the Sobol' sensitivity indices and coefficients of the QRS-HDMR method necessitates the computation of multidimensional integrals. We apply recently developed a high-dimensional Sobol' sequence generator with advanced uniformity properties. In the second part of our talk we consider a relationship between Properties A and A' for Sobol' sequences. We show that for dimensions higher than 3 it is possible to construct the Sobol' sequence satisfying Property A' but not Property A.

*Reliability Engineering & Systems Safety*, 2008.

*Reliability Engineering & System Safety*, 92(7):957–960, 2007.

## Tractability through increasing smoothness

We prove that some multivariate linear
tensor product problems
are tractable in the worst case setting
if they are defined as tensor products
of univariate problems with logarithmically
increasing smoothness. This is demonstrated for
the approximation problem defined over Korobov spaces
of periodic functions
and the approximation problem of certain diagonal operators.
For these two problems we show necessary and
sufficient conditions on the smoothness parameters
of the univariate problems to obtain strong polynomial tractability.
We prove
that polynomial tractability is equivalent to strong polynomial
tractability, and that weak tractability always holds for these problems.
The choice of Korobov spaces is crucial since the approximation problem
defined over Sobolev spaces of non-periodic functions
with a special choice of the norm is
*not* polynomially tractable for *all* smoothness parameters
no matter how fast they go to infinity.
Furthermore, depending on the choice of the norm
we can even lose weak tractability.

*Journal of Approximation Theory*(to appear), 2009.

*Constructive Approximation*(to appear), 2009.

## Liberating the dimension

Many recent papers considered the problem of multivariate integration, and studied the "tractability" of the problem as the dimensionality d increases. The typical question is: can we find an algorithm for which the error is bounded polynomially in d, or even independently of d? And the general answer is: yes, if we have a suitably "weighted" function space.

Here we take one step further: we consider the integration problem with
*infinitely* many variables – thus liberating the dimension – and we seek
algorithms with small error and low "cost". In particular, we assume that the
cost for evaluating a function depends on the number of "active" variables.
The cost plays a crucial role in the infinite dimensional setting and forced us
to be creative when designing the algorithms.

## Geometrical properties of generalized nets and sequences

Generalized digital nets and sequences have been introduced for the numerical integration of smooth functions using quasi-Monte Carlo rules. In this talk we present geometrical properties of such nets and sequences. The definition of our nets and sequences does not depend on linear algebra over finite fields, it only requires that the point set or sequence satisfies certain distribution properties. Generalized digital nets and sequences appear as special cases. We prove some propagation rules and give bounds on the quality parameter t.

## On the approximation of smooth functions using generalized digital nets

In this talk, we study the approximation of functions using an algorithm which firstly approximates certain Walsh coefficients of the function and secondly approximates the function using a Walsh polynomial. A similar approach has previously been used for approximating periodic functions using lattice rules (and Fourier polynomials) and Walsh series using digital nets for functions in Walsh spaces. Here we use generalized digital nets (which have recently been shown to achieve higher order convergence rates for the integration errors of smooth functions) with which we approximate functions which have mixed partial derivatives of order \alpha > 1 in each variable which are square integrable.

The approximation error is studied in the worst-case setting in the \mathcal{L}_2 norm. We also discuss tractability of our proposed approximation algorithm and present numerical examples.

## Algorithmic construction of low-discrepancy point sets via dependent randomized rounding

The star discrepancy is intimately related to high-dimensional numerical integration of certain function classes as expressed by the well-known Koksma-Hlawka inequality. A sharp version of this inequality states that the worst-case error of approximating the integral of functions from the unit ball of a distinguished Sobolev space by an equal-weight cubature is exactly the star discrepancy of the set of sample points.

The classical bounds on the smallest possible star discrepancy show a good asymptotic behavior, but the number of sample points one needs to reach the asymptotic range is exponential in the dimension d. For high dimensional numerical integration such a large number of sample points is infeasible. For this reason one is interested in pre-asymptotic bounds that exhibit a better behavior in the dimension d and in sets of sample points that satisfy these bounds.

In this talk we want to present two algorithms that generate low-discrepancy sample sets of small size. The first algorithm uses first a deterministic step based on the derandomization of dependent randomized rounding and uses randomness afterwards to vary the actual local placement of the generated sample points. The second algorithm is again based on the derandomization of dependent randomized rounding, but this time the derandomization is carried out successively for each component of the resulting sample points.

We present numerical results and discuss relations to other methods. So e.g., the proof method we used to establish that the points generated by our second algorithm satisfy a certain discrepancy bound showed in addition to be useful to prove good probabilistic bounds for point sets that result from padding quasi-Monte Carlo (QMC) by Monte Carlo (MC). Since our algorithm can be used to extend given point sets in the dimension, it may be used to produce interesting variants of padding QMC by MC.

The talk is based on joint work with Benjamin Doerr, Peter Kritzer, Friedrich Pillichshammer, and Magnus Wahlström.

## Determining the rank of a lattice rule

In some computer searches for good s-dimensional lattice rules, one makes use of the s\times s generator matrix B of the dual lattice. In such searches, one might be interested in the rank of the current lattice rule being considered. (The rank of a lattice rule may be considered to be the minimum number of sums required to write it down.)

Here we present recent results showing how this lattice rank may be determined from this integer matrix B without necessarily having to reduce B to Smith normal form. These results make use of standard matrix ranks of modulo p versions of B, where p is a prime number.

This presentation is based on work that appeared in [1].

*BIT Numerical Mathematics*

**48**, pp. 79–93, 2008.

## On higher order of convergence using lattice sequences

School of Mathematics and Statistics, University of New South Wales, Australia

We study the worst case integration error of combinations of quadrature rules in a reproducing kernel Hilbert space. We show that the error, with respect to the total number N of function evaluations used, cannot decrease faster than O(N^{-1}) in the case where several quasi-Monte Carlo rules are combined to a compound quasi-Monte Carlo rule. However, if the errors of the quadrature rules constituting the compound rule have an order of convergence O(N^{-\alpha}) for \alpha>1 then, by introducing weights, this order of convergence can be shown to be recovered for the compound rule. We apply our results to the case of lattice sequences.

## Wave functions – high-dimensional objects of (theoretically) low complexity

This talk considers the electronic Schrödinger equation of quantum theory that describes the motion of N electrons under Coulomb interaction forces in a field of clamped nuclei. Solutions of this equation depend on 3N variables, three spatial dimensions for each electron. Approximating the solutions is thus inordinately challenging, and it is conventionally believed that a reduction to simplified models, such as those of the Hartree-Fock method or density functional theory, is the only tenable approach. We indicate why this conventional wisdom need not be ironclad: the regularity of the solutions, which increases with the number of electrons, the decay behavior of their mixed derivatives, and the antisymmetry enforced by the Pauli principle contribute properties that allow these functions to be approximated with an order of complexity which comes arbitrarily close to that for a system of only two electrons.

## Size consistent solution of the electronic Schrödinger equation by the coupled cluster method

The electronic Schrödinger equation describes the stationary nonrelativistic behaviour of an quantum mechanical N electron system in the electric field generated by the nuclei. The (Projected) Coupled Cluster Method has been developed for the numerical computation of the ground state energy and wave function. It provides a powerful tool for high accuracy electronic structure calculations. The present talk aims is about a rigorous analytical treatment and convergence analysis of this method. If the discrete Hartree Fock solution is sufficiently good, the quasi-optimal convergence of the projected coupled cluster solution to the full CI solution is shown. Under reasonable assumptions also the convergence to the exact wave function can be shown in the Sobolev H1 norm. The error of the ground state energy computation is estimated by an Aubin Nitsche type approach. Although the Projected Coupled Cluster method is nonvariational it shares advantages with the Galerkin or CI method. In addition it provides size consistency, which is considered as a fundamental property in many particle quantum mechanics.

## BOSSANOVA: A bond order dissection approach for efficient electronic structure calculations

Many problems in nature sciences are of high-dimensional complexity and thus plagued by the curse of dimension. Dimension-wise decomposition approaches have long been known. One of them is the ANalysis Of VAriance (Hoeffding). One of the possible fields of employment are electronic structure calculations, where in the Born-Oppenheimer approximation a high-dimensional electronic Schrödinger equation remains to be solved. Early attempts have been called additivity models, more recent advances are the so-called High-Dimensional Model Representation, see references in [1], or general fragmentation schemes such as [2]. Usually, molecular systems are solved by the density functional or Hartree-Fock theory approximately. We propose a dimension-wise decomposition scheme that would act as a preconditioner in a Krylov subspace sense, regarding a molecular system as an ANOVA decomposition (see [3]) of fragments, of one-body, two-body, three-body fragments and so on in a superposition sense. We discuss this new method which results not only in linear scaling complexity but also in eventually an exact solution of the Schrödinger equation and demonstrate its qualities for a wide range of molecules.

*Journal of Physical Chemistry*, 110:264–272, 2006.

*Journal of Chemical Physics*, 125:104104–1 – 104104–15, 2006.

*Foundations of Computational Mathematics (FoCM05), Santander*, pages 106–161, 2006.

## On the approximation of some multi- and high-dimensional probability distributions arrising in molecular biology

In my talk I will discuss some of the computational challenges related to high dimensionality and statistical dependence. Tensor product based ideas like sparse grids, optimal combination techniques and their connections with graphical models and ANOVA decompositions will be covered.

## An adaptive wavelet method for chemical master equations

Biochemical reaction systems are traditionally modelled by ordinary differential equations representing the concentrations of the species. The reaction-rate approach, however, is insufficient if some of the species are present in a very low number of copies and small-scale stochastic fluctuations can have large-scale effects. This is the case, e.g., in gene regulatory networks where gene expression is regulated by a few activators or repressors, in viral kinetics where the fate of very few infectious individuals decides whether the infection spreads over large parts of the population, or in predator-prey systems where the presence of a few predators keeps the entire ecosystem in equilibrium. In all these examples, small changes in the population numbers of the critical species due to inherent stochastic noise can cause large-scale effects. For such applications, the appropriate description is provided by the chemical master equation. Its solution is a time-dependent probability distribution p(t,x_1, \ldots, x_d) which indicates the probability that at time t exactly x_i particles of the i-th species exist. Unfortunately, solving the chemical master equation is often a formidable task, because many applications involve multiple scales in space and time, low regularity and multi-modal, metastable solutions. The main challenge, however, is the tremendous number of degrees of freedom, which originates from the fact that the solution p(t,x) has to be approximated in each state x of a large and multi-dimensional state space.

In this talk, we present an adaptive method for chemical master equations which allows to reduce the large number of degrees of freedom considerably. The method is based on a spatial representation of the probability distribution in a sparse Haar wavelet basis, time discretization by Rothe's method, and a strategy which in each time-step selects the basis vectors required to approximate the solution. The accuracy of the method is discussed, and its efficiency is illustrated by a numerical example

## High-precision, high-dimension numerical integration (part I)

Recently numerous scientific applications have arisen that require more than the standard 64-bit (15-digit) floating-point arithmetic available in hardware on modern computer systems. These applications range from supernova simulations to the theory of nonlinear dynamics. One particularly fruitful area is in the evaluation of definite integrals to high precision. These high-precision values can then be combined with constant recognition techniques of "experimental mathematics" to produce analytic evaluations. Dozens of new and intriguing results have been found in this manner, including results in Ising theory of mathematical physics, quantum field theory and even computational biology. In spite of these successes, progress is often hampered by the extremely long run times of evaluating high-dimensional integrals.

## High-precision, high-dimension numerical integration (part II)

Faculty of Computer Science, Dalhousie University, Canada

Recently numerous scientific applications have arisen that require more than the standard 64-bit (15-digit) floating-point arithmetic available in hardware on modern computer systems. These applications range from supernova simulations to the theory of nonlinear dynamics. One particularly fruitful area is in the evaluation of definite integrals to high precision. These high-precision values can then be combined with constant recognition techniques of "experimental mathematics" to produce analytic evaluations. Dozens of new and intriguing results have been found in this manner, including results in Ising theory of mathematical physics, quantum field theory and even computational biology. In spite of these successes, progress is often hampered by the extremely long run times of evaluating high-dimensional integrals.

## Quadrature for Galerkin approximation of integral equations

Galerkin methods for integral equations in \mathbb{R}^{d} require the evaluation of integrals

where K^{(1)},K^{(2)} are d-dimensional simplices or cubes and g has a singularity at x=y. We assume that g is Gevrey smooth for x\ne y and satisfies bounds for the derivatives which allow algebraic singularities at x=y. This holds for kernel functions commonly occurring in integral equations. We construct a quadrature rule Q_{N} using N function evaluations of g which achieves exponential convergence \left|I-Q_{N}\right|\le C\exp(-rN^{\gamma}) with r,\gamma>0.

It turns out in practice that the efficient approximation of the singular integrals is a hard problem, because standard (e.g. Gaussian) quadrature methods deteriorate near the singularity, and the convergence of the resulting Galerkin method is very sensitive to quadrature errors.

Most methods in the literature rely on a very specific form of the kernel function g(x,y). Our proposed method has the advantage that it only uses pointwise evaluations of g(x,y) (no antiderivatives needed), works for integrands with noninteger singularity orders and logarithmic singularities, and uses the same algorithm in all dimensions d and all possible cases how K^{(1)} and K^{(2)} may touch.

*Computing*, 53(2):173–194, 1994.

## Multidimensional pseudo-spectral methods on lattice grids

When multidimensional functions are approximated by a truncated Fourier series, the number of terms typically increase exponentially with the dimension s. In this paper we argue for a way of reducing this number by omitting insignificant terms. We then show how lattice rules can be used for approximating the associated Fourier coefficients, allowing a similar reduction in grid-points as in expansion terms. We also show that using a lattice grids permits efficient computation of the Fourier coefficient by the FFT-algorithm. Finally we put all this together in a pseudo-spectral method and demonstrate its efficiency on the Poisson equation.