Balancing three matrices in control theory

Several problems from control theory are presented which are sensitive to badly scaled matrices. We were specially concerned with the algorithms involving three matrices, thus we extended the Ward\u27s balancing algorithm for two matrices. Numerical experiments confirmed that balancing three matrices can produce an accurate frequency response matrix for descriptor linear systems, it can also improve the solution of the pole assignment problem via state feedback and the determination of the controllable part of the system

Michael James David Powell:29 July 1936-19 April 2015

Michael James David Powell was a British numerical analyst who was among the pioneers of computational mathematics. During a long and distinguished career, first at the Atomic Energy Research Establishment (AERE) Harwell and subsequently as the John Humphrey Plummer Professor of Applied Numerical Analysis in Cambridge, he contributed decisively towards establishing optimization theory as an effective tool of scientific enquiry, replete with highly effective methods and mathematical sophistication. He also made crucial contributions to approximation theory, in particular to the theory of spline functions and of radial basis functions. In a subject that roughly divides into practical designers of algorithms and theoreticians who seek to underpin algorithms with solid mathematical foundations, Mike Powell refused to follow this dichotomy. His achievements span the entire range from difficult and intricate convergence proofs to the design of algorithms and production of software. He was among the leaders of a subject area that is at the nexus of mathematical enquiry and applications throughout science and engineering.</jats:p

Computing the singular value decomposition with high relative accuracy

AbstractWe analyze when it is possible to compute the singular values and singular vectors of a matrix with high relative accuracy. This means that each computed singular value is guaranteed to have some correct digits, even if the singular values have widely varying magnitudes. This is in contrast to the absolute accuracy provided by conventional backward stable algorithms, which in general only guarantee correct digits in the singular values with large enough magnitudes. It is of interest to compute the tiniest singular values with several correct digits, because in some cases, such as finite element problems and quantum mechanics, it is the smallest singular values that have physical meaning, and should be determined accurately by the data. Many recent papers have identified special classes of matrices where high relative accuracy is possible, since it is not possible in general. The perturbation theory and algorithms for these matrix classes have been quite different, motivating us to seek a common perturbation theory and common algorithm. We provide these in this paper, and show that high relative accuracy is possible in many new cases as well. The briefest way to describe our results is that we can compute the SVD of G to high relative accuracy provided we can accurately factor G=XDYT where D is diagonal and X and Y are any well-conditioned matrices; furthermore, the LDU factorization frequently does the job. We provide many examples of matrix classes permitting such an LDU decomposition

Robusne numeričke metode za nelinearne probleme svojstvenih vrijednosti

In this thesis we study numerical methods for solving nonlinear eigenvalue problems of polynomial type, i.e. $P(\lambda)x \equiv (\Sigma_{\ell=0}^k \lambda^{\ell} A_{\ell})x = 0$, where $A_{\ell} \in \mathbb{C}^{n \times n}, \lambda \in \mathbb{C}, 0 \neq x \in \mathbb{C}^n$. In particular, we are interested in the quadratic $(k = 2)$ and the quartic $(k = 4)$ eigenvalue problems. The methods are based on the corresponding linearization – the nonlinear problem is replaced with an equivalent linear problem of the type $(A - \lambda B)y = 0$, of dimension $kn$. We propose several modifications and improvements of the existing methods for both the complete and partial solution; this results in new numerical algorithms that are a substantial improvement over the existing ones. In particular, as an improvement of the state of the art quadeig method of Hammarling, Munro and Tisseur, we develop a scheme to deflate all zero and infinite eigenvalues before calling the QZ algorithm for the linear problem. This provides numerically more robust procedure, which we illustrate by numerical examples. Further, we supplement the parameter scaling (designed to equilibrate the norms of the coefficient matrices) with a two–sided diagonal scaling to nearly equilibrate (in modulus) the nonzero matrix entries. In addition, we analyze the fine details of the rank revealing factorization used in the deflation process. We advocate to use complete pivoting in the QR factorization, and we also propose a LU based approach, which is shown to be competitive, or even better than the one based on the QR factorization. The new method is extended to the quartic problem. For the partial quadratic eigenvalue problem (computing only a part of the spectrum), the iterative Arnoldi–like methods are studied, especially the implicitly restarted two level orthogonal Arnoldi algorithm (TOAR). We propose several improvements of the method. In particular, new shift selection strategy is proposed for the implicit restart for the class of overdamped quadratic eigenvalue problems. Also, we show the benefit of choosing the starting vector for TOAR, based on spectral information of a nearby proportionally damped pencil. Finally, we provide some new ideas for the development of a Krylov–Schur like methods that is capable of using arbitrary polynomial filters in the implicit restarting.Nelinearni problemi svojstvenih vrijednosti se javljaju u mnogim primjenama kako u prirodnim znanostima, tako i u inženjerstvu. Jedna od najpoznatijih klasa nelinearnih svojstvenih problema su polinomni svojstveni problemi. Tako se, na primjer, kvadratični svojstveni problem $(\lambda^2 M + \lambda C + K)x = 0$ pojavljuje u dinamičkoj analizi mehaničkih i električnih struktura, u vibro–akustici, mehanici fluida, obradi signala. S druge strane, polinomni se problem četvrtog reda $(\lambda^4 A + \lambda^3 B + \lambda^2 C + \lambda D + K)x = 0$ pojavljuje u analizi stabilnosti Poiseuilleovog toka u cijevi. Za razliku od linearnih problema svojstvenih vrijednosti, numeričke metode za nelinearne probleme još uvijek nisu dovoljno razrađene, niti numerički pouzdane, iako je algebarska teorija za polinomne probleme svojstvenih vrijednosti dobro razvijena. Naglasak ove disertacije je na numeričkom rješavanju kvadratičnog svojstvenog problema. Cilj je razviti nove, robusnije numeričke metode koje se mogu koristiti u praksi kao pouzdan numerički softver. U disertaciji se proučavaju dvije vrste metoda: direktne i iterativne. Direktne metode se razvijaju za računanje svih svojstvenih vrijednosti i odgovarajućih svojstvenih vektora zadanog problema. Kada nas zanima samo dio spektra, recimo one svojstvene vrijednosti koje su najveće po modulu ili one koje se nalaze u lijevoj kompleksnoj poluravnini, tada koristimo iterativne metode. Ovdje je najšešće slučaj da je dimenzija originalnog problema mnogo veća od broja svojstvenih vrijednosti koje želimo izračunati. Ideja iterativnih metoda je konstruirati potprostor mnogo manje dimenzije od originalnog problema koji sadrži informaciju o traženom dijelu spektra, a aproksimacija traženog dijela spektra se onda izračuna koristeći projekciju problema na nađeni potprostor. Osnova većine metoda za rješavanje polinomnih svojstvenih problema je linearizacija, to jest polinomni problem se zamijeni ekvivalentnim linearnim problemom koji se onda rješava koristeći već razvijene metode za linearne probleme. Međutim, naivno direktno korištenje linearnih metoda ne garantira zadovoljavajuće rezultate za originalni problem. Čak i ako izračunati svojstveni par ima malu grešku unazad za odgovarajuću linearizaciju, greška unazad za rekonstruirani svojstveni par originalnog problema može biti puno veća. Prije razvijanja metoda, u Poglavlju 2 je predstavljena analiza grešaka unazad za polinomni svojstveni problem, bazirana na radu F. Tisseur [66]. Ideja analize grešaka unazad je da se izračunate aproksimacije interpretiraju kao egzaktna rješenja problema koji je blizu originalnom problemu, i čiji matrični koeficijenti su definirani kao $A_\ell + \Delta A_\ell$ pri čemu je $\Delta A_\ell$ malo. Međutim, u mnogim primjenama matrice $A_\ell$ imaju određenu strukturu, npr. hermitske su, ili anti hermitske. Prema tome, bilo bi prirodno zahtijevati da greška unazad $\Delta A_\ell$ čuva ovu strukturu. U slučaju kad je ta struktura hermitska i anti hermitska, postojeći rezultati za realne svojstvene vrijednosti su prošireni na općenite svojstvene vrijednosti. U poglavlju 3 se proučavaju direktne metode za rješavanje kvadratičnog svojstvenog problema. Standardni pristup je korištenje QZ algoritma na odgovarajućoj linearizaciji. Međutim, ako originalni problem ima svojstvene vrijednosti koje su nula ili beskonačno, ovakav pristup je sklon numeričkim poteškoćama. 2011. Hammarling, Munro i Tisseur [37] su razvili quadeig algoritam koji prije korištenja QZ metode za linearni problem skalira originalni problem kako bi norme matričnih koeficijenata bile ujednačene te pokuša detektirati postojanje svojstvenih vrijednosti nula i beskonačno koje ona procesom deflacije ukloni iz linearizacije. Deflacija se temelji na određivanju ranga matrica M i K. Kod quadeiga se koristi QR faktorizacija pivotiranjem stupaca. Koristeći ortogonalne transformacije $n-rank(M)$ beskonačnih i $n-rank(K)$ svojstvenih vrijednosti nula je uklonjeno iz odgovarajuće linearizacije. Glavni doprinos ovog poglavlja je novi algoritam za nalaženje svih svojstvenih vrijednosti kvadratično problema kojeg zovemo KVADeig. Kao motivacija za potrebu poboljšanja quadeiga je predstavljen primjer kod kojeg quadeig nije uspio detektirati sve beskonačne svojstvene vrijednosti. Štoviše, nakon što je uklonjen određen broj ovih svojstvenih vrijednosti, preostale izračunate svojstvene vrijednosti koje su konačne čak nemaju ni veliku apsolutnu vrijednost koja bi nas možda mogla nagnati na zaključak da bi one trebale biti proglašene beskonačnim. Problem nastane kada postoji više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. Naime, deflacija u quadeigu ukloni samo jedan Jordanov blok. Kako bismo riješili ovaj problem razvili smo test koji služi za provjeru postoji li više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. On je baziran na Van Doorenovom algoritmu za određivanje Kroneckerove strukture generaliziranog svojstvenog problema. Dodatno se analizira utjecaj metoda koje se koriste kao faktorizacije za određivanje ranga te utjecaj kriterija po kojem se rang određuje. Pored skaliranja koje je predloženo u quadeigu uvodimo i dvostrano dijagonalno balansiranje čiji je cilj ujednačavanje elemenata u matricama koje definiraju problem. Na kraju razvijamo metodu baziranu na LU faktorizaciji potpunim pivotiranjem za određivanje ranga. Numerički eksperimenti u Sekciji 3.7 ilustriraju prednosti predložene metode. U poglavlju 4 je razvijen novi algoritam KVARTeig za rješavanje polinomnog svojstvenog problema stupnja četiri. Umjesto direktne linearizacije koristimo kvadratifikaciju koja je uvedena u [17], tj. definiramo ekvivalentan kvadratični problem. Novi algoritam je baziran na KVADeigu, s tim da je skaliranje definirano na matricama originalnog problema i proces deflacije je prilagođen tako da što više iskoristi strukturu originalnog problema. Kao i za kvadratični problem, i ovdje je razvijen test za provjeru postojanja više od jednog Jordanovog bloka za svojstvene vrijednosti nula i beskonačno. Numerički primjeri u Sekciji 4.5 prikazuju prednost nove metode nad quadeigom i polyeigom koji je implementiran u MATLABu. U Poglavlju 5 se proučavaju iterativne metode Arnoldijevog tipa za kvadratični svojstveni problem. Bai i Su [3] su prvi primijetili da je u slučaju iterativnih metoda Arnoldijevog tipa bolje primijeniti Rayleigh–Ritzovu projekciju direktno na originalni kvadratični problem. U tu svrhu su definirani Krilovljev potprostor drugog reda i odgovarajući algoritam SOAR (Second Order Arnoldi) za računanje odgovarajuće baze. Ovaj algoritam je dodatno modificiran te je razvijen takozvani TOAR (Two level orthogonal Arnoldi) algoritam [49]. U ovom poglavlju predlažemo nekoliko modifikacija implicitno restartanog TOAR algoritma koje su temeljene na činjenici da algoritam koristimo za rješavanje kvadratičnog problema svojstvenih vrijednosti. Pod implicitnim restartanjem se misli na korištenje polinomih filtera kako bi se definirao novi početni vektor koji uvelike utječe na konvergenciju metode. Za posebnu klasu pregušenih problema svojstvenih vrijednosti predlažemo novi način definiranja polinomnih filtera. Također, za općenite probleme, predlažemo novi izbor početnog vektora koji se temelji na aproksimaciji kvadratičnog svojstvenog problema problemom čije je gušenje linearno. Numerički primjeri pokazuju da predložene modifikacije rezultiraju manjim brojem restartanja potrebnih za nalaženje svojstvenih parova sa zadovoljavajućom greškom unatrag. U drugom dijelu Poglavlja 5 dajemo pregled implicitno restartanog Krylov–Schurovog algoritma kojeg je uveo Stewart [64]. Ideja ovog algoritma je da se definira faktorizacija koja ne zahtijeva posebnu strukturu kao Arnoldijeva, i na koju će se lakše primijeniti implicitno restartanje. Međutim, prilikom ovakvog restartanja moguće je koristiti samo egzaktne pomake za definiranje polinomnog filtera. Drmač i Bujanović su razvili metodu koja omogućava korištenje proizvoljnih pomaka kod implicitno restartanog Krylov–Schurovog algoritma. U ovom poglavlju generaliziramo predloženi proces u svrhu korištenja Krylov–Schurovog algoritma za rješavanje kvadratičnog svojstvenog problema

Berechnung und Anwendungen Approximativer Randbasen

This thesis addresses some of the algorithmic and numerical challenges associated with the computation of approximate border bases, a generalisation of border bases, in the context of the oil and gas industry. The concept of approximate border bases was introduced by D. Heldt, M. Kreuzer, S. Pokutta and H. Poulisse in "Approximate computation of zero-dimensional polynomial ideals" as an effective mean to derive physically relevant polynomial models from measured data. The main advantages of this approach compared to alternative techniques currently in use in the (hydrocarbon) industry are its power to derive polynomial models without additional a priori knowledge about the underlying physical system and its robustness with respect to noise in the measured input data. The so-called Approximate Vanishing Ideal (AVI) algorithm which can be used to compute approximate border bases and which was also introduced by D. Heldt et al. in the paper mentioned above served as a starting point for the research which is conducted in this thesis. A central aim of this work is to broaden the applicability of the AVI algorithm to additional areas in the oil and gas industry, like seismic imaging and the compact representation of unconventional geological structures. For this purpose several new algorithms are developed, among others the so-called Approximate Buchberger Möller (ABM) algorithm and the Extended-ABM algorithm. The numerical aspects and the runtime of the methods are analysed in detail - based on a solid foundation of the underlying mathematical and algorithmic concepts that are also provided in this thesis. It is shown that the worst case runtime of the ABM algorithm is cubic in the number of input points, which is a significant improvement over the biquadratic worst case runtime of the AVI algorithm. Furthermore, we show that the ABM algorithm allows us to exercise more direct control over the essential properties of the computed approximate border basis than the AVI algorithm. The improved runtime and the additional control turn out to be the key enablers for the new industrial applications that are proposed here. As a conclusion to the work on the computation of approximate border bases, a detailed comparison between the approach in this thesis and some other state of the art algorithms is given. Furthermore, this work also addresses one important shortcoming of approximate border bases, namely that central concepts from exact algebra such as syzygies could so far not be translated to the setting of approximate border bases. One way to mitigate this problem is to construct a "close by" exact border bases for a given approximate one. Here we present and discuss two new algorithmic approaches that allow us to compute such close by exact border bases. In the first one, we establish a link between this task, referred to as the rational recovery problem, and the problem of simultaneously quasi-diagonalising a set of complex matrices. As simultaneous quasi-diagonalisation is not a standard topic in numerical linear algebra there are hardly any off-the-shelf algorithms and implementations available that are both fast and numerically adequate for our purposes. To bridge this gap we introduce and study a new algorithm that is based on a variant of the classical Jacobi eigenvalue algorithm, which also works for non-symmetric matrices. As a second solution of the rational recovery problem, we motivate and discuss how to compute a close by exact border basis via the minimisation of a sum of squares expression, that is formed from the polynomials in the given approximate border basis. Finally, several applications of the newly developed algorithms are presented. Those include production modelling of oil and gas fields, reconstruction of the subsurface velocities for simple subsurface geometries, the compact representation of unconventional oil and gas bodies via algebraic surfaces and the stable numerical approximation of the roots of zero-dimensional polynomial ideals