Higher order derivatives and norms of certain matrix functions

Tese de doutoramento, Matemática (Álgebra Lógica e Fundamentos), Universidade de Lisboa, Faculdade de Ciências, 2014

Krylov subspace methods for approximating functions of symmetric positive definite matrices with applications to applied statistics and anomalous diffusion

Matrix function approximation is a current focus of worldwide interest and finds application in a variety of areas of applied mathematics and statistics. In this thesis we focus on the approximation of A^(-α/2)b, where A ∈ ℝ^(n×n) is a large, sparse symmetric positive definite matrix and b ∈ ℝ^n is a vector. In particular, we will focus on matrix function techniques for sampling from Gaussian Markov random fields in applied statistics and the solution of fractional-in-space partial differential equations. Gaussian Markov random fields (GMRFs) are multivariate normal random variables characterised by a sparse precision (inverse covariance) matrix. GMRFs are popular models in computational spatial statistics as the sparse structure can be exploited, typically through the use of the sparse Cholesky decomposition, to construct fast sampling methods. It is well known, however, that for sufficiently large problems, iterative methods for solving linear systems outperform direct methods. Fractional-in-space partial differential equations arise in models of processes undergoing anomalous diffusion. Unfortunately, as the fractional Laplacian is a non-local operator, numerical methods based on the direct discretisation of these equations typically requires the solution of dense linear systems, which is impractical for fine discretisations. In this thesis, novel applications of Krylov subspace approximations to matrix functions for both of these problems are investigated. Matrix functions arise when sampling from a GMRF by noting that the Cholesky decomposition A = LL^T is, essentially, a `square root' of the precision matrix A. Therefore, we can replace the usual sampling method, which forms x = L^(-T)z, with x = A^(-1/2)z, where z is a vector of independent and identically distributed standard normal random variables. Similarly, the matrix transfer technique can be used to build solutions to the fractional Poisson equation of the form ϕn = A^(-α/2)b, where A is the finite difference approximation to the Laplacian. Hence both applications require the approximation of f(A)b, where f(t) = t^(-α/2) and A is sparse. In this thesis we will compare the Lanczos approximation, the shift-and-invert Lanczos approximation, the extended Krylov subspace method, rational approximations and the restarted Lanczos approximation for approximating matrix functions of this form. A number of new and novel results are presented in this thesis. Firstly, we prove the convergence of the matrix transfer technique for the solution of the fractional Poisson equation and we give conditions by which the finite difference discretisation can be replaced by other methods for discretising the Laplacian. We then investigate a number of methods for approximating matrix functions of the form A^(-α/2)b and investigate stopping criteria for these methods. In particular, we derive a new method for restarting the Lanczos approximation to f(A)b. We then apply these techniques to the problem of sampling from a GMRF and construct a full suite of methods for sampling conditioned on linear constraints and approximating the likelihood. Finally, we consider the problem of sampling from a generalised Matern random field, which combines our techniques for solving fractional-in-space partial differential equations with our method for sampling from GMRFs.

GPU accelerated algorithms for computing matrix function vector products with applications to exponential integrators and fractional diffusion

The efficient computation of matrix function vector products has become an important area of research in recent times, driven in particular by two important applications: the numerical solution of fractional partial differential equations and the integration of large systems of ordinary differential equations. In this work we consider a problem that combines these two applications, in the form of a numerical solution algorithm for fractional reaction diffusion equations that after spatial discretisation, is advanced in time using the exponential Euler method. We focus on the efficient implementation of the algorithm on Graphics Processing Units (GPU), as we wish to make use of the increased computational power available with this hardware. We compute the matrix function vector products using the contour integration method in [N. Hale, N. Higham, and L. Trefethen. Computing Aα, log(A), and related matrix functions by contour integrals. SIAM J. Numer. Anal., 46(5):2505–2523, 2008]. Multiple levels of preconditioning are applied to reduce the GPU memory footprint and to further accelerate convergence. We also derive an error bound for the convergence of the contour integral method that allows us to pre-determine the appropriate number of quadrature points. Results are presented that demonstrate the effectiveness of the method for large two-dimensional problems, showing a speedup of more than an order of magnitude compared to a CPU-only implementation.

Novel analytical and numerical methods for solving fractional dynamical systems

During the past three decades, the subject of fractional calculus (that is, calculus of integrals and derivatives of arbitrary order) has gained considerable popularity and importance, mainly due to its demonstrated applications in numerous diverse and widespread fields in science and engineering. For example, fractional calculus has been successfully applied to problems in system biology, physics, chemistry and biochemistry, hydrology, medicine, and finance. In many cases these new fractional-order models are more adequate than the previously used integer-order models, because fractional derivatives and integrals enable the description of the memory and hereditary properties inherent in various materials and processes that are governed by anomalous diffusion. Hence, there is a growing need to find the solution behaviour of these fractional differential equations. However, the analytic solutions of most fractional differential equations generally cannot be obtained. As a consequence, approximate and numerical techniques are playing an important role in identifying the solution behaviour of such fractional equations and exploring their applications. The main objective of this thesis is to develop new effective numerical methods and supporting analysis, based on the finite difference and finite element methods, for solving time, space and time-space fractional dynamical systems involving fractional derivatives in one and two spatial dimensions. A series of five published papers and one manuscript in preparation will be presented on the solution of the space fractional diffusion equation, space fractional advectiondispersion equation, time and space fractional diffusion equation, time and space fractional Fokker-Planck equation with a linear or non-linear source term, and fractional cable equation involving two time fractional derivatives, respectively. One important contribution of this thesis is the demonstration of how to choose different approximation techniques for different fractional derivatives. Special attention has been paid to the Riesz space fractional derivative, due to its important application in the field of groundwater flow, system biology and finance. We present three numerical methods to approximate the Riesz space fractional derivative, namely the L1/ L2-approximation method, the standard/shifted Gr¨unwald method, and the matrix transform method (MTM). The first two methods are based on the finite difference method, while the MTM allows discretisation in space using either the finite difference or finite element methods. Furthermore, we prove the equivalence of the Riesz fractional derivative and the fractional Laplacian operator under homogeneous Dirichlet boundary conditions – a result that had not previously been established. This result justifies the aforementioned use of the MTM to approximate the Riesz fractional derivative. After spatial discretisation, the time-space fractional partial differential equation is transformed into a system of fractional-in-time differential equations. We then investigate numerical methods to handle time fractional derivatives, be they Caputo type or Riemann-Liouville type. This leads to new methods utilising either finite difference strategies or the Laplace transform method for advancing the solution in time. The stability and convergence of our proposed numerical methods are also investigated. Numerical experiments are carried out in support of our theoretical analysis. We also emphasise that the numerical methods we develop are applicable for many other types of fractional partial differential equations.

Exponential integrators and a dual-scale model for wood drying

For the timber industry, the ability to simulate the drying of wood is invaluable for manufacturing high quality wood products. Mathematically, however, modelling the drying of a wet porous material, such as wood, is a diffcult task due to its heterogeneous and anisotropic nature, and the complex geometry of the underlying pore structure. The well{ developed macroscopic modelling approach involves writing down classical conservation equations at a length scale where physical quantities (e.g., porosity) can be interpreted as averaged values over a small volume (typically containing hundreds or thousands of pores). This averaging procedure produces balance equations that resemble those of a continuum with the exception that effective coeffcients appear in their deffnitions. Exponential integrators are numerical schemes for initial value problems involving a system of ordinary differential equations. These methods differ from popular Newton{Krylov implicit methods (i.e., those based on the backward differentiation formulae (BDF)) in that they do not require the solution of a system of nonlinear equations at each time step but rather they require computation of matrix{vector products involving the exponential of the Jacobian matrix. Although originally appearing in the 1960s, exponential integrators have recently experienced a resurgence in interest due to a greater undertaking of research in Krylov subspace methods for matrix function approximation. One of the simplest examples of an exponential integrator is the exponential Euler method (EEM), which requires, at each time step, approximation of φ(A)b, where φ(z) = (ez - 1)/z, A E Rnxn and b E Rn. For drying in porous media, the most comprehensive macroscopic formulation is TransPore [Perre and Turner, Chem. Eng. J., 86: 117-131, 2002], which features three coupled, nonlinear partial differential equations. The focus of the first part of this thesis is the use of the exponential Euler method (EEM) for performing the time integration of the macroscopic set of equations featured in TransPore. In particular, a new variable{ stepsize algorithm for EEM is presented within a Krylov subspace framework, which allows control of the error during the integration process. The performance of the new algorithm highlights the great potential of exponential integrators not only for drying applications but across all disciplines of transport phenomena. For example, when applied to well{ known benchmark problems involving single{phase liquid ow in heterogeneous soils, the proposed algorithm requires half the number of function evaluations than that required for an equivalent (sophisticated) Newton{Krylov BDF implementation. Furthermore for all drying configurations tested, the new algorithm always produces, in less computational time, a solution of higher accuracy than the existing backward Euler module featured in TransPore. Some new results relating to Krylov subspace approximation of '(A)b are also developed in this thesis. Most notably, an alternative derivation of the approximation error estimate of Hochbruck, Lubich and Selhofer [SIAM J. Sci. Comput., 19(5): 1552{1574, 1998] is provided, which reveals why it performs well in the error control procedure. Two of the main drawbacks of the macroscopic approach outlined above include the effective coefficients must be supplied to the model, and it fails for some drying configurations, where typical dual{scale mechanisms occur. In the second part of this thesis, a new dual{scale approach for simulating wood drying is proposed that couples the porous medium (macroscale) with the underlying pore structure (microscale). The proposed model is applied to the convective drying of softwood at low temperatures and is valid in the so{called hygroscopic range, where hygroscopically held liquid water is present in the solid phase and water exits only as vapour in the pores. Coupling between scales is achieved by imposing the macroscopic gradient on the microscopic field using suitably defined periodic boundary conditions, which allows the macroscopic ux to be defined as an average of the microscopic ux over the unit cell. This formulation provides a first step for moving from the macroscopic formulation featured in TransPore to a comprehensive dual{scale formulation capable of addressing any drying configuration. Simulation results reported for a sample of spruce highlight the potential and flexibility of the new dual{scale approach. In particular, for a given unit cell configuration it is not necessary to supply the effective coefficients prior to each simulation.

Geometric elasticity for graphics, simulation, and computation

We develop new algorithms which combine the rigorous theory of mathematical elasticity with the geometric underpinnings and computational attractiveness of modern tools in geometry processing. We develop a simple elastic energy based on the Biot strain measure, which improves on state-of-the-art methods in geometry processing. We use this energy within a constrained optimization problem to, for the first time, provide surface parameterization tools which guarantee injectivity and bounded distortion, are user-directable, and which scale to large meshes. With the help of some new generalizations in the computation of matrix functions and their derivative, we extend our methods to a large class of hyperelastic stored energy functions quadratic in piecewise analytic strain measures, including the Hencky (logarithmic) strain, opening up a wide range of possibilities for robust and efficient nonlinear elastic simulation and geometry processing by elastic analogy.

Measuring the significance of inconsistency in the Viewpoints framework

Measuring inconsistency is crucial to effective inconsistency management in software development. A complete measurement of inconsistency should focus on not only the degree but also the significance of inconsistency. However, most of the approaches available only take the degree of inconsistency into account. The significance of inconsistency has not yet been given much needed consideration. This paper presents an approach for measuring the significance of inconsistency arising from different viewpoints in the Viewpoints framework. We call an individual set of requirements belonging to different viewpoints a combined requirements collection in this paper. We argue that the
significance of inconsistency arising in a combined requirements collection is closely associated with global priority levels of requirements involved in the inconsistency. Here we assume that the global priority level of an individual requirement captures the relative importance of every viewpoint including this requirement as well as the local priority level of the requirement within the viewpoint. Then we use the synthesis of global priority levels of all the requirements in a combined collection to measure the significance of the
collection. Following this, we present a scoring matrix function to measure the significance of inconsistency in an inconsistent combined requirements collection, which describes the contribution made by each subset of the requirements collection to the significance of the set of requirements involved in the inconsistency. An ordering relationship between inconsistencies of two combined requirements collections, termed more significant than, is also presented by comparing their significance scoring matrix functions. Finally, these techniques were implemented in a prototype tool called IncMeasurer, which we developed as a proof of concept.

QRT: A Qr-based tridiagonalization algorithm for nonsymmetric matrices

The stable similarity reduction of a nonsymmetric square matrix to tridiagonal form has been a long-standing problem in numerical linear algebra. The biorthogonal Lanczos process is in principle a candidate method for this task, but in practice it is confined to sparse matrices and is restarted periodically because roundoff errors affect its three-term recurrence scheme and degrade the biorthogonality after a few steps. This adds to its vulnerability to serious breakdowns or near-breakdowns, the handling of which involves recovery strategies such as the look-ahead technique, which needs a careful implementation to produce a block-tridiagonal form with unpredictable block sizes. Other candidate methods, geared generally towards full matrices, rely on elementary similarity transformations that are prone to numerical instabilities. Such concomitant difficulties have hampered finding a satisfactory solution to the problem for either sparse or full matrices. This study focuses primarily on full matrices. After outlining earlier tridiagonalization algorithms from within a general framework, we present a new elimination technique combining orthogonal similarity transformations that are stable. We also discuss heuristics to circumvent breakdowns. Applications of this study include eigenvalue calculation and the approximation of matrix functions.

Latent TGF-beta Binding Proteins : Adhesive functions and matrix association of LTBP-2 and potential functions of LTBP-1 and LTBP-3 in mesothelioma

Latent transforming growth factor-beta (TGF-beta) binding proteins (LTBPs) -1, -3 and -4 are ECM components whose major function is to augment the secretion and matrix targeting of TGF-beta, a multipotent cytokine. LTBP-2 does not bind small latent TGF-beta but has suggested functions as a structural protein in ECM microfibrils. In the current work we focused on analyzing possible adhesive functions of LTBP-2 as well as on characterizing the kinetics and regulation of LTBP-2 secretion and ECM deposition. We also explored the role of TGF-beta binding LTBPs in endothelial cells activated to mimic angiogenesis as well as in malignant mesothelioma. We found that, unlike most adherent cells, several melanoma cell lines efficiently adhered to purified recombinant LTBP-2. Further characterization revealed that the adhesion was mediated by alpha3beta1 and alpha6beta1 integrins. Heparin also inhibited the melanoma cell adhesion suggesting a role for heparan sulphate proteoglycans. LTBP-2 was also identified as a haptotactic substrate for melanoma cell migration. We used cultured human embryonic lung fibroblasts to analyze the temporal and spatial association of LTBP-2 into ECM. By We found that LTBP-2 was efficiently assembled to the ECM only in confluent cultures following the deposition of fibronectin (FN) and fibrillin-1. In early, subconfluent cultures it remained primarily in soluble form after secretion. LTBP-2 colocalized transiently with FN and fibrillin-1. Silencing of fibrillin-1 expression by lentiviral shRNAs profoundly disrupted the deposition of LTBP-2 indicating that the ECM association of LTBP-2 depends on a pre-formed fibrillin-1 network. Considering the established role of TGF-beta as a regulator of angiogenesis we induced morphological activation of endothelial cells by phorbol 12-myristate 13-acetate (PMA) and followed the fate of LTBP-1 in the endothelial ECM. This resulted in profound proteolytic processing of LTBP-1 and release of latent TGF-beta complexes from the ECM. The processing was coupled with increased activation of MT-MMPs and specific upregulation of MT1-MMP. The major role of MT1-MMP in the proteolysis of LTBP-1 was confirmed by suppressing the expression with lentivirally induced short-hairpin RNAs as well as by various metalloproteinases inhibitors. TGF-beta can promote tumorigenesis of malignant mesothelioma (MM), which is an aggressive tumor of the pleura with poor prognosis. TGF-beta activity was analyzed in a panel of MM tumors by immunohistochemical staining of phosphorylated Smad-2 (P-Smad2). The tumor cells were strongly positive for P-Smad2 whereas LTBP-1 immunoreactivity was abundant in the stroma, and there was a negative correlation between LTBP-1 and P-Smad2 staining. In addition, the high P-Smad2 immunoreactivity correlated with shorter survival of patients. mRNA analysis revealed that TGF-beta1 was the most highly expressed isoform in both normal human pleura and MM tissue. LTBP-1 and LTBP-3 were both abundantly expressed. LTBP-1 was the predominant isoform in established MM cell lines whereas the expression of LTBP-3 was high in control cells. Suppression of LTBP-3 expression by siRNAs resulted in increased TGF-beta activity in MM cell lines accompanied by decreased proliferation. Our results suggest that decreased expression of LTBP-3 in MM could alter the targeting of TGF-beta to the ECM and lead to its increased activation. The current work emphasizes the coordinated process of the assembly and appropriate targeting of LTBPs with distinct adhesive or cytokine harboring properties into the ECM. The hierarchical assembly may have implications in the modulation of signaling events during morphogenesis and tissue remodeling.

Structure and representation of correlation functions and the density matrix for a statistical wave field in optics

A systematic structure analysis of the correlation functions of statistical quantum optics is carried out. From a suitably defined auxiliary two‐point function we are able to identify the excited modes in the wave field. The relative simplicity of the higher order correlation functions emerge as a byproduct and the conditions under which these are made pure are derived. These results depend in a crucial manner on the notion of coherence indices and of unimodular coherence indices. A new class of approximate expressions for the density operator of a statistical wave field is worked out based on discrete characteristic sets. These are even more economical than the diagonal coherent state representations. An appreciation of the subtleties of quantum theory obtains. Certain implications for the physics of light beams are cited.

Putative functions of extracellular matrix glycoproteins in secondary palate morphogenesis

Cleft palate is a common birth defect in humans. Elevation and fusion of paired palatal shelves are coordinated by growth and transcription factors, and mutations in these can cause malformations. Among the effector genes for growth factor signaling are extracellular matrix (ECM) glycoproteins. These provide substrates for cell adhesion (e.g., fibronectin, tenascins), but also regulate growth factor availability (e.g., fibrillins). Cleft palate in Bmp7 null mouse embryos is caused by a delay in palatal shelf elevation. In contrast, palatal shelves of Tgf-β3 knockout mice elevate normally, but a cleft develops due to their failure to fuse. However, nothing is known about a possible functional interaction between specific ECM proteins and Tgf-β/Bmp family members in palatogenesis. To start addressing this question, we studied the mRNA and protein distribution of relevant ECM components during secondary palate development, and compared it to growth factor expression in wildtypewild type and mutant mice. We found that fibrillin-2 (but not fibrillin-1) mRNA appeared in the mesenchyme of elevated palatal shelves adjacent to the midline epithelial cells, which were positive for Tgf-β3 mRNA. Moreover, midline epithelial cells started expressing fibronectin upon contact of the two palatal shelves. These findings support the hypothesis that fibrillin-2 and fibronectin are involved in regulating the activity of Tgf-β3 at the fusing midline. In addition, we observed that tenascin-W (but not tenascin-C) was misexpressed in palatal shelves of Bmp7-deficient mouse embryos. In contrast to tenascin-C, tenascin-W secretion was strongly induced by Bmp7 in embryonic cranial fibroblasts in vitro. These results are consistent with a putative function for tenascin-W as a target of Bmp7 signaling during palate elevation. Our results indicate that distinct ECM proteins are important for morphogenesis of the secondary palate, both as downstream effectors and as regulators of Tgf-β/Bmp activity.