A hierarchically blocked Jacobi SVD algorithm for single and multiple graphics processing units

We present a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD), targeting both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of GPU's memory hierarchy. The algorithm may outperform MAGMA's dgesvd, while retaining high relative accuracy. To this end, we developed a family of parallel pivot strategies on GPU's shared address space, but applicable also to inter-GPU communication. Unlike common hybrid approaches, our algorithm in a single GPU setting needs a CPU for the controlling purposes only, while utilizing GPU's resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a Tesla S2050 system with four GPUs vs. a single Fermi card.Comment: Accepted for publication in SIAM Journal on Scientific Computin

High performance Cholesky and symmetric indefinite factorizations with applications

The process of factorizing a symmetric matrix using the Cholesky (LLT ) or indefinite (LDLT ) factorization of A allows the efficient solution of systems Ax = b when A is symmetric. This thesis describes the development of new serial and parallel techniques for this problem and demonstrates them in the setting of interior point methods. In serial, the effects of various scalings are reported, and a fast and robust mixed precision sparse solver is developed. In parallel, DAG-driven dense and sparse factorizations are developed for the positive definite case. These achieve performance comparable with other world-leading implementations using a novel algorithm in the same family as those given by Buttari et al. for the dense problem. Performance of these techniques in the context of an interior point method is assessed

Basis preconditioning in interior point methods

Solving normal equations AAᵀx = b, where A is an m x n matrix, is a common task in numerical optimization. For the efficient use of iterative methods, this thesis studies the class of preconditioners of the form BBᵀ , where B is a nonsingular "basis" matrix composed of m columns of A. It is known that for any matrix A of full row rank B can be chosen so that the entries in [B⁻¹A] are bounded by 1. Such a basis is said to have "maximum volume" and its preconditioner bounds the spectrum of the transformed normal matrix in the interval [1, 1+mn]. The theory is extended to (numerically) rank deficient matrices, yielding a rank revealing variant of Gaussian elimination and a method for computing the minimum norm solution for x from a reduced normal system and a low-rank update. Algorithms for finding a maximum volume basis are discussed. In the linear programming interior point method a sequence of normal equations needs to be solved, in which A changes by a column scaling from one system to the next. A heuristical algorithm is proposed for maintaining a basis of approximate maximum volume by update operations as those in the revised simplex method. Empirical results demonstrate that the approximation means no loss in the effectiveness of the preconditioner, but makes basis selection much more efficient. The implementation of an interior point solver based on the new linear algebra is described. Features of the code include the elimination of free variables during preconditioning and the removal of degenerate variables from the optimization process once sufficiently close to a bound. A crossover method recovers a vertex solution to the linear program, starting from the basis at the end of the interior point solve. A computational study shows that the implementation is robust and of general applicability, and that its average performance is comparable to that of state-of-the-art solvers

Paralelni algoritmi Jacobijeva tipa za singularnu i generaliziranu singularnu dekompoziciju

In this thesis, a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD) is presented. The algorithm targets both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of the GPU’s memory hierarchy. To this end, a family of parallel pivot strategies on the GPU’s shared address space has been developed, but the strategies are applicable to inter-node communication as well, with GPU nodes, CPU nodes, or, in general, any NUMA nodes. Unlike common hybrid approaches, the presented algorithm in a single-GPU setting needs a CPU for the controlling purposes only, while utilizing the GPU’s resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a four-GPU system vs. a single GPU. The subsequent part of the thesis describes how to modify the two-sided Hari–Zimmermann algorithm for computation of the generalized eigendecomposition of a symmetric matrix pair (A; B), where B is positive definite, to an implicit algorithm that computes the generalized singular value decomposition (GSVD) of a pair (F; G). In addition, blocking and parallelization techniques for accelerating both the CPU and the GPU computation are presented, with the GPU approach following the Jacobi SVD algorithm from the first part of the thesis. For triangular matrix pairs of a moderate size, numerical tests show that the double precision sequential pointwise algorithm is several times faster than the established DTGSJA algorithm in LAPACK, while the accuracy is slightly better, especially for the small generalized singular values. Cache-aware blocking increases the performance even further. As with the one-sided Jacobi-type (G)SVD algorithms in general, the presented algorithm is almost perfectly parallelizable and scalable on the shared memory machines, where the speedup almost solely depends on the number of cores used. A distributed memory variant, intended for huge matrices that do not fit into a single NUMA node, as well as a GPU variant, are also sketched. The thesis concludes with the affirmative answer to a question whether the onesided Jacobi-type algorithms can be an efficient and scalable choice for computing the (G)SVD of dense matrices on the massively parallel CPU and GPU architectures. Unless otherwise noted by the inline citations or implied by the context, this thesis is an overview of the original research results, most of which has already been published in [55, 58]. The author’s contributions are the one-sided Jacobi-type GPU algorithms for the ordinary and the generalized SVD, of which the latter has not yet been published, as well as the parallelization technique and some implementation details of the one-sided Hari–Zimmermann CPU algorithm for the GSVD. The rest is joint work with Sanja and Saša Singer.Singularna dekompozicija, katkad zvana prema engleskom originalu i dekompozicija singularnih vrijednosti, ili kraće SVD, jedna je od najkorisnijih matričnih dekompozicija, kako za teorijske, tako i za praktične svrhe. Svaka matrica $G \in \mathbb{C}^{m \times n}$ (zbog jednostavnijeg zapisa, uobičajeno se smatra da je $m \geq n$; u protivnom, traži se SVD matrice $G^\ast$) može se rastaviti u produkt tri matrice $$G = U \Sigma V^\ast,$$ gdje su $U \in \mathbb{C}^{m \times m}$ i $V \in \mathbb{C}^{n \times n}$ unitarne, a $\Sigma \in \mathbb{R}^{m \times n}$ je 'dijagonalna' s nenegativnim dijagonalnim elementima. Osim ovog oblika dekompozicije, koristi se i skraćeni oblik $$G = U'\Sigma'V^\ast,$$ pri čemu je $U' \in \mathbb{C}^{m \times n}$ matrica s ortonormiranim stupcima, a $\Sigma' = diag(\sigma_1, \dots, \sigma_n), \sigma_i \geq 0$ za $i = 0, \dots, n$, je sada stvarno dijagonalna. Izvan matematike, u 'stvarnom' životu, SVD se koristi u procesiranju slika (rekonstrukciji, sažimanju, izoštravanju) i signala, s primjenama u medicini (CT, tj. kompjuterizirana tomografija; MR, tj. magnetna rezonancija), geoznanostima, znanosti o materijalima, kristalografiji, sigurnosti (prepoznavanje lica), izvlačenja informacija iz velike količine podataka (na primjer, LSI, tj. latent semantic indexing), ali i drugdje. Većina primjena koristi svojstvo da se iz SVD-a lako čita najbolja aproksimacija dane matrice matricom fiksnog (niskog) ranga. Čini se da je lakše reći gdje se SVD ne koristi, nego gdje se koristi, stoga se SVD često naziva i "švicarskim nožićem matričnih dekompozicija"$^1$. Prvi počeci razvoja SVD-a sežu u 19. stoljeće, kad su poznati matematičari Eugenio Beltrami, Camille Jordan, James Joseph Sylvester, Erhard Schmidt i Herman Weyl pokazali njezinu egzistenciju i osnovna svojstva (za detalje pogledati [74]). Pioniri u numeričkom računanju SVD-a su Ervand George Kogbetliantz, te Gene Golub i William Kahan, koji su razvili algoritam za računanje (bidijagonalni QR), koji je dvadeset i pet godina vladao scenom numeričkog računanja SVD-a. U to vrijeme, sveučilište Stanford (gdje je Gene Golub radio) bilo je 'glavno sjedište' za razvoj primjena SVD-a. Početkom devedesetih godina, 'sjedište SVD-a' preseljeno je u Europu, nakon objave članka [21] o relativnoj točnosti računanja svojstvenih vrijednosti simetričnih pozitivno definitnih matrica korištenjem Jacobijeve metode. Naime, problem računanja svojstvene dekompozicije pozitivno definitne matrice i problem računanja SVD-a usko su vezani. Ako je poznata dekompozicija singularnih vrijednosti matrice $G$ punog stupčanog ranga, $G \in \mathbb{C}^{m \times n} = U \Sigma V^\ast$, pri čemu je $G$ faktor matrice $A$, $A = G \ast G$, onda je $A$ simetrična i pozitivno definitna i vrijedi $$A = G \ast G = V \Sigma^T U^\ast U \Sigma V^\ast = V diag(\sigma_1^2, \dots, \sigma_m^2)V^\ast .$$ Matrica $V$ je matrica svojstvenih vektora, a svojstvene vrijednosti su kvadrati singularnih vrijednosti. Stoga se algoritmi za računanje svojstvenih vrijednosti, kod kojih se transformacija vrši dvostranim (i slijeva i zdesna) djelovanjem na matricu $A$, mogu napisati implicitno, tako da se transformacija vrši ili zdesna na faktor $G$ ili slijeva na faktor $G^\ast$. U svojoj doktorskoj disertaciji Drmač [24] je napravio daljnju analizu, ne samo singularne dekompozicije računate Jacobijevim algoritmom, nego i generalizirane singularne dekompozicije (GSVD). Temeljem tih istraživanja, SVD baziran na Jacobijevim rotacijama ušao je i u numeričku biblioteku LAPACK. U međuvremenu, gotovo sva računala postala su višejezgrena, a moderni klasteri računala za znanstveno računanje sastoje se od nekoliko tisuća do nekoliko stotina tisuća višejezgrenih procesora$^2$, pa standardni sekvencijalni algoritmi nipošto više nisu primjereni za numeričko računanje. Stoga se ubrzano razvijaju paralelni algoritmi koji poštuju i hijerarhijsku memorijsku strukturu odgovarajućih računala, težeći iskoristiti brzu cache memoriju za procesiranje potproblema u blokovima, na koje je moguće primijeniti BLAS-3 operacije. Ideja blokiranja je u primjeni što više (tipično, kubično u dimenziji matrice) numeričkih operacija nad podacima u brzoj memoriji. Nadalje, pojavom grafičkih procesnih jedinica namijenjenih znanstvenom računanju, kao i drugih visokoparalelnih numeričkih akceleratora (npr. Intel Xeon Phi), otvorio se novi segment istraživanja, koji poštuje njihov masivni paralelizam, s pojedinačno slabašnom snagom svake dretve u odnosu na središnji procesor. Generaliziranu singularnu dekompoziciju (GSVD) uveli su Van Loan [77], te Paige i Saunders [62]. Definicija GSVD-a nešto je manje poznata. Ako su zadane matrice $F \in \mathbb{C}^{m \times n}$ i $G \in \mathbb{C}^{p \times n}$, za koje vrijedi $$K = {F \brack G} , k = rank(K),$$ tad postoje unitarne matrice $U \in \mathbb{C}^{m \times m}, V \in \mathbb{C}^{p \times p}$, i matrica $X \in \mathbb{C}^{k \times n}$, takve da je $$F = U \Sigma_F X, \qquad G = V \Sigma_G X, \qquad \Sigma_F \in \mathbb{R}^{m \times k}, \qquad \Sigma_G \in \mathbb{R}^{p \times k}.$$ Elementi matrica $\Sigma_F$ i $\Sigma_G$ su nula, osim dijagonalnih elemenata, koji su realni i nenegativni. Nadalje, $\Sigma_F$ i $\Sigma_G$ zadovoljavaju $$\Sigma_F^T\Sigma_F + \Sigma_G^T\Sigma_G = I.$$ Omjeri $(\Sigma_F)_{ii} / (\Sigma_G)_{ii}$ su generalizirane singularne vrijednosti para $(F, G)$. Ako je $G$ punog stupčanog ranga, tada je $rank(K) = n$ i generalizirane singularne vrijednosti su konačni brojevi. Ako je par $(F, G)$ realan, onda su realne sve matrice u dekompoziciji. Odavde nadalje, zbog jednostavnoti pretpostavlja se da je par realan. Može se pokazati da, ako je $k = n$, tada se relacija između GSVD-a i reducirane forme CS (kosinus-sinus) dekompozicije (vidjeti, na primjer, [26]) može iskoristiti za njezino računanje (pogledati, na primjer članke Stewarta [72, 73] i Suttona [75]). Slično kao i SVD, generalizirana singularna dekompozicija ima primjene u mnogim područjima, kao što je usporedna analiza podataka vezanih uz genome [1], nepotpuna singularna metoda rubnih elemeneata [47], ionosferna tomografija [9], ali i mnogo drugih. GSVD para matrica $(F, G)$ blisko je vezana s hermitskim generaliziranim svojstvenim problemom za par $(A, B) := (F^\ast F, G^\ast G)$, tako da se metode za istovremenu dijagonalizaciju para $(A, B)$ mogu modificirati za računanje GSVD-a para $(F, G)$. U ovoj radnji razvijen je brzi i efikasan algoritam za računanje generalizirane singularne dekompozicije realnog para $(F, G)$. Metoda razvijena u radnji bazirana je na algoritmu za računanje generalizirane svojstvene dekompozicije, $$Ax = \lambda Bx; \quad x \neq 0; \qquad (1)$$ gdje su $A$ i $B$ simetrične matrice, a par je definitan, tj. postoji realna konstanta $\mu$ takva da je matrica $A-\mu B$ pozitivno definitna. Članke s metodom objavili su 1960. Falk i Langemeyer [31, 32] u slabo poznatom priručniku. Kad je paralelna verzija metode testirana, pokazalo se da pati zbog problema rastuće skale stupaca matrice tijekom procesa ortogonalizacije. Treba još primijetiti da pozitivna definitnost matrice $B$ odmah znači da je definitan i par $(A, B)$. Gotovo desetljeće nakon Falka i Langemeyera, Katharina Zimmermann je u svojoj doktorskoj disertaciji [81] grubo skicirala metodu za rješavanje generaliziranog svojstvenog problema (1) ako je B pozitivno definitna. Gose [34] je predložio optimalnu ne-cikličku pivotnu strategiju i dokazao globalnu konvergenciju originalne metode. Hari je u svojoj disertaciji [37], potaknut Zimmermanninom skicom metode, izveo algoritam i pokazao njegovu globalnu i kvadratičnu konvergenciju uz cikličke pivotne strategije. Kvadratičnu konvergenciju originalne Falk–Langemeyerove metode dokazao je 1988. Slapničar u svojem magisteriju, četiri godine nakon dokaza konvergencije Hari–Zimmermann metode. Hari je u [37] pokazao ključnu vezu između Hari–Zimmermannine i Falk–Langemeyerove varijante algoritma. Ako je matrica $B$ obostrano skalirana dijagonalnom matricom $D$, tako da su joj dijagonalni elementi jednaki 1 prije svakog koraka poništavanja u Falk–Langemeyerovoj metodi, dobiva se Hari–Zimmermannina metoda. Dakle, nova metoda imala je ključno svojstvo normiranosti stupaca barem jedne matrice, što se pokazalo iznimno bitnim za uspjeh algoritma (izbjegavanje skaliranja matrica tijekom procesa ortogonalizacije). Treba reći da se GSVD može računati i na druge načine. Drmač je u [26] izveo algoritam za računanje GSVD-a para $(F, G)$, kad je $G$ punog stupčanog ranga. Algoritam transformira problem na samo jednu matricu, a nakon toga primjenjuje jednostrani Jacobijev SVD algoritam. Taj algoritam računa generalizirane singularne vrijednosti s malom relativnom greškom. Algoritam svođenja na jednu matricu sastoji se od tri koraka: skaliranje stupaca matrica $F$ i $G$, QR faktorizacije sa stupčanim pivotiranjem već skalirane matrice $G$, i konačno, rješavanjem trokutastog linearnog sustava s $k$ desnih strana. Posljednja dva koraka su sekvencijalna i vrlo ih je teško paralelizirati. Sama ideja korištenja implicitne (tj. jednostrane) Falk–Langemeyerove metode za GSVD para $(F, G)$, s $G$ punog stupčanog ranga, sreće se u disertaciji Annette Deichmöller [17], međutim, tamo se ne spominju usporedbe te metode s drugim metodama. S druge strane, algoritam za računanje GSVD-a u biblioteci LAPACK (potprogram xGGSVD), je modificirani Kogbetliantzov algoritam (vidjeti Paige [61]) s obveznim pretprocesiranjem (vidjeti Bai i Demmel [5]). Algoritam pretprocesiranja [6] transformira zadani matrični par $(F_0, G_0)$ u par $(F, G)$, takav da su $F$ i $G$ gornjetrokutaste, a $G$ je i nesingularna. Ako se unaprijed zna da je $G$ punog stupčanog ranga, i implicitna Falk–Langemeyerova i implicitna Hari–Zimmermannina metoda će raditi i bez pretprocesiranja. Ako su $F$ i $G$ vitke (engl. "tall and skinny"), QR factorizacija obje matrice će ubrzati ortogonalizaciju. Ako $G$ nije punog ranga, onda treba koristiti isto pretprocesiranje kao u LAPACK-u, budući da puni stupčani rang matrice $G$ garantira pozitivnu definitnost matrice $B := G^T G$. U ovoj radnji razvijen je i hijerarhijski, blokirani jednostrani algoritam za računanje SVD-a. Opisani algoritam može raditi na višeprocesorskom računalu, računalnim klasterima, jednoj ili više grafičkih procesnih jedinica. Princip rada algoritma na svim arhitekturama je sličan. Posebno je opisan algoritam koji radi na grafičkim procesnim jedinicama. Struktura blokiranja reflektira razine memorijske strukture grafičke procesne jedninice. Da bi se to postiglo, razvijene su familije paralelnih pivotnih strategija za dijeljenu (engl. shared) memoriju grafičkih procesnih jedinica. Uz dodatak rasporeda po procesima, strategije se mogu koristiti i kao strategije za komuniciranje među računalnim čvorovima (bili oni grafičke procesne jedinice, jezgre procesora ili tzv. NUMA čvorovi). Razvijeni algoritam nije hibridni, tj. centralnu procesnu jedinicu koristi samo za kontrolne svrhe, a cjelokupno računanje odvija se na grafičkoj procesnoj jedinici. Kad je zbog veličine problema potrebno, algoritam se može rasprostrijeti (skalirati) na proizvoljan broj grafičkih procesnih jedinica. Na dovoljno velikim matricama, skalabilnost je pokazana ubrzanjem od preko dva puta na četiri grafičke procesne jedinice, obzirom na jednu. U drugom dijelu radnje opisuje se jedan način modifikacije dvostranog Hari–Zimmermanninog algoritma za računanje generalizirane svojstvene dekompozicije matričnog para $(A, B)$, gdje su obje matrice simetrične, a $B$ je pozitivno definitna. Implicitni algoritam računa GSVD para $(F, G)$, pri čemu je $(A, B) := (F^T F, G^T G)$. Nadalje, pokazuje se kako treba blokirati algoritam, te kako ga paralelizirati, i u slučaju standardnih, i u slučaju grafičkih procesora. Za trokutaste matrične parove srednje velikih dimenzija (približno 5 000), pokazano je da je već sekvencijalni, neblokirani algoritam u dvostrukoj točnosti, predložen u radnji, nekoliko desetaka puta brži no što je to LAPACK potprogram DTGSJA i pritom ima nešto bolju točnost, posebno za male generalizirane singularne vrijednosti. Blokiranje algoritma koje odgovara cacheima znatno ubrzava algoritam. Pokazuje se da je i ovaj algoritam, slično kao jednostrani Jacobijev algoritam za SVD, gotovo idealno paralelizabilan i skalabilan na računalima s dijeljenom memorijom, te da njegovo ubrzanje gotovo isključivo ovisi o broju korištenih jezgara. U vrijeme testiranja, pokazalo se da je paralelizirani i blokirani Hari–Zimmermannin algoritam preko sto puta brži od LAPACK potprograma DTGESJA s višedretvenim BLAS potprogramima. Varijanta algoritma za razdijeljenu (engl. distributed) memoriju namijenjena je ogromnim matricama koje ne stanu u jedan NUMA čvor. Također, skicirana je i GPU varijanta algoritma, koja je vrlo slična jednostranom Jacobijevom algoritmu za SVD. Disertacija završava zaključkom da su ovi algoritmi Jacobijevog tipa efikasni i skalabilni i izvrstan su izbor za računanje (G)SVD-a punih matrica na masivno paralelnim standardnim arhitekturama i na grafičkim procesnim jedinicama. Ova doktorska disertacija bazirana je na originalnim znanstvenim radovima [55, 58], te proširena nekim novim rezultatima. Autorov doprinos u ovoj disertaciji su novi paralelni algoritmi za (G)SVD za grafičke procesne jedinice, tehnike paralelizacije, te detalji implementacije jednostranog Hari–Zimmermannina algoritma. Ostatak je zajednički rad sa Sanjom Singer i Sašom Singerom. $^1$ Diane O’Leary, 2006. $^2$ https://www.top500.or

Paralelni algoritmi Jacobijeva tipa za singularnu i generaliziranu singularnu dekompoziciju

In this thesis, a hierarchically blocked one-sided Jacobi algorithm for the singular value decomposition (SVD) is presented. The algorithm targets both single and multiple graphics processing units (GPUs). The blocking structure reflects the levels of the GPU’s memory hierarchy. To this end, a family of parallel pivot strategies on the GPU’s shared address space has been developed, but the strategies are applicable to inter-node communication as well, with GPU nodes, CPU nodes, or, in general, any NUMA nodes. Unlike common hybrid approaches, the presented algorithm in a single-GPU setting needs a CPU for the controlling purposes only, while utilizing the GPU’s resources to the fullest extent permitted by the hardware. When required by the problem size, the algorithm, in principle, scales to an arbitrary number of GPU nodes. The scalability is demonstrated by more than twofold speedup for sufficiently large matrices on a four-GPU system vs. a single GPU. The subsequent part of the thesis describes how to modify the two-sided Hari–Zimmermann algorithm for computation of the generalized eigendecomposition of a symmetric matrix pair (A; B), where B is positive definite, to an implicit algorithm that computes the generalized singular value decomposition (GSVD) of a pair (F; G). In addition, blocking and parallelization techniques for accelerating both the CPU and the GPU computation are presented, with the GPU approach following the Jacobi SVD algorithm from the first part of the thesis. For triangular matrix pairs of a moderate size, numerical tests show that the double precision sequential pointwise algorithm is several times faster than the established DTGSJA algorithm in LAPACK, while the accuracy is slightly better, especially for the small generalized singular values. Cache-aware blocking increases the performance even further. As with the one-sided Jacobi-type (G)SVD algorithms in general, the presented algorithm is almost perfectly parallelizable and scalable on the shared memory machines, where the speedup almost solely depends on the number of cores used. A distributed memory variant, intended for huge matrices that do not fit into a single NUMA node, as well as a GPU variant, are also sketched. The thesis concludes with the affirmative answer to a question whether the onesided Jacobi-type algorithms can be an efficient and scalable choice for computing the (G)SVD of dense matrices on the massively parallel CPU and GPU architectures. Unless otherwise noted by the inline citations or implied by the context, this thesis is an overview of the original research results, most of which has already been published in [55, 58]. The author’s contributions are the one-sided Jacobi-type GPU algorithms for the ordinary and the generalized SVD, of which the latter has not yet been published, as well as the parallelization technique and some implementation details of the one-sided Hari–Zimmermann CPU algorithm for the GSVD. The rest is joint work with Sanja and Saša Singer.Singularna dekompozicija, katkad zvana prema engleskom originalu i dekompozicija singularnih vrijednosti, ili kraće SVD, jedna je od najkorisnijih matričnih dekompozicija, kako za teorijske, tako i za praktične svrhe. Svaka matrica $G \in \mathbb{C}^{m \times n}$ (zbog jednostavnijeg zapisa, uobičajeno se smatra da je $m \geq n$; u protivnom, traži se SVD matrice $G^\ast$) može se rastaviti u produkt tri matrice $$G = U \Sigma V^\ast,$$ gdje su $U \in \mathbb{C}^{m \times m}$ i $V \in \mathbb{C}^{n \times n}$ unitarne, a $\Sigma \in \mathbb{R}^{m \times n}$ je 'dijagonalna' s nenegativnim dijagonalnim elementima. Osim ovog oblika dekompozicije, koristi se i skraćeni oblik $$G = U'\Sigma'V^\ast,$$ pri čemu je $U' \in \mathbb{C}^{m \times n}$ matrica s ortonormiranim stupcima, a $\Sigma' = diag(\sigma_1, \dots, \sigma_n), \sigma_i \geq 0$ za $i = 0, \dots, n$, je sada stvarno dijagonalna. Izvan matematike, u 'stvarnom' životu, SVD se koristi u procesiranju slika (rekonstrukciji, sažimanju, izoštravanju) i signala, s primjenama u medicini (CT, tj. kompjuterizirana tomografija; MR, tj. magnetna rezonancija), geoznanostima, znanosti o materijalima, kristalografiji, sigurnosti (prepoznavanje lica), izvlačenja informacija iz velike količine podataka (na primjer, LSI, tj. latent semantic indexing), ali i drugdje. Većina primjena koristi svojstvo da se iz SVD-a lako čita najbolja aproksimacija dane matrice matricom fiksnog (niskog) ranga. Čini se da je lakše reći gdje se SVD ne koristi, nego gdje se koristi, stoga se SVD često naziva i "švicarskim nožićem matričnih dekompozicija"$^1$. Prvi počeci razvoja SVD-a sežu u 19. stoljeće, kad su poznati matematičari Eugenio Beltrami, Camille Jordan, James Joseph Sylvester, Erhard Schmidt i Herman Weyl pokazali njezinu egzistenciju i osnovna svojstva (za detalje pogledati [74]). Pioniri u numeričkom računanju SVD-a su Ervand George Kogbetliantz, te Gene Golub i William Kahan, koji su razvili algoritam za računanje (bidijagonalni QR), koji je dvadeset i pet godina vladao scenom numeričkog računanja SVD-a. U to vrijeme, sveučilište Stanford (gdje je Gene Golub radio) bilo je 'glavno sjedište' za razvoj primjena SVD-a. Početkom devedesetih godina, 'sjedište SVD-a' preseljeno je u Europu, nakon objave članka [21] o relativnoj točnosti računanja svojstvenih vrijednosti simetričnih pozitivno definitnih matrica korištenjem Jacobijeve metode. Naime, problem računanja svojstvene dekompozicije pozitivno definitne matrice i problem računanja SVD-a usko su vezani. Ako je poznata dekompozicija singularnih vrijednosti matrice $G$ punog stupčanog ranga, $G \in \mathbb{C}^{m \times n} = U \Sigma V^\ast$, pri čemu je $G$ faktor matrice $A$, $A = G \ast G$, onda je $A$ simetrična i pozitivno definitna i vrijedi $$A = G \ast G = V \Sigma^T U^\ast U \Sigma V^\ast = V diag(\sigma_1^2, \dots, \sigma_m^2)V^\ast .$$ Matrica $V$ je matrica svojstvenih vektora, a svojstvene vrijednosti su kvadrati singularnih vrijednosti. Stoga se algoritmi za računanje svojstvenih vrijednosti, kod kojih se transformacija vrši dvostranim (i slijeva i zdesna) djelovanjem na matricu $A$, mogu napisati implicitno, tako da se transformacija vrši ili zdesna na faktor $G$ ili slijeva na faktor $G^\ast$. U svojoj doktorskoj disertaciji Drmač [24] je napravio daljnju analizu, ne samo singularne dekompozicije računate Jacobijevim algoritmom, nego i generalizirane singularne dekompozicije (GSVD). Temeljem tih istraživanja, SVD baziran na Jacobijevim rotacijama ušao je i u numeričku biblioteku LAPACK. U međuvremenu, gotovo sva računala postala su višejezgrena, a moderni klasteri računala za znanstveno računanje sastoje se od nekoliko tisuća do nekoliko stotina tisuća višejezgrenih procesora$^2$, pa standardni sekvencijalni algoritmi nipošto više nisu primjereni za numeričko računanje. Stoga se ubrzano razvijaju paralelni algoritmi koji poštuju i hijerarhijsku memorijsku strukturu odgovarajućih računala, težeći iskoristiti brzu cache memoriju za procesiranje potproblema u blokovima, na koje je moguće primijeniti BLAS-3 operacije. Ideja blokiranja je u primjeni što više (tipično, kubično u dimenziji matrice) numeričkih operacija nad podacima u brzoj memoriji. Nadalje, pojavom grafičkih procesnih jedinica namijenjenih znanstvenom računanju, kao i drugih visokoparalelnih numeričkih akceleratora (npr. Intel Xeon Phi), otvorio se novi segment istraživanja, koji poštuje njihov masivni paralelizam, s pojedinačno slabašnom snagom svake dretve u odnosu na središnji procesor. Generaliziranu singularnu dekompoziciju (GSVD) uveli su Van Loan [77], te Paige i Saunders [62]. Definicija GSVD-a nešto je manje poznata. Ako su zadane matrice $F \in \mathbb{C}^{m \times n}$ i $G \in \mathbb{C}^{p \times n}$, za koje vrijedi $$K = {F \brack G} , k = rank(K),$$ tad postoje unitarne matrice $U \in \mathbb{C}^{m \times m}, V \in \mathbb{C}^{p \times p}$, i matrica $X \in \mathbb{C}^{k \times n}$, takve da je $$F = U \Sigma_F X, \qquad G = V \Sigma_G X, \qquad \Sigma_F \in \mathbb{R}^{m \times k}, \qquad \Sigma_G \in \mathbb{R}^{p \times k}.$$ Elementi matrica $\Sigma_F$ i $\Sigma_G$ su nula, osim dijagonalnih elemenata, koji su realni i nenegativni. Nadalje, $\Sigma_F$ i $\Sigma_G$ zadovoljavaju $$\Sigma_F^T\Sigma_F + \Sigma_G^T\Sigma_G = I.$$ Omjeri $(\Sigma_F)_{ii} / (\Sigma_G)_{ii}$ su generalizirane singularne vrijednosti para $(F, G)$. Ako je $G$ punog stupčanog ranga, tada je $rank(K) = n$ i generalizirane singularne vrijednosti su konačni brojevi. Ako je par $(F, G)$ realan, onda su realne sve matrice u dekompoziciji. Odavde nadalje, zbog jednostavnoti pretpostavlja se da je par realan. Može se pokazati da, ako je $k = n$, tada se relacija između GSVD-a i reducirane forme CS (kosinus-sinus) dekompozicije (vidjeti, na primjer, [26]) može iskoristiti za njezino računanje (pogledati, na primjer članke Stewarta [72, 73] i Suttona [75]). Slično kao i SVD, generalizirana singularna dekompozicija ima primjene u mnogim područjima, kao što je usporedna analiza podataka vezanih uz genome [1], nepotpuna singularna metoda rubnih elemeneata [47], ionosferna tomografija [9], ali i mnogo drugih. GSVD para matrica $(F, G)$ blisko je vezana s hermitskim generaliziranim svojstvenim problemom za par $(A, B) := (F^\ast F, G^\ast G)$, tako da se metode za istovremenu dijagonalizaciju para $(A, B)$ mogu modificirati za računanje GSVD-a para $(F, G)$. U ovoj radnji razvijen je brzi i efikasan algoritam za računanje generalizirane singularne dekompozicije realnog para $(F, G)$. Metoda razvijena u radnji bazirana je na algoritmu za računanje generalizirane svojstvene dekompozicije, $$Ax = \lambda Bx; \quad x \neq 0; \qquad (1)$$ gdje su $A$ i $B$ simetrične matrice, a par je definitan, tj. postoji realna konstanta $\mu$ takva da je matrica $A-\mu B$ pozitivno definitna. Članke s metodom objavili su 1960. Falk i Langemeyer [31, 32] u slabo poznatom priručniku. Kad je paralelna verzija metode testirana, pokazalo se da pati zbog problema rastuće skale stupaca matrice tijekom procesa ortogonalizacije. Treba još primijetiti da pozitivna definitnost matrice $B$ odmah znači da je definitan i par $(A, B)$. Gotovo desetljeće nakon Falka i Langemeyera, Katharina Zimmermann je u svojoj doktorskoj disertaciji [81] grubo skicirala metodu za rješavanje generaliziranog svojstvenog problema (1) ako je B pozitivno definitna. Gose [34] je predložio optimalnu ne-cikličku pivotnu strategiju i dokazao globalnu konvergenciju originalne metode. Hari je u svojoj disertaciji [37], potaknut Zimmermanninom skicom metode, izveo algoritam i pokazao njegovu globalnu i kvadratičnu konvergenciju uz cikličke pivotne strategije. Kvadratičnu konvergenciju originalne Falk–Langemeyerove metode dokazao je 1988. Slapničar u svojem magisteriju, četiri godine nakon dokaza konvergencije Hari–Zimmermann metode. Hari je u [37] pokazao ključnu vezu između Hari–Zimmermannine i Falk–Langemeyerove varijante algoritma. Ako je matrica $B$ obostrano skalirana dijagonalnom matricom $D$, tako da su joj dijagonalni elementi jednaki 1 prije svakog koraka poništavanja u Falk–Langemeyerovoj metodi, dobiva se Hari–Zimmermannina metoda. Dakle, nova metoda imala je ključno svojstvo normiranosti stupaca barem jedne matrice, što se pokazalo iznimno bitnim za uspjeh algoritma (izbjegavanje skaliranja matrica tijekom procesa ortogonalizacije). Treba reći da se GSVD može računati i na druge načine. Drmač je u [26] izveo algoritam za računanje GSVD-a para $(F, G)$, kad je $G$ punog stupčanog ranga. Algoritam transformira problem na samo jednu matricu, a nakon toga primjenjuje jednostrani Jacobijev SVD algoritam. Taj algoritam računa generalizirane singularne vrijednosti s malom relativnom greškom. Algoritam svođenja na jednu matricu sastoji se od tri koraka: skaliranje stupaca matrica $F$ i $G$, QR faktorizacije sa stupčanim pivotiranjem već skalirane matrice $G$, i konačno, rješavanjem trokutastog linearnog sustava s $k$ desnih strana. Posljednja dva koraka su sekvencijalna i vrlo ih je teško paralelizirati. Sama ideja korištenja implicitne (tj. jednostrane) Falk–Langemeyerove metode za GSVD para $(F, G)$, s $G$ punog stupčanog ranga, sreće se u disertaciji Annette Deichmöller [17], međutim, tamo se ne spominju usporedbe te metode s drugim metodama. S druge strane, algoritam za računanje GSVD-a u biblioteci LAPACK (potprogram xGGSVD), je modificirani Kogbetliantzov algoritam (vidjeti Paige [61]) s obveznim pretprocesiranjem (vidjeti Bai i Demmel [5]). Algoritam pretprocesiranja [6] transformira zadani matrični par $(F_0, G_0)$ u par $(F, G)$, takav da su $F$ i $G$ gornjetrokutaste, a $G$ je i nesingularna. Ako se unaprijed zna da je $G$ punog stupčanog ranga, i implicitna Falk–Langemeyerova i implicitna Hari–Zimmermannina metoda će raditi i bez pretprocesiranja. Ako su $F$ i $G$ vitke (engl. "tall and skinny"), QR factorizacija obje matrice će ubrzati ortogonalizaciju. Ako $G$ nije punog ranga, onda treba koristiti isto pretprocesiranje kao u LAPACK-u, budući da puni stupčani rang matrice $G$ garantira pozitivnu definitnost matrice $B := G^T G$. U ovoj radnji razvijen je i hijerarhijski, blokirani jednostrani algoritam za računanje SVD-a. Opisani algoritam može raditi na višeprocesorskom računalu, računalnim klasterima, jednoj ili više grafičkih procesnih jedinica. Princip rada algoritma na svim arhitekturama je sličan. Posebno je opisan algoritam koji radi na grafičkim procesnim jedinicama. Struktura blokiranja reflektira razine memorijske strukture grafičke procesne jedninice. Da bi se to postiglo, razvijene su familije paralelnih pivotnih strategija za dijeljenu (engl. shared) memoriju grafičkih procesnih jedinica. Uz dodatak rasporeda po procesima, strategije se mogu koristiti i kao strategije za komuniciranje među računalnim čvorovima (bili oni grafičke procesne jedinice, jezgre procesora ili tzv. NUMA čvorovi). Razvijeni algoritam nije hibridni, tj. centralnu procesnu jedinicu koristi samo za kontrolne svrhe, a cjelokupno računanje odvija se na grafičkoj procesnoj jedinici. Kad je zbog veličine problema potrebno, algoritam se može rasprostrijeti (skalirati) na proizvoljan broj grafičkih procesnih jedinica. Na dovoljno velikim matricama, skalabilnost je pokazana ubrzanjem od preko dva puta na četiri grafičke procesne jedinice, obzirom na jednu. U drugom dijelu radnje opisuje se jedan način modifikacije dvostranog Hari–Zimmermanninog algoritma za računanje generalizirane svojstvene dekompozicije matričnog para $(A, B)$, gdje su obje matrice simetrične, a $B$ je pozitivno definitna. Implicitni algoritam računa GSVD para $(F, G)$, pri čemu je $(A, B) := (F^T F, G^T G)$. Nadalje, pokazuje se kako treba blokirati algoritam, te kako ga paralelizirati, i u slučaju standardnih, i u slučaju grafičkih procesora. Za trokutaste matrične parove srednje velikih dimenzija (približno 5 000), pokazano je da je već sekvencijalni, neblokirani algoritam u dvostrukoj točnosti, predložen u radnji, nekoliko desetaka puta brži no što je to LAPACK potprogram DTGSJA i pritom ima nešto bolju točnost, posebno za male generalizirane singularne vrijednosti. Blokiranje algoritma koje odgovara cacheima znatno ubrzava algoritam. Pokazuje se da je i ovaj algoritam, slično kao jednostrani Jacobijev algoritam za SVD, gotovo idealno paralelizabilan i skalabilan na računalima s dijeljenom memorijom, te da njegovo ubrzanje gotovo isključivo ovisi o broju korištenih jezgara. U vrijeme testiranja, pokazalo se da je paralelizirani i blokirani Hari–Zimmermannin algoritam preko sto puta brži od LAPACK potprograma DTGESJA s višedretvenim BLAS potprogramima. Varijanta algoritma za razdijeljenu (engl. distributed) memoriju namijenjena je ogromnim matricama koje ne stanu u jedan NUMA čvor. Također, skicirana je i GPU varijanta algoritma, koja je vrlo slična jednostranom Jacobijevom algoritmu za SVD. Disertacija završava zaključkom da su ovi algoritmi Jacobijevog tipa efikasni i skalabilni i izvrstan su izbor za računanje (G)SVD-a punih matrica na masivno paralelnim standardnim arhitekturama i na grafičkim procesnim jedinicama. Ova doktorska disertacija bazirana je na originalnim znanstvenim radovima [55, 58], te proširena nekim novim rezultatima. Autorov doprinos u ovoj disertaciji su novi paralelni algoritmi za (G)SVD za grafičke procesne jedinice, tehnike paralelizacije, te detalji implementacije jednostranog Hari–Zimmermannina algoritma. Ostatak je zajednički rad sa Sanjom Singer i Sašom Singerom. $^1$ Diane O’Leary, 2006. $^2$ https://www.top500.or

Uncertainty quantification in ocean state estimation

Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy at the Massachusetts Institute of Technology and the Woods Hole Oceanographic Institution February 2013Quantifying uncertainty and error bounds is a key outstanding challenge in ocean state estimation and climate research. It is particularly difficult due to the large dimensionality of this nonlinear estimation problem and the number of uncertain variables involved. The “Estimating the Circulation and Climate of the Oceans” (ECCO) consortium has developed a scalable system for dynamically consistent estimation of global time-evolving ocean state by optimal combination of ocean general circulation model (GCM) with diverse ocean observations. The estimation system is based on the "adjoint method" solution of an unconstrained least-squares optimization problem formulated with the method of Lagrange multipliers for fitting the dynamical ocean model to observations. The dynamical consistency requirement of ocean state estimation necessitates this approach over sequential data assimilation and reanalysis smoothing techniques. In addition, it is computationally advantageous because calculation and storage of large covariance matrices is not required. However, this is also a drawback of the adjoint method, which lacks a native formalism for error propagation and quantification of assimilated uncertainty. The objective of this dissertation is to resolve that limitation by developing a feasible computational methodology for uncertainty analysis in dynamically consistent state estimation, applicable to the large dimensionality of global ocean models. Hessian (second derivative-based) methodology is developed for Uncertainty Quantification (UQ) in large-scale ocean state estimation, extending the gradient-based adjoint method to employ the second order geometry information of the model-data misfit function in a high-dimensional control space. Large error covariance matrices are evaluated by inverting the Hessian matrix with the developed scalable matrix-free numerical linear algebra algorithms. Hessian-vector product and Jacobian derivative codes of the MIT general circulation model (MITgcm) are generated by means of algorithmic differentiation (AD). Computational complexity of the Hessian code is reduced by tangent linear differentiation of the adjoint code, which preserves the speedup of adjoint checkpointing schemes in the second derivative calculation. A Lanczos algorithm is applied for extracting the leading rank eigenvectors and eigenvalues of the Hessian matrix. The eigenvectors represent the constrained uncertainty patterns. The inverse eigenvalues are the corresponding uncertainties. The dimensionality of UQ calculations is reduced by eliminating the uncertainty null-space unconstrained by the supplied observations. Inverse and forward uncertainty propagation schemes are designed for assimilating observation and control variable uncertainties, and for projecting these uncertainties onto oceanographic target quantities. Two versions of these schemes are developed: one evaluates reduction of prior uncertainties, while another does not require prior assumptions. The analysis of uncertainty propagation in the ocean model is time-resolving. It captures the dynamics of uncertainty evolution and reveals transient and stationary uncertainty regimes. The system is applied to quantifying uncertainties of Antarctic Circumpolar Current (ACC) transport in a global barotropic configuration of the MITgcm. The model is constrained by synthetic observations of sea surface height and velocities. The control space consists of two-dimensional maps of initial and boundary conditions and model parameters. The size of the Hessian matrix is O(1010) elements, which would require O(60GB) of uncompressed storage. It is demonstrated how the choice of observations and their geographic coverage determines the reduction in uncertainties of the estimated transport. The system also yields information on how well the control fields are constrained by the observations. The effects of controls uncertainty reduction due to decrease of diagonal covariance terms are compared to dynamical coupling of controls through off-diagonal covariance terms. The correlations of controls introduced by observation uncertainty assimilation are found to dominate the reduction of uncertainty of transport. An idealized analytical model of ACC guides a detailed time-resolving understanding of uncertainty dynamics.This thesis was supported in part by the National Science Foundation (NSF) Collaboration in Mathematical Geosciences (CMG) grant ARC-0934404, and the Department of Energy (DOE) ISICLES initiative under LANL sub-contract 139843-1. Partial funding was provided by the department of Mechanical Engineering at MIT and by the Academic Programs Office at WHOI. My participation in the IMA "Large-scale Inverse Problems and Quantification of Uncertainty" workshop was partially funded by IMA NSF grants

Improving multifrontal solvers by means of algebraic Block Low-Rank representations

We consider the solution of large sparse linear systems by means of direct factorization based on a multifrontal approach. Although numerically robust and easy to use (it only needs algebraic information: the input matrix A and a right-hand side b, even if it can also digest preprocessing strategies based on geometric information), direct factorization methods are computationally intensive both in terms of memory and operations, which limits their scope on very large problems (matrices with up to few hundred millions of equations). This work focuses on exploiting low-rank approximations on multifrontal based direct methods to reduce both the memory footprints and the operation count, in sequential and distributed-memory environments, on a wide class of problems. We first survey the low-rank formats which have been previously developed to efficiently represent dense matrices and have been widely used to design fast solutions of partial differential equations, integral equations and eigenvalue problems. These formats are hierarchical (H and Hierarchically Semiseparable matrices are the most common ones) and have been (both theoretically and practically) shown to substantially decrease the memory and operation requirements for linear algebra computations. However, they impose many structural constraints which can limit their scope and efficiency, especially in the context of general purpose multifrontal solvers. We propose a flat format called Block Low-Rank (BLR) based on a natural blocking of the matrices and explain why it provides all the flexibility needed by a general purpose multifrontal solver in terms of numerical pivoting for stability and parallelism. We compare BLR format with other formats and show that BLR does not compromise much the memory and operation improvements achieved through low-rank approximations. A stability study shows that the approximations are well controlled by an explicit numerical parameter called low-rank threshold, which is critical in order to solve the sparse linear system accurately. Details on how Block Low-Rank factorizations can be efficiently implemented within multifrontal solvers are then given. We propose several Block Low-Rank factorization algorithms which allow for different types of gains. The proposed algorithms have been implemented within the MUMPS (MUltifrontal Massively Parallel Solver) solver. We first report experiments on standard partial differential equations based problems to analyse the main features of our BLR algorithms and to show the potential and flexibility of the approach; a comparison with a Hierarchically SemiSeparable code is also given. Then, Block Low-Rank formats are experimented on large (up to a hundred millions of unknowns) and various problems coming from several industrial applications. We finally illustrate the use of our approach as a preconditioning method for the Conjugate Gradient

Seventh Copper Mountain Conference on Multigrid Methods

The Seventh Copper Mountain Conference on Multigrid Methods was held on 2-7 Apr. 1995 at Copper Mountain, Colorado. This book is a collection of many of the papers presented at the conference and so represents the conference proceedings. NASA Langley graciously provided printing of this document so that all of the papers could be presented in a single forum. Each paper was reviewed by a member of the conference organizing committee under the coordination of the editors. The multigrid discipline continues to expand and mature, as is evident from these proceedings. The vibrancy in this field is amply expressed in these important papers, and the collection shows its rapid trend to further diversity and depth