Eldad Haber

The goal of this thesis is to develop new numerical techniques for the solution ofnoisy optimization problems. Noisy optimization problems occur in diverse fields,such as engineering design, machine learning, operations research, and geophysics. Thus, effective numerical techniques for solving noisy optimization problems are of practical interest in a variety of industrial applications. This thesis presents three numerical techniques for solving noisy optimization problems: Secant penalized Broyden-Fletcher-Goldfarb-Shanno (SP-BFGS), Neural Network Accelerated Implicit Filtering (NNAIF), and Linear Model Based Gradient Estimation (LMBGE).SP-BFGS is a new variant of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method designed to perform well when gradient measurements are corrupted by noise. SP-BFGS provides a means of incrementally updating the new inverse Hessian approximation with a controlled amount of bias towards the previous inverse Hessian approximation, which allows one to replace the overwriting nature of the original BFGS update with an averaging nature that resists the destructive effects of noise and can cope with negative curvature measurements. I study the theoretical and empirical behaviour of SP-BFGS in this thesis.NNAIF intelligently combines the established literature on implicit filtering (IF)optimization methods with a neural network surrogate model of the objective function, resulting in accelerated derivative free methods for unconstrained optimization problems. I show that NNAIF directly inherits the convergence properties of IF optimization methods, and study the performance of NNAIF via numerical experiments with noisy problems from the Constrained and Unconstrained Testing Environment with safe threads (CUTEst) optimization benchmark set.LMBGE leverages techniques from linear inverse problems to estimate the gradient of a noisy objective function, and can be viewed as a noise robust extension of standard finite difference based approaches for gradient estimation. I demonstrate the use of LMBGE in the solution of the coverage directed generation (CDG) problem, an industrial problem encountered in the verification subfield of electrical and computer engineering.

View record

Machine learning is used daily in areas such as security, medical care, and financial systems. Failures in such institutions can have dire consequences. Adversarial attacks on deep neural networks exploit instabilities in the network with regard to noise and thus are important problems in the field of machine learning. In recent years, much attention has been given to the design of architectures and procedures that help mitigate the effect of such attacks. In this work, we take a different route and present a novel approach for dealing with adversarial attacks by soft-labeling data.In this work, we extend an existing popular image dataset: CIFAR-10. CIFAR-10 has disjoint, well-defined, and exhaustive labels with 10 classes. We design and execute a survey and graphical user interface to collect individual soft labels for 10,000 images. In the survey, the participants are shown one image at a time and asked to rate how much characteristics from each of the classes were in the image. These answers are then turned into soft labels. To our knowledge, our novel dataset is the only soft-labeled dataset with both unambiguous images and soft labels that are independent of the hard labels.Because soft labeled datasets are more complex than hard labeled datasets, using them will require larger datasets for training. For this, we present a novel algorithm that uses the alternating direction method of multipliers to solve multinomial logistic regression in parallel. The algorithm shows improved generalization. We train the residual neural network with our novel dataset and compare it to models trained with the original CIFAR-10 with and without uniform label smoothing regularization. We show that models trained with soft-labeled data resist targeted adversarial attacks better than these other models.

View record

The work described in this thesis examines transient geophysical electromagneticforward modelling and inversion in the presence of induced polarization (IP) effects. The thesis introduces a new method of modelling IP using stretched exponential relaxation. A three-dimensional (3D) forward modelling algorithm takingfull account of the coupling of IP and electromagnetic induction is developed. Thestretched exponential modelling algorithm has been implemented using efficientnumerical methods that allow it to tackle large-scale problems and are amenableto use in inversion. In particular, a parallel time-stepping technique has been de-veloped that allows transient electric fields at multiple time steps to be computedsimultaneously. The behavior of the stretched exponential model is demonstratedby applying it to synthetic numerical examples that simulate a grounded sourceIP survey with significant electromagnetic induction effects and a concentric-loopairborne electromagnetic sounding over a polarizable body.An inversion algorithm using the stretched exponential model was developedthat is able to recover the 3D structure of physical properties of the earth relatedto IP from transient geophysical electromagnetic data. The method is tested on asimple synthetic example problem. The thesis finishes with the development ofa novel stochastic parametric level-set inversion algorithm, which could be usefulin applying stretched exponential inversion to real world problems in the future.The algorithm addresses some of the shortcomings of the simple inversion approach used for stretched exponential inversion earlier in the thesis. The stochasticparametric inversion algorithm is used to solve shape reconstruction inverse problems in which the object of interest is embedded in a heterogeneous backgroundmedium that is known only approximately. Shape reconstruction is posed as a stochastic programming problem, in which the background medium is treated asa random field with a known probability distribution. It is demonstrated that byusing accelerated stochastic gradient descent the method can be applied to large-scale problems. The capabilities of the method are demonstrated on a simple 2Dmodel problem and in a more demanding application to a 3D inverse conductivityproblem in geophysical imaging.

View record

Inverse modeling is a powerful tool for extracting information about the subsurface from geophysical and hydrologic data. Geophysical inverse problems are inherently multidisciplinary, requiring elements from the relevant physics, numerical simulation, and optimization, as well as knowledge of the geologic setting, hydrologic processes, and a comprehension of the interplay between all of these elements. Increasingly geoscientists are tackling complex problems that require integration of multiple types of information in order to better characterize the subsurface. However, many of the sub-fields of geophysics are developing simulation and inversion approaches, algorithms, and supporting software in isolation. This isolation is a barrier to quantitative integration and leads to inefficiencies in advancing interdisciplinary research. Greater efficiencies, and higher quality outcomes, could be achieved if (hydro)geophysicists had a common framework to accelerate an integrated approach. The main goal of my thesis is to organize the components of (hydro)geophysical simulations and inverse problems, and synthesize these into a comprehensive, modular framework.The development of a geophysical framework requires considering a number of disciplines and geophysical problems (e.g. electromagnetics and potential fields) to ensure generality as well as extensibility. However, the goal is also to have the framework work outside of geophysics and most notably in hydrogeology; vadose zone fluid flow is used as a model problem. Fluid flow in the vadose zone is governed by the Richards equation; it is parameterized by hydraulic conductivity, which is a nonlinear function of pressure head. The computational scalability of the Richards equation inversion is a significant challenge for three dimensional inversions in hydrogeophysics. Existing work explicitly calculates the sensitivity matrix using finite difference or automatic differentiation, however, for large-scale problems these methods are constrained by computation and memory. This dissertation provides an implicit sensitivity algorithm that enables large-scale inversion problems for distributed parameters in the Richards equation to become tractable on modest computational resources.

View record

In the modern era of diminishing returns on fixed exploration budgets, challenging targets, and ever-increasing numbers of multi-parameter datasets, proper management and integration of available data is a crucial component of any mineral exploration program. Machine learning algorithms have successfully been used for years by the technology sector to accomplish just this task on their databases, and recent developments aim at appropriating these successes to the field of mineral exploration. Framing the exploration task as a supervised learning problem, the geological, geochemical and geophysical information can be used as training data, and known mineral occurences can be used as training labels. The goal is to parameterize the complex relationships between the data and the labels such that mineral potential can be estimated in under-explored regions using available geoscience data. Numerous models and algorithms have been attempted for mineral prospectivity mapping in the past, and in this thesis we propose two new approaches. The first is a modified support vector machine algorithm which incorporates uncertainties on both the data and the labels. Due to the nature of geoscience data and the characteristics of the mineral prospectivity mapping problem, uncertainties are known to be very important. The algorithm is demonstrated on a synthetic dataset to highlight this importance, and then used to generate a prospectivity map for copper-gold porphyry targets in central British Columbia using the QUEST dataset as a case study. The second approach, convolutional neural networks, was selected due to its inherent sensitivity to spatial patterns. Though neural networks have been used for mineral prospectivity mapping, convolutional neural nets have yet to be applied to the problem. Having gained extreme popularity in the computer vision field for tasks involving image segmentation, identification and anomaly detection, the algorithm is ideally suited to handle the mineral prospectivity mapping problem. A CNN code is developed in Julia, then tested on a synthetic example to illustrate its effectiveness at identifying coincident structures in a multi-modal dataset. Finally, a subset of the QUEST dataset is used to generate a prospectivity map using CNNs.

View record

Parameter and state estimation for groundwater models within a coupled hydrogeophysical framework has become common in the last few years as it has been shown that such estimates are usually better than those from a single data inversion. Different approaches have been suggested in literature to combine the essentially two different modalities in order to obtain better estimates for groundwater models, and improve monitoring of processes such as solute transport. However, the coupled approaches usually come at a price of higher computational cost and difficulties in coupling the geophysical and groundwater inverse problems. Unlike in other studies, we developed both the groundwater and geophysical models in the same computational environment in order to test different minimization strategies. When solving the coupled inverse problem, the objective function consists of data misfit and regularization terms as well as a coupling term that relates groundwater and geophysical states. We present a novel approach to solve the inverse problem using an Alternating Direction Method of Multipliers (ADMM) to minimize the coupled objective function. ADMM enables us to treat the groundwater and geophysical part separately and thus use existing software with minor changes.However, ADMM as well as many other coupled approaches relies on implementing some petrophysical relationship to couple the groundwater and geophysical variable. Such relationships are usually uncertain and hard to parametrize for a large region and can potentially produce solute mass errors in the final model estimates. Therefore, in this thesis we examine coupled approaches that replace the fixed petrophysical relationship by a more loose structure similarity constraint. Besides, we propose efficient computational methods to minimize the objective function when there is no explicit petrophysical constraint. All approaches were tested on 3D synthetic examples. In the solute tracer test we estimated hydraulic conductivity or solute distribution using a structure coupled inversion, and were able to reduce the errors compared to a single data inversion alone. For a more complex example of seawater intrusion we implemented the ADMM method, and obtained better estimates for the solute distribution compared to just considering each data separately, or solving the problem with a simple coupled approach.

View record

Airborne electromagnetic (AEM) data is commonly collected for detecting buried natural resources, and this technique is sensitive to subsurface electrical resistivity distributions. The subsequent process of 3D AEM inversion constructs a physical property model from this data in order to better understand the size and shape of potential hidden resources. This thesis is designed to develop practical strategies to improve 3D AEM inversion accuracy in geologic settings where AEM data sets produce inconsistent or unsatisfactory inversion results.In this research, two overarching problematic scenarios are examined. First, in regions where an AEM survey overlaps with other electromagnetic data sets, a novel cooperative approach is introduced. This method is first tested on synthetic data where instead of producing an inversion model from each data set, the cooperative algorithm finds one consistent 3D resistivity model with improved resolution. The approach is then applied to field data over a high-sulfidation epithermal gold deposit where similar improvements are observed.The second scenario relates to improving 3D AEM inversions over thin conductive anomalies, a common geophysical target for copper and gold deposits. A new parametric inversion is developed using skewed Gaussian ellipsoids to represent target bodies. The approach is general but applied to frequency and time-domain AEM data with one or multiple anomalies. Combined with a voxel algorithm, the parametric inversion forms a hybrid approach capable of resolving thin conductive targets with smooth surrounding features. This hybrid technique is tested on synthetic data over conductive targets in a resistive background, and consistently produces models that are easier to interpret compared to voxel inversions alone. Field examples from a volcanogenic massive sulfide and an orogenic gold deposit highlight the practical nature of this method to image conductive mineralization with increased precision.The thesis concludes by analyzing a setting where multiple spatially overlapping AEM data sets exist over narrow conductive anomalies. Here, parametric, cooperative and voxel inversions are applied together to generate a consistent 3D resistivity model with thin targets and smooth background features. This section includes a discussion about potential pitfalls of such an approach when incompatible field measurements are encountered.

View record

Accurate and efficient simulation of electromagnetic responses in realistic geophysical settings is crucial to the exploration, imaging, and characterization of buried natural resources, such as mineral and hydrocarbon deposits. However, in practice, these simulations are computationally expensive. The geophysical settings consider highly heterogeneous media and features at multiple spatial scales that require a very large mesh to be accurately represented. This results in a system of equations to be solved that often exceeds the limits of average computers. Thus, the key is to reduce the problem size but retain the accuracy of the electromagnetic responses. Upscaling and multiscale techniques have been successfully applied to the problem of simulating fluid flow through heterogeneous porous media, where they are able to drastically reduce the size of the resulting fine-mesh system by casting it into a coarse-mesh system that is much cheaper to solve, while achieving a level of accuracy similar to that obtained with conventional discretization schemes. Recognizing the success that such techniques have had in fluid flow applications, this dissertation extends their use for application to electromagnetic modeling. In this dissertation, two new parallel simulation methods for the quasi-static Maxwell’s equations in the frequency domain are proposed: an upscaling framework for the electrical conductivity, and a multiscale finite volume with oversampling method. Both methods are combined with an adaptive mesh refinement technique (OcTree) to boost their computational performance. The performance of these methods is demonstrated by using field-inspired and synthetic examples that include a large electrical conductivity contrast. This investigation shows that both proposed methods are feasible to tackle geophysical electromagnetic problems, where being able to reduce the size of the problem can be particularly advantageous when extended domains are considered or when the mesh must capture the spatial distribution of the media heterogeneity outside the region where the electromagnetic responses are measured. Furthermore, both methods are new contributions to the literature in the field of computational methods in geophysical electromagnetics. Finally, both methods increase the current predictive and analytic capabilities by making the simulation of electromagnetic responses in larger and more complex geophysical settings more feasible than currently is possible.

View record

This thesis examines the derivation and implementation of the discrete adjoint method for several time-stepping methods. Our results are important for gradient-based numerical optimization in the context of large-scale model calibration problems that are constrained by nonlinear time-dependent PDEs. To this end, we discuss finding the gradient and the action of the Hessian of the data misfit function with respect to three sets of parameters: model parameters, source parameters and the initial condition. We also discuss the closely related topic of computing the action of the sensitivity matrix on a vector, which is required when performing a sensitivity analysis. The gradient and Hessian of the data misfit function with respect to these parameters requires the derivatives of the misfit with respect to the simulated data, and we give the procedures for computing these derivatives for several data misfit functions that are of use in seismic imaging and elsewhere. The methods we consider can be divided into two categories, linear multistep (LM) methods and Runge-Kutta (RK) methods, and several variants of these are discussed. Regular LM and RK methods can be used for ODE systems arising from the semi-discretization of general nonlinear time-dependent PDEs, whereas implicit-explicit and staggered variants can be applied when the PDE has a more specialized form. Exponential time-differencing RK methods are also discussed. The implementation of the associated adjoint time-stepping methods is discussed in detail. Our motivation is the application of the discrete adjoint method to high-order time-stepping methods, but the approach taken here does not exclude lower-order methods. All of the algorithms have been implemented in MATLAB using an object-oriented design and are written with extensibility in mind. For exponential RK methods it is illustrated numerically that the adjoint methods have the same order of accuracy as their corresponding forward methods, and for linear PDEs we give a simple proof that this must always be the case. The applicability of some of the methods developed here to pattern formation problems is demonstrated using the Swift-Hohenberg model.

View record

Imaging and prediction of fluid flow within the subsurface provides information crucial to decision making processes in fields such as groundwater management and enhanced oil recovery. The flow of a fluid through a reservoir depends primarily on the permeability of the subsurface rock; a quantity that is often unknown throughout the entire domain of the reservoir. One method for predicting flow is to estimate the permeability of the reservoir and simulate flow through a mathematical subsurface flow model. Given the model, flow data can be inverted to estimate the permeability. However, this inversion approach can lead to inaccurate results due to the sparse sampling of flow data, and thus inaccurate predictions. To acquire a higher sampling of data, geophysical survey techniques are applied in order to efficiently collect a higher density of data sampled at the surface. These data are sensitive to changes to the geophysical properties of the reservoir due to flow. Inversion of geophysical data then provides images of changes to the geophysical properties of the reservoir. In order to estimate the flow parameters using geophysical data, the two mathematical models require coupling. The thesis therefore proposes two approaches to improve the imaging and prediction of flow. First, a novel coupled inverse problem for estimating the fluid velocity field and the initial geophysical property model from geophysical data is developed. Second, a new method of optimally designing the geophysical survey for the coupled inverse problem is developed. The new adaptive design approach builds on traditional A-Optimal design methods such that historic data are included in the design algorithm. This produces designs that adapt with flow in the subsurface and reduce the collection of unnecessary data. Both the coupled inverse problem and adaptive survey design method are demonstrated using a seismic tomography geophysical survey and a tracer advection fluid flow model. Numerical examples show that the coupled approach yields an improved flow estimate as well as improved image quality, while the adaptive optimal designs provide sufficient geophysical data.

View record

In this thesis, we introduce a new class of embarrassingly parallel parameter learning algorithms for Markov random fields (MRFs) with untied parameters, which are efficient for a large class of practical models. The algorithms parallelize naturally over cliques and, for graphs of bounded degree, have complexity that is linear in the number of cliques. We refer to these algorithms with the acronym LAP, which stands for Linear And Parallel. Unlike their competitors, the marginal versions of the proposed algorithms are fully parallel and for log-linear models they are also data efficient, requiring only the local sufficient statistics of the data to estimate parameters. LAP algorithms are ideal for parameter learning in big graphs and big data applications.The correctness of the newly proposed algorithms relies heavily on the existence and uniqueness of the normalized potential representation of an MRF. We capitalize on this theoretical result to develop a new theory of correctness and consistency of LAP estimators corresponding to different local graph neighbourhoods. This theory also establishes a general condition on composite likelihood decompositions of MRFs that guarantees the global consistency of distributed estimators, provided the local estimators are consistent. We introduce a conditional variant of LAP that enables us to attack parameter estimation of fully-observed models of arbitrary connectivity, including fully connected Boltzmann distributions. Once again, we show consistency for this distributed estimator, and relate it to distributed pseudo-likelihood estimators.Finally, for linear and non-linear inverse problems with a sparse forward operator, we present a new algorithm, named iLAP, which decomposes the inverse problem into a set of smaller dimensional inverse problems that can be solved independently. This parallel estimation strategy is also memory efficient.

View record