| Hlavní výstupy projektu |
Metoda sdružených gradientů (CG) [1] je jednou z nejpopulárnějších metod pro řešení velkých lineárních systémů Mx=b se symetrickou pozitivně definitivní (SPD) maticí M. V naší práci jsme zkoumali použití metody CG pro řešení tzv. Lyapunovovy rovnice AX+XA^T = K, kde A, X, K jsou čtvercové matice stejných rozměrů a A^T označuje transpozi matice. Tyto rovnice hrají důležitou roli v dynamických systémech a teorii kontroly. V aplikacích typicky máme K = BB^T s maticí B obsahující pouze několik sloupců [2]. Pro zjednodušení navíc předpokládáme, že A je symetrická pozitivně definitivní.
Pomocí Kroneckerových součinů se může Lyapunovova rovnice AX+XA^T = K přepsat jako standardní lineární systém tvaru Mx = b s M = kron(A,I)+kron(I,A), x = vec(X) a b = vec(K), kde kron označuje Kroneckerův součin dvou matic a vec seřazuje sloupce matice do dlouhého vektoru. Za předpokladů uvedených výše M je matice SPD, a proto CG lze použít, alespoň v zásadě. Pro větší velikost matice n se však metoda CG rychle stává nepraktickou, jednoduše z důvodu potřeby uložení vektorů délky n*n. Problémy objevující se v praxi často obsahují n v řádu milionů nebo miliard. Na základě myšlenek a algoritmů v [3] jsme vyvinuli variantu metody CG s nízkou hodností, která tento problém eliminuje. Základní myšlenkou za naší variantou je to, že řešení X Lyapunovovy rovnice lze zobrazit tak, že má exponenciální rozklad ve svých singulárních hodnotách. Tedy X lze aproximovat velmi přesně maticí s nízkou hodností vyžadující značně méně uložení. Naše varianta metody CG přistupuje k této aproximace s nízkou hodností přes iteraci s nízkou hodností, což implikuje, že požadavky na skladování rostou lineárně nebo logaritmicko-lineárně (místo kvadraticky) jako n nárůsty. Ukázali jsme, že metoda CG i naše varianta s nízkou hodností si zachovávají symetrii řešení. Zůstává ovšem obtížným otevřeným problémem také kontrolovat rozklad singulárních hodnot (a tedy kvalitu aproximace s nízkou hodností) pomocných iterací. Máme silný numerický důkaz tohoto silného rozkladu za předpokladu, že metoda CG konverguje dostatečně rychle. K dosažení rychlé konvergence je potřeba použít velmi efektivní předpoklady.
Pro získání smysluplné podmíněné CG metody s nízkou hodností musí předpoklad respektovat struktury s nízkou hodností. Testovali jsme několik předpokladů, které splňují tento požadavek. Dva z našich předpokladů jsou založeny na stacionární SSOR iteraci a metodě ADI, jak je navrženo v [4] pro jinou metodu. Třetí předpoklad, který jsme testovali, je založen na iteraci funkce signum. Implementace tohoto předpokladu je obtížná, protože vyžaduje vytvořit následné částky matice A a její inverzi. I když A je řídký, tyto sumy jsou obvykle plně zaplněny. Proto není možné je skladovat explicitně pro velké problémy. Pro matice vznikající z diskretizací konečného elementu se mohou tyto sumy často účinně skladovat jako hierarchické matice nebo jejich varianty. Použili jsme dvě existující implementace, HLib c-library [5, 6] a Matlab toolbox pro hierarchické polooddělitelné (HSS) matice [7]. Na základě numerických výsledků jsme se rozhodli preferovat jen jeden krok metody ADI jako předpoklad a jeden krok iteraci signum. Ukázalo se, že SSOR předpoklad se dobře hodí pro aritmetiku s nízkou hodností. Podrobné výsledky výkonu pro ADI a metody funkce signum jsou uvedeny ve zprávě [8].
Při výrazném rozšíření metod diskutovaných výše jsme vyvinuli tenzorovou CG metodu s nízkou hodností pro řešení Lyapunovovy rovnice závislé na parametrech A(p)X(p)+X(p)A(p)^T = – B(p)B(p)^T , kde p je parametr nebo vektor několika parametrů. Pro případ jednoho parametru iterace v metodě CG může být reprezentován tenzory řádu 3. Místo matic s nízkou hodností se používá symetrická varianta Tuckerova formátu k uložení a manipulování s iterací. Kvalitu aproximace seříznutí s nízkou hodností v Tuckerově formátu lze kvantifikovat rozkladem singulární hodnoty matricizací. Ve větě 3.2 v [8] jsme prokázali ohraničení na tomto rozkladu pro tenzor obsahující příklady řešení Lyapunovovy rovnice závislé na parametrech. Důkaz používá novou kombinaci existujících ohraničení rozkladu pro Lyapunovovy rovnice a lineární systémy závislé na parametrech. Kvalitativní rozklad predikovaný našimi ohraničeními se dobře hodí pro numerické experimenty. Tento výsledek poskytuje teoretický základ pro numericky pozorovanou účinnost kombinování seříznutí s nízkou hodností v Tuckerově formátu s metodou CG. Předpoklad používá metodu ADI pro průměr příkladů parametrů.
Na závěr jsme vytvořili metodu CG pro řešení Lyapunovovy rovnice v závislosti na několika, řekněme d>=2, parametrech. Generalizace případu o jednom parametru se ukázala být relativně přímá, a to kombinováním vhodného formátu tenzoru s nízkou hodností a vyššího řádu (low-rank hihger-order) s metodou CG. V naší implementaci jsme použili tzv. hierarchický Tuckerův formát [12, 13], ale formát –tenzor-řada [14, 15] by byl stejně vhodný. Numerické experimenty v [8] odhalují, že tento přístup může být účinně použit k řešení Lyapunovových rovnic v závislosti na několika parametrech pro velký počet vzorků parametrů. Řešení Lyapunovových rovnic závislých na parametrech hraje významnou úlohu při interpolačních přístupech k redukci modelů parametrizovaných lineárních kontrolních systémů [16]. Metody vytvořené v tomto projektu představují první důležitý krok směrem k účinným algoritmům pro tyto techniky redukce modelů.
Všechny vytvořené algoritmy byly implementovány v Matlab a jsou k dispozici na požádání. Během práce na projektu příjemce grantu úzce spolupracoval s Prof. Danielem Kressnerem a členy jeho skupiny (Cedric Effenberger, Dr. Effrosyni Kokiopoulou, Christine Tobler, Dr. Bart Vandereycken).
[1] M. R. Hestenes and E. Stiefel: Methods of Conjugate Gradients for Solving Linear Systems, Journal of Research of the National Bureau of Standards 49 (6), pp. 409-436 (1952).
[2] A. C. Antoulas: Approximation of large-scale dynamical systems, SIAM publications, Philadelphia PA, 2005.
[3] D. Kressner and C. Tobler: Low-rank tensor Krylov subspace methods for parametrized linear systems, SIAM J. Matrix Anal. Appl. 32, pp. 1288-1316 (2011).
[4] M. Hochbruck and G. Starke: Preconditioned Krylov subspace methods for Lyapunov matrix equation, SIAM J. Matrix Anal. Appl. 16, pp. 156-171 (1995).
[5] HLib library, http://www.hlib.org.
[6] W. Hackbusch, Hierarchische Matrizen: Algorithmen und Analysis, Springer, Berlin, 2009.
[7] S. Pauli: A numerical solver for Lyapunov equations based on the matrix sign function iteration in HSS arithmetic, 2010, Semester thesis, SAM, ETH Zurich.
[8] D. Kressner, M. Plesinger, and C. Tobler: A preconditioned low-rank CG method for parameter-dependent Lyapunov matrix equations, submitted to Numerical Linear Algebra with Applications, available as MATHICSE technical report 18.2012
(http://mathicse.epfl.ch/page-68906-en.html).
[9] L. R. Tucker: Implications of factor analysis of three-way matrices for measurement of change, in Problem in Measuring Change, C. W. Harris, ed., University of Wiskonsin, New York, pp. 122-137, 1963.
[10] L. R. Tucker: The extension of factor analysis to three-dimensional matrices, in Contribution to Mathematical Psychology, H. Gulliksen and N. Frederiksen, eds., Holt, Rinehardt, & Winston, New York, pp. 110-127, 1964.
[11] L. R. Tucker: Some mathematical notes on three-mode factor analysis, Psychometrika 31, pp. 279-311 (1966).
[12] L. Grasedyck: Hierarchical singular value decomposition of tensors, SIAM J. on Matrix Analysis and Applications 31, pp. 2029-2054 (2010).
[13] D. Kressner and C. Tobler: htucker-a Matlab toolbox for tensors in hierarchical Tucker format, SAM Research Report No. 2012-2.
[14] I. V. Oseledets and E. E. Tyrtyshnikov: Recursive decomposition of multidimensional tensors, Doklady Mathematics 80, pp. 460-462 (2009).
[15] I. V. Oseledets: Tensor train decomposition, SIAM J. on Scientific Computing 33, pp. 2295-2317.
|
|---|
