TECHNICAL REPORTS: (2005)
A Reduced Basis Method for Molecular Dynamics Simulation
Rachel Elisabeth Vincent-Finley
In this dissertation, we develop a method for molecular simulation based on principal component analysis (PCA) of a molecular dynamics trajectory and least squares approximation of a potential energy function. Molecular dynamics (MD) simulation is a computational tool used to study molecular systems as they evolve through time. With respect to protein dynamics, local motions, such as bond stretching, occur within femtoseconds, while rigid body and large-scale motions, occur within a range of nanoseconds to seconds. To capture motion at all levels, time steps on the order of a femtosecond are employed when solving the equations of motion and simulations must continue long enough to capture the desired large-scale motion.
To date, simulations of solvated proteins on the order of nanoseconds have been reported. It is typically the case that simulations of a few nanoseconds do not provide adequate information for the study of large-scale motions. Thus, the development of techniques that allow longer simulation times can advance the study of protein function and dynamics. In this dissertation we use principal component analysis (PCA) to identify the dominant characteristics of an MD trajectory and to represent the coordinates with respect to these characteristics. We augment PCA with an updating scheme based on a reduced representation of a molecule and consider equations of motion with respect to the reduced representation. We apply our method to butane and BPTI and compare the results to standard MD simulations of these molecules. Our results indicate that the molecular activity with respect to our simulation method is analogous to that observed in the standard MD simulation with simulations on the order of picoseconds.
December 2005
A Quantitative Study of Neuronal Calcium Signaling
Jane Hartsfield
Neurons have both a fast and slow mode of signaling. Fast signals are communicated by transmembrane voltage changes, while calcium levels within the cell communicate information on a much slower time scale. Calcium acts as a second messenger responsible for modulating neuronal excitability in many ways including the mediation of gene transcription in the cell and the sensitivity of the cell to further stimulus. I propose a means of determining calcium conductance density of the cell membrane from intracellular calcium concentration measurement data using a two step process. The first step is the inference of calcium current density from calcium concentration measurements using a least squares fit to the data. Once an estimate of the calcium current density is determined, the minimum value over time is used to determine the calcium conductance density.
I develop a voltage model of the neuron's electrical signal with ion diffusion and drift which includes voltage-gated calcium currents and calcium-dependent potassium currents. The influx of calcium resulting from the voltage model will prime the endoplasmic reticulum with calcium. A model of the dynamics of calcium induced calcium release from the endoplasmic reticulum via IP3 receptors which includes diffusion of calcium and IP3 as well as calcium buffering by the mitochondria results in a calcium wave similar to what has been observed experimentally. Finally, I use a branch structure together with IP3 generation, calcium buffers in the cytosol and ER, cell membrane calcium transports (voltage-gated calcium channels, pumps, exchangers, and store-operated channels), and ER calcium transports (IP3 receptors, ryanodine receptors, pumps, leak channels) to show that calcium waves initiate in the apical trunk at the point where the stimulated oblique branches off.
December 2005
Nonlinear Programming by Mesh Adaptive Direct Searches
M. A. Abramson, C. Audet, & J.E. Dennis, Jr.
This paper is intended not as a survey, but as an introduction to some ideas behind the class of mesh adaptive direct search (MADS) methods.
Space limitations dictate a brief description of various key topics be provided along with several references, which themselves provide further references.
The convergence theory for the methods presented here make a case for closing the gap between nonlinear optimizers and nonsmooth analysts.
However these methods are certainly not of purely theoretical interest; they are successful on difficult practical problems. To encourage further use, we give references to available implementations. MADS is implemented in the direct search portion of the MathWorks MATLAB Genetic Algorithm and Direct Search (GADS) Toolbox.
October 2005
A Software Framework for the Abstract Expression of Coordinate-Free Linear Algebra and Optimization Algorithms
William W. Symes, Anthony D. Padula, & Shannon D. Scott
Object oriented design solves a fundamental programming problem arising in scientific and engineering applications of linear algebra and optimization: the separation in code of multiple levels of abstraction naturally appearing in solution algorithms for such problems. The Rice Vector Library provides C++ classes expressing core concepts (vector, function,...) of calculus in Hilbert space with minimal implementation dependence, and standardized interfaces behind which to hide application-dependent implementation details (data containers, function objects). A variety of coordinate free algorithms from linear algebra and optimization, including Krylov subspace methods and various relatives of Newton's method for nonlinear equations and constrained and unconstrained optimization, may be expressed purely in terms of this system of classes. The resulting code may be used {\em without alteration} in a wide range of control, design, and parameter esti mation applications, in serial and parallel computing environments.
October 2005
Software Design for Simulation Driven Optimization
Anthony D. Padula
This thesis describes a flexible framework for abstract numerical algorithms which treats algorithms as objects and makes them reusable, composable, and modifiable. These algorithm objects are implemented using the Rice Vector Library (RVL) interface, decoupling the algorithmic code from the details of linear algebra and calculus in Hilbert Space. I made many improvements to the RVL design, including abstract return types for reductions. These improvements allowed me to demonstrate the breadth of this design by incorporating semantically similar objects from other packages which had significant syntatic differences to the RVL objects. By adapting other libraries, I gain access to a variety of tools, including parallel linear algebra implementations. The benefits of the algorithm framework can be seen when abstract numerical algorithms are coupled with parallel simulators without needing to modify either the algorithm or the simulator.
October 2005
A Simple Proof for Recoverability of L1-Minimization (II): the Nonnegativity Case
Yin Zhang
When using L1 minimization to recover a sparse, nonnegative solution to a under-determined linear system of equations, what is the highest sparsity level at which recovery can still be guaranteed? Recently, Donoho and Tanner discovered, by invoking classic results from the theory of convex polytopes that the highest sparsity level equals half of the number of equations. In this report, we provide a completely self-contained, yet short and elementary, proof for this remarkable result. We also connect dots for different recoverability conditions obtained from different spaces.
September 2005
A Simple Proof for Recoverability of L1-Minimization: Go Over or Under?
Yin Zhang
It is well-known by now that L1 minimization can help recover sparse solutions to under-determined linear equations or sparsely corrupted solutions to over-determined equations, and the two problems are equivalent under appropriate conditions. So far almost all theoretic results have been obtained through studying the ``under-determined side'' of the problem. In this note, we take a different approach from the ``over-determined side'' and show that a recoverability result (with the best available order) follows almost immediately from an inequality of Garnaev and Gluskin. We also connect dots with recoverability conditions obtained from different spaces.
August 2005
Convergence of Mesh Adaptive Direct Search to Second-Order Stationary Points
Mark A. Abramson & Charles Audet
A previous analysis of second-order behavior of pattern search algorithms for unconstrained and linearly constrained minimization is extended to the more general class of mesh adaptive direct search (MADS) algorithms for general constrained optimization. Because of the ability of MADS to generate an asymptotically dense set of search directions, we are able to establish reasonable conditions under which a subsequence of MADS iterates converges to a limit point satisfying second-order necessary or sufficient optimality conditions for general set-constrained optimization problems.
August 2005
Optimal Stepsize Selection for Verifying Adifor Derivatives
M. Fagan
Divided differences are the most natural way to verify derivative values computed by Adifor. Unfortunately, the accuracy of a divided difference depends heavily on the stepsize (or stepsizes). Consequently, to reliably verify Adifor derivatives, we seek a method for selecting stepsizes that minimize the error. This technical report gives a method for selecting a stepsize that will give the 'most accurate' finite difference approximation. The mathematics that justifies the selction method is based on the work of Gill, Murray, and Wright.
July 2005
Solution-Recovery in L1-norm for Non-square Linear Systems: Deterministic Conditions and Open Questions
Yin Zhang
Consider an over-determined linear system A'x = b and an under-determined linear system By = c. Given b = A'x* + h, under what conditions x* will minimize the residual A'x - b in L1-norm? On the other hand, given c = Bh, under what conditions h will be the minimum L1-norm solutio to By = c? These two ``solution-recovery'' problems have been the focus of a number of recent works. Moreover, these two problems are equivalent under appropriate conditions on the data sets (A,b) and (B,c). In this paper, we give deterministic conditions for these solution-recoverary problems and raise a few open questuions. Some of the results in this paper are already known or partially known, but our derivations are different and thus may provide different perspectives.
A Sparse, Bound-Respecting Parametrization of Velocity Models
Eric Dussaud & William W. Symes
We present a parsimonious representation of velocity models which allows for user-defined placement of nodes. Mild restrictions are imposed on the data structure so that a computationally efficient algorithm can be used to smoothly approximate nodal values on finely sampled regular grids. The building block of the algorithm operates on one-dimensional arrays, allows for user-defined control of smoothness and guarantees that bounds are preserved. The specific data structure allows to carry out this process recursively to obtain multi-dimension smooth and regularly gridded velocity models.
April 2005
Kinematics of Shot-Geophone Migration
Christiaan C. Stolk, Maarten V. de Hoop, & William W. Symes
Prestack migration methods based on data binning produce {\em kinematic artifacts}, i.e. coherent events not corresponding to actual reflectors, in the prestack image volume. Shot-geophone migration, on the other hand, generally does not produce such artifacts when events to be migrated arrive in the data along non-turning rays. This condition is required for successful implementation via wavefield depth extrapolation (``survey sinking''). In contrast to prestack migration methods based on data binning, common image gathers produced by shot-geophone migration exhibit the appropriate semblance property in either offset domain (focussing at zero offset) or angle domain (focussing at zero slope), when the migration velocity is kinematically correct. Thus shot-geophone migration may be a particularly appropriate tool for migration velocity analysis of data exhibiting structural complexity.
April 2005
A Spatial Domain Decomposition Method for Parabolic Optimal Control Problems
M. Heinkenschloss & M. Herty
We present a non-overlapping spatial domain decomposition method for the solution of linear-quadratic parabolic optimal control problems. The spatial domain is decomposed into non-overlapping subdomains. The original parabolic optimal control problem is decomposed into smaller problems posed on space-time cylinder subdomains with auxiliary state and adjoint variables imposed as Dirichlet boundary conditions on the space-time interface boundary. The subdomain problems are coupled through Robin transmission conditions. This leads to a Schur complement equation in which the unknowns are the auxiliary state adjoint variables on the space-time interface boundary. The Schur complement operator is the sum of space-time subdomain Schur complement operators. The application of these subdomain Schur complement operators is equivalent to the solution of an subdomain parabolic optimal control problem. The subdomain Schur complement operators are shown to be invertible and the application of their inverses is equivalent to the solution of a related subdomain parabolic optimal control problem. We introduce a new family of Neumann-Neumann type preconditioners for the Schur complement system including several different coarse grid corrections. We compare the numerical performance of our preconditioners with an alternative approach recently introduced by Benamou.
Key words: Optimal control, parabolic equations, domain decomposition, Neumann-Neumann methods.
March 2005
Accelerating the Lee-Seung Algorithm for Nonnegative Matrix Factorization
Edward F. Gonzalez & Yin Zhang
Approximate nonnegative matrix factorization is an emerging technique with a wide spectrum of potential applications in data analysis. Currently, the most-used algorithms for this problem are those proposed by Lee and Seung. In this paper we present a variation of one of the Lee-Seung algorithms with a notably improved performance. We also show that algorithms of this type do not necessarily converge to local minima.
March 2005
A Symmetry Preserving Singular Value Decomposition
M. Shah & D.C. Sorensen
A reduced order representation of a large data set is often realized through a principle component analysis based upon a singular value decomposition (SVD) of the data. The left singular vectors of a truncated SVD provide the reduced basis. In several applications such as facial analysis and protein dynamics, structural symmetry is inherent in the data. Typically, reflective or rotational symmetry is expected to be present in these applications. In protein dynamics, determining this symmetry allows one to provide SVD major modes of motion that best describe the symmetric movements of the protein. In face detection, symmetry in the SVD allows for more efficient compression algorithms. Here, we present a method to compute the plane of reflective symmetry or the axis of rotational symmetry of a large set of points. Moreover, we develop a symmetry preserving singular value decomposition (SPSVD) that best approximates the given set while respecting the symmetry.
Interesting subproblems arise in the presence of noisy data or in situations where most, but not all, of the structure is symmetric. An important part of the determination of the axis of rotational symmetry or the plane of reflective symmetry is an iterative re-weighting scheme. This scheme is rapidly convergent in practice and seems to be very effective in ignoring outliers (points that do not respect the symmetry).
January 2005