Bilinear systems are a class of nonlinear dynamical systems that arise in a variety of applications. In order to obtain a sufficiently accurate representation of the underlying physical phenomenon, these models frequently have state-spaces of very large dimension, resulting in the need for model reduction. In this work, we introduce two new methods for the model reduction of bilinear systems in an interpolation framework. Our first approach is to construct reduced models that satisfy multipoint interpolation constraints defined on the Volterra kernels of the full model. We show that this approach can be used to develop an asymptotically optimal solution to the H_2 model reduction problem for bilinear systems. In our second approach, we construct a solution to a bilinear system realization problem posed in terms of constructing a bilinear realization whose kth-order transfer functions satisfy interpolation conditions in k complex variables. The solution to this realization problem can be used to construct a bilinear system realization directly from sampling data on the kth-order transfer functions, without requiring the formation of the realization matrices for the full bilinear system. / Ph. D.

In science and engineering, many applications require the solution of a sequence of linear systems. There are many ways to solve linear systems and we always look for methods that are faster and/or require less storage. In this dissertation, we focus on solving these systems with Krylov subspace methods and how to obtain effective preconditioners inexpensively. We first present an application for electronic structure calculation. A sequence of slowly changing linear systems is produced in the simulation. The linear systems change by rank-one updates. Properties of the system matrix are analyzed. We use Krylov subspace methods to solve these linear systems. Krylov subspace methods need a preconditioner to be efficient and robust. This causes the problem of computing a sequence of preconditioners corresponding to the sequence of linear systems. We use recycling preconditioners, which is to update and reuse existing preconditioner. We investigate and analyze several preconditioners, such as ILU(0), ILUTP, domain decomposition preconditioners, and inexact matrix-vector products with inner-outer iterations. Recycling preconditioners produces cumulative updates to the preconditioner. To reduce the cost of applying the preconditioners, we propose approaches to truncate the cumulative preconditioner updates, which is a low-rank matrix. Two approaches are developed. The first one is to truncate the low-rank matrix using the best approximation given by the singular value decomposition (SVD). This is effective if many singular values are close to zero. If not, based on the ideas underlying GCROT and recycling, we use information from an Arnoldi recurrence to determine which directions to keep. We investigate and analyze their properties. We also prove that both truncation approaches work well under suitable conditions. We apply our truncation approaches on two applications. One is the Quantum Monte Carlo (QMC) method and the other is a nonlinear second order partial differential equation (PDE). For the QMC method, we test both truncation approaches and analyze their results. For the PDE problem, we discretize the equations with finite difference method, solve the nonlinear problem by Newton's method with a line-search, and utilize Krylov subspace methods to solve the linear system in every nonlinear iteration. The preconditioner is updated by Broyden-type rank-one updates, and we truncate the preconditioner updates by using the SVD finally. We demonstrate that the truncation is effective. In the last chapter, we develop a matrix reordering algorithm that improves the diagonal dominance of Slater matrices in the QMC method. If we reorder the entire Slater matrix, we call it global reordering and the cost is O(N^3), which is expensive. As the change is geometrically localized and impacts only one row and a modest number of columns, we propose a local reordering of a submatrix of the Slater matrix. The submatrix has small dimension, which is independent of the size of Slater matrix, and hence the local reordering has constant cost (with respect to the size of Slater matrix). / Ph. D.

Thesis (MSc)--Stellenbosch University, 2011. / ENGLISH ABSTRACT: In several numerical approaches to PDEs shifted linear systems of the form (zI - A)x = b, need to be solved for several values of the complex scalar z. Often, these linear systems are large and sparse. This thesis investigates efficient numerical methods for these systems that arise from a contour integral approximation to PDEs and compares these methods with direct solvers. In the first part, we present three model PDEs and discuss numerical approaches to solve them. We use the first problem to demonstrate computations with a dense matrix, the second problem to demonstrate computations with a sparse symmetric matrix and the third problem for a sparse but nonsymmetric matrix. To solve the model PDEs numerically we apply two space discrerization methods, namely the finite difference method and the Chebyshev collocation method. The contour integral method mentioned above is used to integrate with respect to the time variable. In the second part, we study a Hessenberg reduction method for solving shifted linear systems with a dense matrix and present numerical comparison of it with the built-in direct linear system solver in SciPy. Since both are direct methods, in the absence of roundoff errors, they give the same result. However, we find that the Hessenberg reduction method is more efficient in CPU-time than the direct solver. As application we solve a one-dimensional version of the heat equation. In the third part, we present efficient techniques for solving shifted systems with a sparse matrix by Krylov subspace methods. Because of their shift-invariance property, the Krylov methods allow one to obtain approximate solutions for all values of the parameter, by generating a single approximation space. Krylov methods applied to the linear systems are generally slowly convergent and hence preconditioning is necessary to improve the convergence. The use of shift-invert preconditioning is discussed and numerical comparisons with a direct sparse solver are presented. As an application we solve a two-dimensional version of the heat equation with and without a convection term. Our numerical experiments show that the preconditioned Krylov methods are efficient in both computational time and memory space as compared to the direct sparse solver. / AFRIKAANSE OPSOMMING: In verskeie numeriese metodes vir PDVs moet geskuifde lineêre stelsels van die vorm (zI − A)x = b, opgelos word vir verskeie waardes van die komplekse skalaar z. Hierdie stelsels is dikwels groot en yl. Hierdie tesis ondersoek numeriese metodes vir sulke stelsels wat voorkom in kontoerintegraalbenaderings vir PDVs en vergelyk hierdie metodes met direkte metodes vir oplossing. In die eerste gedeelte beskou ons drie model PDVs en bespreek numeriese benaderings om hulle op te los. Die eerste probleem word gebruik om berekenings met ’n vol matriks te demonstreer, die tweede probleem word gebruik om berekenings met yl, simmetriese matrikse te demonstreer en die derde probleem vir yl, onsimmetriese matrikse. Om die model PDVs numeries op te los beskou ons twee ruimte-diskretisasie metodes, naamlik die eindige-verskilmetode en die Chebyshev kollokasie-metode. Die kontoerintegraalmetode waarna hierbo verwys is word gebruik om met betrekking tot die tydveranderlike te integreer. In die tweede gedeelte bestudeer ons ’n Hessenberg ontbindingsmetode om geskuifde lineêre stelsels met ’n vol matriks op te los, en ons rapporteer numeriese vergelykings daarvan met die ingeboude direkte oplosser vir lineêre stelsels in SciPy. Aangesien beide metodes direk is lewer hulle dieselfde resultate in die afwesigheid van afrondingsfoute. Ons het egter bevind dat die Hessenberg ontbindingsmetode meer effektief is in terme van rekenaartyd in vergelyking met die direkte oplosser. As toepassing los ons ’n een-dimensionele weergawe van die hittevergelyking op. In die derde gedeelte beskou ons effektiewe tegnieke om geskuifde stelsels met ’n yl matriks op te los, met Krylov subruimte-metodes. As gevolg van hul skuifinvariansie eienskap, laat die Krylov metodes mens toe om benaderde oplossings te verkry vir alle waardes van die parameter, deur slegs een benaderingsruimte voort te bring. Krylov metodes toegepas op lineêre stelsels is in die algemeen stadig konvergerend, en gevolglik is prekondisionering nodig om die konvergensie te verbeter. Die gebruik van prekondisionering gebasseer op skuif-en-omkeer word bespreek en numeriese vergelykings met direkte oplossers word aangebied. As toepassing los ons ’n twee-dimensionele weergawe van die hittevergelyking op, met ’n konveksie term en daarsonder. Ons numeriese eksperimente dui aan dat die Krylov metodes met prekondisionering effektief is, beide in terme van berekeningstyd en rekenaargeheue, in vergelyking met die direkte metodes.

Today the multisensor array signal processing of noisy measurements has received much attention. The classical problem in array signal processing is determining the location of an energy-radiating source relative to the location of the array, in other words, direction-of-arrival (DOA) estimation. One is considering the signal estimation problem when together with the signal(s) of interest some noise and interfering signals are present. In this report a direction-of-arrival estimation system is described based on an antenna array for detecting arrival angles in azimuth plane of signals pitched by the antenna array. For this, the Multiple Signal Classication (MUSIC) algorithmis first of all considered. Studies show that in spite of its good reputation and popularity among researches, it has a certain limit of its performance. In this subspace-based method for DOA estimation of signal wavefronts, the term corresponding to additive noise is initially assumed spatially white. In our paper, we address the problem of DOA estimation of multiple target signals in a particular noise situation - in correlated spherically isotropic noise, which, in many practical cases, models a more real context than under the white noise assumption. The purpose of this work is to analyze the behaviour of the MUSIC algorithm and compare its performance with some other algorithms (such as the Capon and the Classical algorithms) and, uppermost, to explore the quality of the detected angles in terms of precision depending on different parameters, e.g. number of samples, noise variance, number of incoming signals. Some modifications of the algorithms are also done is order to increase their performance. Program MATLAB is used to conduct the studies. The simulation results on the considered antenna array system indicate that in complex conditions the algorithms in question (and first of all, the MUSIC algorithm) are unable to automatically detect and localize the DOA signals with high accuracy. Other algorithms andways for simplification the problem (for example, procedure of denoising) exist and may provide more precision but require more computation time.

We propose a new framework for understanding block Krylov subspace methods, which hinges on a matrix-valued inner product. We can recast the ``classical" block Krylov methods, such as O'Leary's block conjugate gradients, global methods, and loop-interchange methods, within this framework. Leveraging the generality of the framework, we develop an efficient restart procedure and error bounds for the shifted block full orthogonalization method (Sh-BFOM(m)). Regarding BFOM as the prototypical block Krylov subspace method, we propose another formalism, which we call modified BFOM, and show that block GMRES and the new block Radau-Lanczos method can be regarded as modified BFOM. In analogy to Sh-BFOM(m), we develop an efficient restart procedure for shifted BGMRES with restarts (Sh-BGMRES(m)), as well as error bounds. Using this framework and shifted block Krylov methods with restarts as a foundation, we formulate block Krylov subspace methods with restarts for matrix functions acting on multiple vectors f(A)B. We obtain convergence bounds for \bfomfom (BFOM for Functions Of Matrices) and block harmonic methods (i.e., BGMRES-like methods) for matrix functions. With various numerical examples, we illustrate our theoretical results on Sh-BFOM and Sh-BGMRES. We also analyze the matrix polynomials associated to the residuals of these methods. Through a variety of real-life applications, we demonstrate the robustness and versatility of B(FOM)^2 and block harmonic methods for matrix functions. A particularly interesting example is the tensor t-function, our proposed definition for the function of a tensor in the tensor t-product formalism. Despite the lack of convergence theory, we also show that the block Radau-Lanczos modification can reduce the number of cycles required to converge for both linear systems and matrix functions. / Mathematics

Les travaux de ce doctorat concernent le développement de méthodes itératives pour la résolution de systèmes linéaires creux de grande taille comportant de nombreux seconds membres. L’application visée est la résolution d’un problème inverse en géophysique visant à reconstruire la vitesse de propagation des ondes dans le sous-sol terrestre. Lorsque de nombreuses sources émettrices sont utilisées, ce problème inverse nécessite la résolution de systèmes linéaires complexes non symétriques non hermitiens comportant des milliers de seconds membres. Dans le cas tridimensionnel ces systèmes linéaires sont reconnus comme difficiles à résoudre plus particulièrement lorsque des fréquences élevées sont considérées. Le principal objectif de cette thèse est donc d’étendre les développements existants concernant les méthodes de Krylov par bloc. Nous étudions plus particulièrement les techniques de déflation dans le cas multiples seconds membres et recyclage de sous-espace dans le cas simple second membre. Des gains substantiels sont obtenus en terme de temps de calcul par rapport aux méthodes existantes sur des applications réalistes dans un environnement parallèle distribué. / This PhD thesis concerns the development of flexible Krylov subspace iterative solvers for the solution of large sparse linear systems of equations with multiple right-hand sides. Our target application is the solution of the acoustic full waveform inversion problem in geophysics associated with the phenomena of wave propagation through an heterogeneous model simulating the subsurface of Earth. When multiple wave sources are being used, this problem gives raise to large sparse complex non-Hermitian and nonsymmetric linear systems with thousands of right-hand sides. Specially in the three-dimensional case and at high frequencies, this problem is known to be difficult. The purpose of this thesis is to develop a flexible block Krylov iterative method which extends and improves techniques already available in the current literature to the multiple right-hand sides scenario. We exploit the relations between each right-hand side to accelerate the convergence of the overall iterative method. We study both block deflation and single right-hand side subspace recycling techniques obtaining substantial gains in terms of computational time when compared to other strategies published in the literature, on realistic applications performed in a parallel environment.

Miniaturisation of electronic chips poses challenges at the design stage. The progressively decreasing circuit dimensions result in complex electrical behaviour that necessitates complex models. Simulation of complex circuit models involves extraordinarily large compu- tational complexity. Such complexity is better managed through Model Order Reduction. Model order reduction has been successful in large reductions in system order for most types of circuits, at high levels of accuracy. However, multiport circuits with large number of inputs/outputs, pose an additional computational challenge. A strategy based on exible clustering of interconnects results in more e cient reduction of multiport circuits. Clustering methods traditionally use Krylov-subspace methods such as PRIMA for the actual model reduction step. These clustering methods are unable to reduce the model order to the optimum extent. SVD-based methods like Truncated Balanced Realization have shown higher reduction potential than Krylov-subspace methods. In this thesis, the di erences in reduction potential and computational cost thereof between SVD-based methods and Krylov-subspace methods are identi ed, analyzed and quanti ed. A novel algorithm has been developed, utilizing a particular combination of both these methods to achieve better results. It enhances the clustering method for model reduction using Truncated Balanced Realization as a second-level reduction technique. The algorithm is tested and signi cant gains are illustrated. The proposed novel algorithm preserves the other advantages of the current clustering algorithm.

La performance d'un algorithme sur une architecture donnée dépend à la fois de la vitesse à laquelle le processeur effectue des opérations à virgule flottante (flops) et de la vitesse d'accès à la mémoire et au disque. Etant donné que le coût de la communication est beaucoup plus élevé que celui des opérations arithmétiques, celle-là forme un goulot d'étranglement dans les algorithmes numériques. Récemment, des méthodes de sous-espace de Krylov basées sur les méthodes 's-step' ont été développées pour réduire les communications. En effet, très peu de préconditionneurs existent pour ces méthodes, ce qui constitue une importante limitation. Dans cette thèse, nous présentons le préconditionneur nommé ''Communication-Avoiding ILU0'', pour la résolution des systèmes d’équations linéaires (Ax=b) de très grandes tailles. Nous proposons une nouvelle renumérotation de la matrice A ('alternating min-max layers'), avec laquelle nous montrons que le préconditionneur en question réduit la communication. Il est ainsi possible d’effectuer « s » itérations d’une méthode itérative préconditionnée sans communication. Nous présentons aussi deux nouvelles méthodes itératives, que nous nommons 'multiple search direction with orthogonalization CG' (MSDO-CG) et 'long recurrence enlarged CG' (LRE-CG). Ces dernières servent à la résolution des systèmes linéaires d’équations de très grandes tailles, et sont basées sur l’enrichissement de l’espace de Krylov par la décomposition du domaine de la matrice A. / The performance of an algorithm on any architecture is dependent on the processing unit’s speed for performing floating point operations (flops) and the speed of accessing memory and disk. As the cost of communication is much higher than arithmetic operations, and since this gap is expected to continue to increase exponentially, communication is often the bottleneck in numerical algorithms. In a quest to address the communication problem, recent research has focused on communication avoiding Krylov subspace methods based on the so called s-step methods. However there are very few communication avoiding preconditioners, and this represents a serious limitation of these methods. In this thesis, we present a communication avoiding ILU0 preconditioner for solving large systems of linear equations (Ax=b) by using iterative Krylov subspace methods. Our preconditioner allows to perform s iterations of the iterative method with no communication, by applying a heuristic alternating min-max layers reordering to the input matrix A, and through ghosting some of the input data and performing redundant computation. We also introduce a new approach for reducing communication in the Krylov subspace methods, that consists of enlarging the Krylov subspace by a maximum of t vectors per iteration, based on the domain decomposition of the graph of A. The enlarged Krylov projection subspace methods lead to faster convergence in terms of iterations and to parallelizable algorithms with less communication, with respect to Krylov methods. We discuss two new versions of Conjugate Gradient, multiple search direction with orthogonalization CG (MSDO-CG) and long recurrence enlarged CG (LRE-CG).

Title: Krylov Subspace Methods - Analysis and Application Author: Tomáš Gergelits Department: Department of Numerical Mathematics Supervisor: prof. Ing. Zdeněk Strakoš, DrSc., Department of Numerical Mathematics Abstract: Convergence behavior of Krylov subspace methods is often studied for linear algebraic systems with symmetric positive definite matrices in terms of the condition number of the system matrix. As recalled in the first part of this thesis, their actual convergence behavior (that can be in practice also substantially affected by rounding errors) is however determined by the whole spectrum of the system matrix, and by the projections of the initial residual to the associated invariant subspaces. The core part of this thesis investigates the spectra of infinite dimensional operators −∇ · (k(x)∇) and −∇ · (K(x)∇), where k(x) is a scalar coefficient function and K(x) is a symmetric tensor function, preconditioned by the Laplace operator. Subsequently, the focus is on the eigenvalues of the matrices that arise from the discretization using conforming finite elements. Assuming continuity of K(x), it is proved that the spectrum of the preconditi- oned infinite dimensional operator is equal to the convex hull of the ranges of the diagonal function entries of Λ(x) from the spectral decomposition K(x) =...

In this thesis, a method to calculate eigenpairs is implemented for the Multipsi library. While the standard implemtentations use the Davidson method with Rayleigh-Ritz extraction to calculate the eigenpairs with the lowest eigenvalues, the new method uses the harmonic Davidson method with the harmonic Rayleigh-Ritz extraction to calculate eigenpairs with eigenvalues near a chosen target. This is done for Configuration Interaction calculations and for Multiconfigurational methods. From calculations, it seems the new addition to the Multipsi library is worth investigating further as convergence for difficult systems with a lot of near-degeneracy was improved. / I denna avhandling implementeras en metod för att beräkna egenpar för Multipsi-biblioteket. Medan standardimplementeringarna använder Davidson-metoden med Rayleigh-Ritz-extraktion för att beräkna egenparen med de lägsta egenvärdena, använder den nya metoden den harmoniska Davidson-metoden med den harmoniska Rayleigh-Ritz-extraktionen för att beräkna egenparen med egenvärden nära ett valt mål. Detta görs för konfigurationsinteraktionsberäkningar och för multikonfigurationsmetoder. Utifrån beräkningarna verkar det nya tillskottet till Multipsi-biblioteket vara värt att undersöka vidare eftersom konvergensen för svåra system med mycket nära degenerering förbättrades.