In general I work at the intersection of numerical analysis and mathematical physics. In much of my work physical insight has allowed radical improvements in numerical methods; in turn, the development of efficient numerical methods has allowed a deeper investigation and understanding of physical phenomena, particularly in quantum systems. I am also working on a medical imaging inverse problem, where due to the complexity of the problem, many levels of approximation are necessary and physical insight is essential. For this noisy real-world application, a statistical understanding of the imaging process is important---this contrasts much existing work in the field of inverse problems.
Also here are some recent (up to Nov 2006) pages I have put up on certain special topics:
Finite Element or finite-difference type methods are impractical when there is a large number of wavelengths across the system, therefore boundary methods are necessary. I have been developing variants of the `Method of Particular Solutions' (MPS), which can be viewed as approximating trial eigenfunctions by a linear sum of basis functions which satisfy the PDE (the Helmholtz equation at some trial energy parameter E) in the domain. Approximations to the eigenvalues are then found by searching for E values where the minimum achievable boundary norm becomes very small. The MPS was largely abandoned by the numerical analysis community in the 1980's, however many exciting independent developments have occured in the physics (quantum chaos) community since then, using similar ideas, upon which I base my work.
In my thesis work I placed the domain norm and boundary norm on an equal footing, thus tackling the problem that there can exist nontrivial basis combinations which are exponentially small everywhere in the domain. I also derived formulae allowing the domain norm to be evaluated using boundary integrals alone, improving efficiency. I also gave the first account of the working of the powerful `scaling method' of Vergini, which can be seen as an accelerated MPS, whose relative efficiency grows like O(k) where k=E1/2. The scaling method produces O(k) eigenfunctions with a single matrix eigenproblem, whereas the usual MPS requires several such eigenproblems to find each eigenfunction. In typical quantum chaotic studies E~106 therefore the method is a thousand times faster than any other known method (this includes boundary integral equation [BIE] methods). Recently I have analysed the E-dependent spectral problem which arises when minimizing the boundary norm while holding the domain norm constant. By studying the spectrum of the resulting boundary operator (similar to the Dirichlet-Neumann map) before numerical discretization, I have found that a special choice of weight function allows the classical eigenvalue inclusion bounds of Kuttler-Sigillito to be improved, again by factor O(k). This relied on converting the spectral problem into a form in which perturbation theory in (E-Ej) could be used.
I am now analysing the intrinsic errors of eigenfunction and eigenvalues computed by the scaling method, using the same tool. Despite its use within the quantum chaos community, such errors have not been yet understood. In both the scaling method, and the improved inclusion bounds, the role of the matrix Qij of inner products of the eigenfunction normal derivative functions, under the boundary weight function r.n is vital. Here r is the position vector and n the unit normal. I have recently proved a quasi-orthogonality result (the matrix Q is close to diagonal, with bounds on the off-diagonal growth), independent of the domain shape. I intend to work on such questions as: can similar acceleration techniques can be applied to BIE methods for the eigenproblem? Can the scaling method be adapted for generalized (mixed) homogeneous boundary conditions? Can a star-shaped restriction on domain shapes be lifted---empirical evidence suggests the answer is yes, but can this be formalized? How can quadrature best be performed in evaluating boundary integrals of basis functions? Can the scaling method be extended to manifolds of constant negative curvature?
Basis sets for the linear space of Helmholtz solutions, necessary for the MPS and scaling methods, are not well understood. In the original MPS work regular Bessel functions were used, and the physics community has based much work on sets of plane waves. I have developed basis sets of Neumann functions (fundamental solutions) placed on a curve exterior to the domain, which appear to represent all eigenfunctions accurately, even for nonconvex domains with corners. I am analysing approximation properties of such bases, in the semiclassical limit. Can bases for re-entrant corners be developed? Are there shapes for which convential BIE will be the only option?
The Quantum Unique Ergodicity (QUE) conjecture of Rudnick-Sarnak, made in the context of negatively-curved manifolds, is that there is no exceptional subset in the QET. There is much recent interest in QUE within pure mathematics (number theory and ergodic theory). QUE has recently been proven by Lindenstrauss for automorphic forms, where the infinite number of symmetries allows an analytic approach. In billiards no such analytic machinery is known, however I have found strong numerical evidence for QUE in a uniformly-hyperbolic billiard system.
To what extent do eigenfunctions behave like random sums of plane waves (the prediction of a Random Matrix Theory [RMT] model)? Many open issues remain in this area. Correlations between eigenfunctions cannot be predicted by RMT, since they involve the dynamics of the system. For instance, the variance of off-diagonal matrix elements Aij is proven to converge to a classical autocorrelation function. I have tested this convergence and again found it suprisingly slow. The prediction for diagonal variance is related by a factor g to nearby off-diagonal variance. The RMT (GOE) prediction is g=2, however, there is little justification for this in actual systems. What if the system possesses other symmetries? This connects to recent work of Sarnak showing that g depends on level number i in automorphic forms, its value being an L-function. When there are slow decay of correlations (`bouncing ball' modes), the diagonal variance would be predicted to diverge---what is fact happens here? Even less is known about matrix element statistics in mixed systems with sticky islands other power-law decays. I have also started work with E. Bogomolny on levels-spacing statistics in triangular billiards with certain symmetries (Veech polygons).
There are several applications of this work. In quantum physics, such matrix elements form the backbone of study of driven systems and transition rates (Fermi Golden Rule). Off-diagonal matrix element variance corresponds to dissipation (heating) rate under periodic driving. In `quantum dots' (cold confined 2D electron gases) chaotic effects have been instrumental in understanding fluctuations in electrical conduction. Quantum dots are now model systems for study of spins ('spintronics') and are promising candidates for quantum computers.
Forward problem: DOT requires numerical simulation of the forward problem, that is, the passage of photons through a medium of complex geometry, in our case the human head. We approximate the time-dependent transport equation by the diffusion equation, since it can be solved much more rapidly. In tests against Monte Carlo simulations we have found diffusion to be adequate. In iterative or Bayesian approaches to inverse problems, the forward model must be called a large number of times during the optimization (fitting process) or exploration of the posterior probability density function. The scalp/skull/brain system is approximately a layered medium, for which we were surprised to find a lack of efficient diffusion solvers. I therefore have developed an efficient layered solver which makes use of the symmetry, a 1D finite-difference method in the z direction, and a logarithmic timestepping scheme motivated by the physical diffusion mode decay, giving ~0.1 sec CPU run time in relevant regimes. There are many pressing goals such as applying this timestepping scheme to more general 3D finite-difference discretizations. Since MRI data is cartesian, and triangulation in 3D is a tricky problem, we expect this to be superior to currently-used Finite Element schemes. Discretization errors due to poor representation of the smooth head surface by a cartesian grid should be overcome by the inclusion of immersed interface methods into the finite-difference scheme. The gradient (Jacobean) of model outputs with respect to inputs is vital to the optimization procedure, therefore I intend to apply adjoint differentiation techniques to these forward models.
Inverse problem: To make extraction of useful brain activity tractable, our key approach, proposed by David Boas and his group in Boston, is to use structural MRI to define the subject's head geometry, then use some kind of segmented optical model with a reduced number of parameters (of order 10 or 100 unknowns) in order to fit the optical signals. Use can be then made of information about each tissue type: optical parameters of the skull are fixed, and the scalp oxygenation correlates strongly with `global' physiological parameters. Thus changes in the brain can be isolated with more certainty. I have performed a simulated feasibility study of this segmented approach, fitting baseline brain parameters, a nontrivial problem in itself. This used a Bayesian approach, enabling the full probability density function (PDF) in the parameter space to be estimated. This gives marginal PDFs on each parameter, and correlations. Ultimately in an imaging problem with many unknowns, low resolution, and noisy (single photon counting) data collection, I argue that such statistical approaches will be the best way forward, despite the higher computational cost.
Currently in my work PDFs are explored using Markov chain Monte Carlo (MCMC) methods, such as Metropolis algorithms, however there is much to be gained by testing hybrid MCMC methods and combining with gradient-based optimization. In the limit of small noise, the posterior PDF peak approaches Gaussian, and we have established what range of imaging parameters this approximation is adequate. With a Gaussian PDF approximation, expensive MCMC can be replaced with much cheaper optimation followed by estimation of the Hessian (inverse of covariance) matrix. I have implemented this method on layered medium fits, with various numbers of layers, where the layer thicknesses themselves are unknowns. Questions such as whether estimation of skull/scalp layer thickness by DOT or via MRI give more accurate DOT optical parameter fits, and the effect of introduction of (unknown) calibration parameters on the fit are vital to answer in order to advance the reliability of the technology.