This thesis studies how the Douglas–Rachford splitting technique can be leveraged for scalable computational optimal transport (OT). By carefully splitting the problem, we derive an algorithm with several advantages. First, the algorithm enjoys global convergence rates comparable to the state-of-the-art while benefiting from accelerated local rates. In contrast to other methods, it does not depend on hyperparameters that can cause numerical instability. This feature is particularly advantageous when low-precision floating points are used or if the data is noisy. Moreover, the updates can efficiently be carried out on GPUs and, therefore, benefit from the high degree of parallelization achieved via GPU computations. Furthermore, we show that the algorithm can be extended to handle a broad family of regularizers and constraints while enjoying the same theoretical and numerical properties. These factors combined result in a fast algorithm that can be applied to large-scale OT problems and regularized versions thereof, which we illustrate in several numerical experiments. In the first part of the main body of the thesis, we present how Douglas-Rachford splitting can be adapted for the unregularized OT problem to derive a fast algorithm. We present two global convergence guarantees for the resulting algorithm: a 1/k-ergodic rate and a linear rate. We also show that the stopping criteria of the algorithm can be computed on the fly with virtually no extra costs. Further, we specify how a GPU kernel can be efficiently implemented to carry out the operations needed for the algorithm. To show that the algorithm is fast, accurate, and robust, we run a series of numerical benchmarks that demonstrate the advantages of our algorithm. We then extend the algorithm to handle regularized OT using sparsity-promoting regularizers. The generalized algorithm will enjoy the same sublinear rate derived for the unregularized formulation. We also complement the global rate with local guarantees, establishing that, under non-degeneracy assumptions on the solution, the algorithm will identify the correct sparsity pattern of the solution in finitely many iterations. When the sparsity pattern is identified, a faster linear rate typically dominates. We also specify how to extend to the GPU implementation and the stopping criteria to handle regularized OT, and we subsequently specify how to backpropagate through the solver. We end this part of the thesis by presenting some numerical results, including performance on quadratically regularized OT and group Lasso regularized OT for domain adaptation, showing a substantial improvement compared to the state-of-the-art. In the last part of the thesis, we provide a more detailed analysis of the local behavior of the algorithm when applied to unregularized OT and quadratically regularized OT. We subsequently outline how to extend this analysis to several other sparsity-promoting regularizers. In the former two cases, we show that the update that constitutes the algorithm converges to a linear operator in finitely many iterations. By analyzing the spectral properties of these linear operators, we gain insights into the local behavior of the algorithm, and specifically, these results suggest how to tune stepsizes to obtain better local rates. / Denna avhandling behandlar hur Douglas–Rachford-splittning kan tillämpas för skalbara beräkningar av optimal transport (OT). Genom en noggrann splittning av problemet härleder vi en algoritm med flera fördelar. För det första åtnjuter algoritmen en global konvergenshastighet som är jämförbara med populära OT-lösare, samtidigt som den drar nytta av accelererade lokalahastigheter. Till skillnad från andra metoder är den inte beroende av hyperparametrar som kan orsaka numerisk instabilitet. Den här egenskapen är särskilt fördelaktig när lågprecisionsaritmetik används eller när data innehåller mycket brus. Uppdateringarna som algoritmen baseras på kan effektivt utföras på GPU:er och dra nytta av dess parallellberäkningar. Vi visar också att algoritmen kan utökas för att hantera en rad regulariseriseringar och bivillkor samtidigt som den åtnjuter liknande teoretiska och numeriska egenskaper. Tillsammans resulterar dessa faktorer i en snabb algoritm som kan tillämpas på storskaliga OT-problem samt flera av dess regulariserade varianter, vilket vi visar i flera numeriska experiment. I den första delen av avhandlingen presenterar vi hur Douglas-Rachford-splittning kan anpassas för det oregulariserade OT-problemet för att härleda en snabb algoritm. För den resulterande algoritmen presenterar vi två globala konvergensgarantier: en 1/k-ergodisk och en linjär konvergenshastighet. Vi presenterar också hur stoppkriterierna för algoritmen kan beräknas utan vidare extra kostnader. Dessutom specificerar vi hur en GPU-kärna kan implementeras för att effektivt utföra de operationer som algoritmen baseras på. För att visa att algoritmen är snabb, exakt och robust utför vi ett flertal numeriska experiment som påvisar flera fördelar över jämförbara algoritmer. Därefter utökar vi algoritmen för att hantera regulariserad OT med s.k. sparsity-promoting regularizers. Den generaliserade algoritmen åtnjuter samma sublinjära konvergenshastighet som härleddes för den oregulariserade fallet. Vi kompletterar också garantierna genom att tillhandahålla lokala garantier genom att fastställa att givet svaga antaganden om lösningen kommer algoritmen att identifiera den korrekta lösningens gleshetsstuktur inom ett begränsat antal iterationer. När glesheten identifierats dominerar typiskt sett en snabbare linjär konvergenshastighet. Vi specificerar också hur man utökar till GPU-kärnan och resultaten av stoppkriterierna för att hantera regulariserad OT, och vi specificerar sedan hur man differentiera genom lösaren. Vi avslutar den här delen av avhandlingen genom att presentera några numeriska resultat för kvadratiskt reglerad OT och group Lasso-reglerad OT för domänanpassning, vilket visar en betydande förbättring jämfört med de mest populära metoderna för dessa tillämpningar. I den sista delen av avhandlingen ger vi en mer detaljerad analys av algoritmens lokala beteende när den tillämpas på oregulerisad OT och kvadratiskt reglerad OT. Vi föreslår även sätt att studera flera andra fallet. I de två första fallen visar vi att uppdateringen som utgör algoritmen konvergerar till en linjär operator inom ett begränsat antal iterationer. Genom att analysera de spektrala egenskaperna hos dessa linjära operatorer får vi ytterligare insikter i algoritmens lokala beteende, och specifikt indikerar dessa resultat hur man kan justera steglängden för att uppnå ännu bättre konvergenshastigheter. / <p>QC 20231123</p>

Operator splitting methods have been applied to nonlinear partial differential equations that involve operators of different nature. The main idea of these methods is to decompose a complex equation into simpler sub-equations, which can be solved separately. The main advantage of the operator splitting methods is that they provide a great flexibility in choosing different numerical methods, depending on the feature of each sub-problem. In this dissertation, we have developed highly accurate and efficient numerical methods for several nonlinear partial differential equations, which involve both linear and nonlinear operators. We first propose a fast explicit operator splitting method for the modified Buckley-Leverett equations which include a third-order mixed derivatives term resulting from the dynamic effects in the pressure difference between the two phases. The method splits the original equation into two equations, one with a nonlinear convective term and the other one with high-order linear terms so that appropriate numerical methods can be applied to each of the split equations: The high-order linear equation is numerically solved using a pseudo-spectral method, while the nonlinear convective equation is integrated using the Godunov-type central-upwind scheme. The spatial order of the central-upwind scheme depends on the order of the piecewise polynomial reconstruction: We test both the second-order minmod-based reconstruction and fifth-order WENO5 one to demonstrate that using higher-order spatial reconstruction leads to more accurate approximation of solutions. We then propose fast and stable explicit operator splitting methods for two phase-field models (the molecular beam epitaxy equation with slope selection and the Cahn-Hilliard equation), numerical simulations of which require long time computations. The equations are split into nonlinear and linear parts. The nonlinear part is solved using a method of lines combined with an efficient large stability domain explicit ODE solver. The linear part is solved by a pseudo-spectral method, which is based on the exact solution and thus has no stability restriction on the time step size. We have verified the numerical accuracy of the proposed methods and demonstrated their performance on extensive one- and two-dimensional numerical examples, where different solution profiles can be clearly observed and are consistent with previous analytical studies. / Zhuolin Qu

Time splitting is a numerical procedure used in solution of partial differential equations whose solutions allow multiple time scales. Numerical schemes are split for handling the stiffness in equations, i.e. when there are multiple time scales with a few time scales being smaller than the others. When there are such terms with smaller time scales, due to the Courant number restriction, the computational cost becomes high if these terms are treated explicitly. In the present work a nonlinear advective equation is solved numerically using different techniques based on a generalised framework for splitting methods. The nonlinear advective equation was chosen because it has an analytical solution making comparisons with numerical schemes amenable and also because its nonlinearity mimics the equations encountered in atmospheric modelling. Using the nonlinear advective equation as a test bed, an analysis of the splitting methods and their influence on the split solutions has been made. An understanding of influence of splitting schemes requires knowledge of behaviour of unsplit schemes beforehand. Hence a study on unsplit methods has also been made. In the present work, using the nonlinear advective equation, it shown that the three time level schemes have high phase errors and underestimate energy (even though they have a higher order of accuracy in time). It is also found that the leap-frog method, which is used widely in atmospheric modelling, is the worst among examined unsplit methods. The semi implicit method, again a popular splitting method with atmospheric modellers is the worst among examined split methods. Three time-level schemes also need explicit filtering to remove the computational mode. This filtering can have a significant impact on the obtained numerical solutions, and hence three-time level schemes appear to be unattractive in the context of the nonlinear convective equation. Based on this experience, splitting methods for the two-time level schemes is proposed. These schemes realistically capture the phase and energy of the nonlinear advective equation.

There are some nonlinear models for pricing financial derivatives which can improve the linear Black-Scholes model introduced by Black, Scholes and Merton. In these models volatility is not constant anymore, but depends on some extra variables. It can be, for example, transaction costs, a risk from a portfolio, preferences of a large trader, etc. In this thesis we focus on these models. In the first chapter we introduce some important theory of financial derivatives. The second chapter is devoted to the volatility models. We derive three models concerning transaction costs (RAPM, Leland's  and Barles-Soner's model) and Frey's model which assumes a large (dominant) trader on the market. In the third and in the forth chapter we derive portfolio and make numerical experiments with a free boundary. We use the first order additive and the second order Strang splitting methods. We also use approximations of Barles-Soner's model using the identity function and introduce an approximation with the logarithm function of Barles-Soner's model. These models we finally compare with models where the volatility includes constant transaction costs.

Convex optimization problems are a class of mathematical problems which arise in numerous applications. Although interior-point methods can in principle solve these problems efficiently, they may become intractable for solving large-scale problems or be unsuitable for real-time embedded applications. Iterations of operator splitting methods are relatively simple and computationally inexpensive, which makes them suitable for these applications. However, some of their known limitations are slow asymptotic convergence, sensitivity to ill-conditioning, and inability to detect infeasible problems. The aim of this thesis is to better understand operator splitting methods and to develop reliable software tools for convex optimization. The main analytical tool in our investigation of these methods is their characterization as the fixed-point iteration of a nonexpansive operator. The fixed-point theory of nonexpansive operators has been studied for several decades. By exploiting the properties of such an operator, it is possible to show that the alternating direction method of multipliers (ADMM) can detect infeasible problems. Although ADMM iterates diverge when the problem at hand is unsolvable, the differences between subsequent iterates converge to a constant vector which is also a certificate of primal and/or dual infeasibility. Reliable termination criteria for detecting infeasibility are proposed based on this result. Similar ideas are used to derive necessary and sufficient conditions for linear (geometric) convergence of an operator splitting method and a bound on the achievable convergence rate. The new bound turns out to be tight for the class of averaged operators. Next, the OSQP solver is presented. OSQP is a novel general-purpose solver for quadratic programs (QPs) based on ADMM. The solver is very robust, is able to detect infeasible problems, and has been extensively tested on many problem instances from a wide variety of application areas. Finally, operator splitting methods can also be effective in nonconvex optimization. The developed algorithm significantly outperforms a common approach based on convex relaxation of the original nonconvex problem.

We present new results on absolute stability for BDF-ADI (Backward differentiation formula Alternating Direction Implicit) methods applied to a linear advection and diffusion equations. Unconditional absolute stability of the BDF2-ADI method is proven for advection and diffusion separately, as well as to the BDF3-ADI method for the purely-diffusive case. Conditional absolute stability of the BDF4-ADI is also proven for the purely-diffusive case, and stability regions for BDF3-ADI and BDF4- ADI are given in terms of the PDE coefficients and numerical parameters. Lastly, numerical experiments are presented to support the theoretical results and conjectures. These experiments also suggest future work.

We study numerical methods for time-dependent partial differential equations describing wave propagation, primarily applied to problems in quantum dynamics governed by the time-dependent Schrödinger equation (TDSE). We consider both methods for spatial approximation and for time stepping. In most settings, numerical solution of the TDSE is more challenging than solving a hyperbolic wave equation. This is mainly because the dispersion relation of the TDSE makes it very sensitive to dispersion error, and infers a stringent time step restriction for standard explicit time stepping schemes. The TDSE is also often posed in high dimensions, where standard methods are intractable. The sensitivity to dispersion error makes spectral methods advantageous for the TDSE. We use spectral or pseudospectral methods in all except one of the included papers. In Paper III we improve and analyse the accuracy of the Fourier pseudospectral method applied to a problem with limited regularity, and in Paper V we construct a matrix-free spectral method for problems with non-trivial boundary conditions. Due to its stiffness, the TDSE is most often solved using exponential time integration. In this thesis we use exponential operator splitting and Krylov subspace methods. We rigorously prove convergence for force-gradient operator splitting methods in Paper IV. One way of making high-dimensional problems computationally tractable is low-rank approximation. In Paper VI we prove that a splitting method for dynamical low-rank approximation is robust to singular values in the approximation approaching zero, a situation which is difficult to handle since it implies strong curvature of the approximation space. / eSSENCE

Dans une première partie, nous nous intéressons aux systèmes dynamiques gradients gouvernés par des fonctions non lisses, mais aussi non convexes, satisfaisant l'inégalité de Kurdyka-Lojasiewicz. Après avoir obtenu quelques résultats préliminaires pour la dynamique de la plus grande pente continue, nous étudions un algorithme de descente général. Nous prouvons, sous une hypothèse de compacité, que tout suite générée par ce schéma général converge vers un point critique de la fonction. Nous obtenons aussi de nouveaux résultats sur la vitesse de convergence, tant pour les valeurs que pour les itérés. Ce schéma général couvre en particulier des versions parallélisées de la méthode forward-backward, autorisant une métrique variable et des erreurs relatives. Cela nous permet par exemple de proposer une version non convexe non lisse de l'algorithme Levenberg-Marquardt. Enfin, nous proposons quelques applications de ces algorithmes aux problèmes de faisabilité, et aux problèmes inverses. Dans une seconde partie, cette thèse développe une dynamique de descente associée à des problèmes d'optimisation vectoriels sous contrainte. Pour cela, nous adaptons la dynamique de la plus grande pente usuelle aux fonctions à valeurs dans un espace ordonné par un cône convexe fermé solide. Cette dynamique peut être vue comme l'analogue continu de nombreux algorithmes développés ces dernières années. Nous avons un intérêt particulier pour les problèmes de décision multi-objectifs, pour lesquels cette dynamique de descente fait décroitre toutes les fonctions objectif au cours du temps. Nous prouvons l'existence de trajectoires pour cette dynamique continue, ainsi que leur convergence vers des points faiblement efficients. Finalement, nous explorons une nouvelle dynamique inertielle pour les problèmes multi-objectif, avec l'ambition de développer des méthodes rapides convergeant vers des équilibres de Pareto. / In a first part, we focus on gradient dynamical systems governed by non-smooth but also non-convex functions, satisfying the so-called Kurdyka-Lojasiewicz inequality.After obtaining preliminary results for a continuous steepest descent dynamic, we study a general descent algorithm. We prove, under a compactness assumption, that any sequence generated by this general scheme converges to a critical point of the function.We also obtain new convergence rates both for the values and the iterates. The analysis covers alternating versions of the forward-backward method, with variable metric and relative errors. As an example, a non-smooth and non-convex version of the Levenberg-Marquardt algorithm is detailed.Applications to non-convex feasibility problems, and to sparse inverse problems are discussed.In a second part, the thesis explores descent dynamics associated to constrained vector optimization problems. For this, we adapt the classic steepest descent dynamic to functions with values in a vector space ordered by a solid closed convex cone. It can be seen as the continuous analogue of various descent algorithms developed in the last years.We have a particular interest for multi-objective decision problems, for which the dynamic make decrease all the objective functions along time.We prove the existence of trajectories for this continuous dynamic, and show their convergence to weak efficient points.Then, we explore an inertial dynamic for multi-objective problems, with the aim to provide fast methods converging to Pareto points.

La simulation d’écoulements incompressibles pose de nombreuses difficultés. Une première est la question de savoir comment traiter la contrainte d’incompressibilité et le couplage vitesse/pression afin d’obtenir une solution précise à moindre coût. Pour cela, nous nous intéressons en particulier à deux méthodes de time splitting : la correction de pression et la correction de vitesse. Une seconde difficulté porte sur des conditions limites de sortie. Nous nous intéressons ici à deux d’entre elles : la condition limite de traction et la condition limite d’Orlanski. Après avoir détaillé les difficultés pouvant apparaître lors de l’implémentation des méthodes de time-splitting, nous proposons une nouvelle implémentation de la condition limite de traction qui permet d’améliorer les ordres de convergence obtenus. Nous nous intéressons ensuite à la condition limite d’Orlanski qui nécessite une certaine vitesse d’advection C dans la direction normale à la limite dont nous proposons ici une nouvelle définition. Nos propositions sont confrontées à de multiples écoulements physiques afin de valider leurs comportements : l’écoulement en aval d’une marche descendante, l’écoulement au niveau d’une bifurcation,l’écoulement autour d’un obstacle et des écoulements de Poiseuille-Rayleigh-Bénard. / One of the understudied difficulties in the simulation of incompressible flows is how to treat the incompressibilityconstraint and the velocity/pressure coupling in order to obtain an accurate solution at low computationnalcost. In this context, we develop two methods: pressure-correction and velocity-correction. An anotherdifficulty is due to the boundary conditions. We study here two of them : the traction boundary condition andthe Orlanski boundary condition. After having developed the difficulties that appears when implementing timesplittingmethods, we propose a new way to enforce the traction boundary condition which improves the orderof convergence. Then we propose a new definition of the advective velocity C which is needed for the Orlanskiboundary condition. Our propositions are validated against multiple physical flows: flow over a backward facingstep, flow around a biffurcation, flow around an obstacle and several Poiseuille-Rayleigh-Bénard flows.

<p>In this work, we study a class of Markovian information theoretic optimization problems motivated by the recent interests in incorporating mutual information as performance metrics which gives evident success in representation learning, feature extraction and clustering problems. In particular, we focus on the information bottleneck (IB) and privacy funnel (PF) methods and their recent multi-view, multi-source generalizations that gain attention because the performance significantly improved with multi-view, multi-source data. Nonetheless, the generalized problems challenge existing IB and PF solves in terms of the complexity and their abilities to tackle large-scale data. </p> <p>To address this, we study both the IB and PF under a unified framework and propose solving it through splitting methods, including renowned algorithms such as alternating directional method of multiplier (ADMM), Peaceman-Rachford splitting (PRS) and Douglas-Rachford splitting (DRS) as special cases. Our convergence analysis and the locally linear rate of convergence results give rise to new splitting method based IB and PF solvers that can be easily generalized to multi-view IB, multi-source PF. We implement the proposed methods with gradient descent and empirically evaluate the new solvers in both synthetic and real-world datasets. Our numerical results demonstrate improved performance over the state-of-the-art approach with significant reduction in complexity. Furthermore, we consider the practical scenario where there is distribution mismatch between training and testing data generating processes under a known bounded divergence constraint. In analyzing the generalization error, we develop new techniques inspired by the input-output mutual information approach and tighten the existing generalization error bounds.</p>