Non-orthogonal joint block diagonalization based on the LU or QR factorizations for convolutive blind source separation

This article addresses the problem of blind source separation, in which the source signals are most often of the convolutive mixtures, and moreover, the source signals cannot satisfy independent identical distribution generally. One kind of prevailing and representative approaches for overcoming these difficulties is joint block diagonalization (JBD) method. To improve present JBD methods, we present a class of simple Jacobi-type JBD algorithms based on the LU or QR factorizations. Using Jacobi-type matrices we can replace high dimensional minimization problems with a sequence of simple one-dimensional problems. The novel methods are more general i.e. the orthogonal, positive definite or symmetric matrices and a preliminary whitening stage is no more compulsorily required, and further, the convergence is also guaranteed. The performance of the proposed algorithms, compared with the existing state-of-the-art JBD algorithms, is evaluated with computer simulations and vibration experimental. The results of numerical examples demonstrate that the robustness and effectiveness of the two novel algorithms provide a significant improvement i.e., yield less convergence time, higher precision of convergence, better success rate of block diagonalization. And the proposed algorithms are effective in separating the vibration signals of convolutive mixtures

A Kogbetliantz-type algorithm for the hyperbolic SVD

In this paper a two-sided, parallel Kogbetliantz-type algorithm for the hyperbolic singular value decomposition (HSVD) of real and complex square matrices is developed, with a single assumption that the input matrix, of order $n$, admits such a decomposition into the product of a unitary, a non-negative diagonal, and a $J$-unitary matrix, where $J$ is a given diagonal matrix of positive and negative signs. When $J=\pm I$, the proposed algorithm computes the ordinary SVD. The paper's most important contribution -- a derivation of formulas for the HSVD of $2\times 2$ matrices -- is presented first, followed by the details of their implementation in floating-point arithmetic. Next, the effects of the hyperbolic transformations on the columns of the iteration matrix are discussed. These effects then guide a redesign of the dynamic pivot ordering, being already a well-established pivot strategy for the ordinary Kogbetliantz algorithm, for the general, $n\times n$ HSVD. A heuristic but sound convergence criterion is then proposed, which contributes to high accuracy demonstrated in the numerical testing results. Such a $J$-Kogbetliantz algorithm as presented here is intrinsically slow, but is nevertheless usable for matrices of small orders.Comment: a heavily revised version with 32 pages and 4 figure

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

Analyse en Composantes Indépendantes Multidimensionnelles via des cumulants d’ordres variés

The author deals with the problem of multidimensional independent component analysis (MICA) which is the natural generalization of the ordinary problem of independent component analysis (ICA). First, in order to facilitate the use of higher-order cumulants, we present new formulas for the cumulant matrices of a real random vector from its moment matrices. In addition to the usual matrix operations, these formulas use only the Kronecker product, the vec operator and some commutation matrices. These formulas lend themselves to examine more closely the specific structures of cumulant matrices and provide results on the ranks of these matrices that characterize the dependence between random variables composing the random vector. The main practical interest of our matrix formulas lies in much easier cumulant evaluation and faster computation than the conventional method based on repeated use of the Leonov and Shiryaev formulas. In the second part of this thesis, we show that under the usual assumptions of the independent multidimensional component analysis, contracted cumulant matrices at any statistical order are all block diagonalizable in the same basis. We derive an algorithm for solving MICA by block diagonalization and compare the results obtained to the orders 3-6, between them and with other methods, on several synthetic signals. Simple examples are developed to justify the need to combine different levels to ensure the best separation. We also prove that the easiest case to deal with is the case of mixtures of sources that have different dimensions. In the last part of this thesis we propose a set of methods that operate only the higher- order statistics. Under certain additional assumptions, these methods are shown to completely solve the noisy MICA problem without second-order whitening by joint block diagonalization of a cumulant matrices set coming from statistics of orders strictly higher than four. A comparison with the second-order based whitening MICA methods for the separation of fetal and maternal electrical activities (measured using three electrodes placed on the mother’s abdomen) shows that this new family is better suited to this application : it allow an almost perfect separation of these two contributions.L’auteur s’intéresse au problème de l’analyse en composantes indépendantes multidimensionnelles (ACIM) qui est la généralisation naturelle du problème ordinaire de l’analyse en composantes indépendantes (ACI). Tout d’abord, afin de faciliter l’utilisation des cumulants des ordres supérieurs, nous présentons de nou- velles formules pour le calcul matriciel des matrices de cumulants d’un vecteur aléatoire réel à partir de ses matrices de moments. Outre les opérations matricielles usuelles, ces formules utilisent uniquement le produit de Kronecker, l’opérateur vec et des matrices de commutation. Nous pouvons immédiatement à partir de ces formules examiner de plus près les structures particulières des matrices de cumulants et ainsi donner des résultats sur les rangs de ces matrices qui caractérisent la dépendance entre les variables aléatoires constituant le vecteur aléatoire. L’intérêt pratique principal de nos formules matricielles réside certainement dans une évaluation des cumulants beaucoup plus aisée et rapide qu’avec la méthode usuelle basée sur une utilisation répétée des formules de Leonov et Shiryaev. Dans la deuxième partie de cette thèse, nous montrons que sous les hypothèses usuelles de l’analyse en composantes indépendantes mul- tidimensionnelles, les matrices de cumulants contractées à un ordre statistique quelconque sont toutes bloc-diagonalisables dans la même base. Nous en déduisons des algorithmes de résolution d’ACIM par bloc-diagonalisation conjointe et comparons les résultats obtenus aux ordres 3 à 6, entre eux et avec d’autres méthodes, sur quelques signaux synthétiques. Des exemples simples ont élaborés afin de justifier la nécessité de combiner des ordres différents pour garantir la meilleure séparation. Nous prouvons aussi que le cas le plus simple à traiter est celui de mélanges de sources qui ont différentes dimensions. Dans la dernière partie de cette thèse nous proposons une famille de méthodes qui exploitent uniquement les sta- tistiques d’ordres supérieurs à deux. Sous certaines hypothèses supplémentaires, ces méthodes permettent après un blanchiment d’ordre quatre des observations de résoudre complètement le problème ACIM bruité en bloc diagonalisant conjointement un ensemble de matrices de cumulants issues des statistiques d’ordres supérieurs strictement à quatre. Une comparaison avec les méthodes ACIM à blanchiment d’ordre deux pour la séparation des activités électriques foetale et maternelle (mesurées via trois électrodes placées sur l’abdomen de la mère) montre que cette nouvelle famille est mieux adaptée à cette application : elles permettent une séparation quasi parfaite de ces deux contributions

Advanced Algebraic Concepts for Efficient Multi-Channel Signal Processing

Unsere moderne Gesellschaft ist Zeuge eines fundamentalen Wandels in der Art und Weise wie wir mit Technologie interagieren. Geräte werden zunehmend intelligenter - sie verfügen über mehr und mehr Rechenleistung und häufiger über eigene Kommunikationsschnittstellen. Das beginnt bei einfachen Haushaltsgeräten und reicht über Transportmittel bis zu großen überregionalen Systemen wie etwa dem Stromnetz. Die Erfassung, die Verarbeitung und der Austausch digitaler Informationen gewinnt daher immer mehr an Bedeutung. Die Tatsache, dass ein wachsender Anteil der Geräte heutzutage mobil und deshalb batteriebetrieben ist, begründet den Anspruch, digitale Signalverarbeitungsalgorithmen besonders effizient zu gestalten. Dies kommt auch dem Wunsch nach einer Echtzeitverarbeitung der großen anfallenden Datenmengen zugute. Die vorliegende Arbeit demonstriert Methoden zum Finden effizienter algebraischer Lösungen für eine Vielzahl von Anwendungen mehrkanaliger digitaler Signalverarbeitung. Solche Ansätze liefern nicht immer unbedingt die bestmögliche Lösung, kommen dieser jedoch häufig recht nahe und sind gleichzeitig bedeutend einfacher zu beschreiben und umzusetzen. Die einfache Beschreibungsform ermöglicht eine tiefgehende Analyse ihrer Leistungsfähigkeit, was für den Entwurf eines robusten und zuverlässigen Systems unabdingbar ist. Die Tatsache, dass sie nur gebräuchliche algebraische Hilfsmittel benötigen, erlaubt ihre direkte und zügige Umsetzung und den Test unter realen Bedingungen. Diese Grundidee wird anhand von drei verschiedenen Anwendungsgebieten demonstriert. Zunächst wird ein semi-algebraisches Framework zur Berechnung der kanonisch polyadischen (CP) Zerlegung mehrdimensionaler Signale vorgestellt. Dabei handelt es sich um ein sehr grundlegendes Werkzeug der multilinearen Algebra mit einem breiten Anwendungsspektrum von Mobilkommunikation über Chemie bis zur Bildverarbeitung. Verglichen mit existierenden iterativen Lösungsverfahren bietet das neue Framework die Möglichkeit, den Rechenaufwand und damit die Güte der erzielten Lösung zu steuern. Es ist außerdem weniger anfällig gegen eine schlechte Konditionierung der Ausgangsdaten. Das zweite Gebiet, das in der Arbeit besprochen wird, ist die unterraumbasierte hochauflösende Parameterschätzung für mehrdimensionale Signale, mit Anwendungsgebieten im RADAR, der Modellierung von Wellenausbreitung, oder bildgebenden Verfahren in der Medizin. Es wird gezeigt, dass sich derartige mehrdimensionale Signale mit Tensoren darstellen lassen. Dies erlaubt eine natürlichere Beschreibung und eine bessere Ausnutzung ihrer Struktur als das mit Matrizen möglich ist. Basierend auf dieser Idee entwickeln wir eine tensor-basierte Schätzung des Signalraums, welche genutzt werden kann um beliebige existierende Matrix-basierte Verfahren zu verbessern. Dies wird im Anschluss exemplarisch am Beispiel der ESPRIT-artigen Verfahren gezeigt, für die verbesserte Versionen vorgeschlagen werden, die die mehrdimensionale Struktur der Daten (Tensor-ESPRIT), nichzirkuläre Quellsymbole (NC ESPRIT), sowie beides gleichzeitig (NC Tensor-ESPRIT) ausnutzen. Um die endgültige Schätzgenauigkeit objektiv einschätzen zu können wird dann ein Framework für die analytische Beschreibung der Leistungsfähigkeit beliebiger ESPRIT-artiger Algorithmen diskutiert. Verglichen mit existierenden analytischen Ausdrücken ist unser Ansatz allgemeiner, da keine Annahmen über die statistische Verteilung von Nutzsignal und Rauschen benötigt werden und die Anzahl der zur Verfügung stehenden Schnappschüsse beliebig klein sein kann. Dies führt auf vereinfachte Ausdrücke für den mittleren quadratischen Schätzfehler, die Schlussfolgerungen über die Effizienz der Verfahren unter verschiedenen Bedingungen zulassen. Das dritte Anwendungsgebiet ist der bidirektionale Datenaustausch mit Hilfe von Relay-Stationen. Insbesondere liegt hier der Fokus auf Zwei-Wege-Relaying mit Hilfe von Amplify-and-Forward-Relays mit mehreren Antennen, da dieser Ansatz ein besonders gutes Kosten-Nutzen-Verhältnis verspricht. Es wird gezeigt, dass sich die nötige Kanalkenntnis mit einem einfachen algebraischen Tensor-basierten Schätzverfahren gewinnen lässt. Außerdem werden Verfahren zum Finden einer günstigen Relay-Verstärkungs-Strategie diskutiert. Bestehende Ansätze basieren entweder auf komplexen numerischen Optimierungsverfahren oder auf Ad-Hoc-Ansätzen die keine zufriedenstellende Bitfehlerrate oder Summenrate liefern. Deshalb schlagen wir algebraische Ansätze zum Finden der Relayverstärkungsmatrix vor, die von relevanten Systemmetriken inspiriert sind und doch einfach zu berechnen sind. Wir zeigen das algebraische ANOMAX-Verfahren zum Erreichen einer niedrigen Bitfehlerrate und seine Modifikation RR-ANOMAX zum Erreichen einer hohen Summenrate. Für den Spezialfall, in dem die Endgeräte nur eine Antenne verwenden, leiten wir eine semi-algebraische Lösung zum Finden der Summenraten-optimalen Strategie (RAGES) her. Anhand von numerischen Simulationen wird die Leistungsfähigkeit dieser Verfahren bezüglich Bitfehlerrate und erreichbarer Datenrate bewertet und ihre Effektivität gezeigt.Modern society is undergoing a fundamental change in the way we interact with technology. More and more devices are becoming "smart" by gaining advanced computation capabilities and communication interfaces, from household appliances over transportation systems to large-scale networks like the power grid. Recording, processing, and exchanging digital information is thus becoming increasingly important. As a growing share of devices is nowadays mobile and hence battery-powered, a particular interest in efficient digital signal processing techniques emerges. This thesis contributes to this goal by demonstrating methods for finding efficient algebraic solutions to various applications of multi-channel digital signal processing. These may not always result in the best possible system performance. However, they often come close while being significantly simpler to describe and to implement. The simpler description facilitates a thorough analysis of their performance which is crucial to design robust and reliable systems. The fact that they rely on standard algebraic methods only allows their rapid implementation and test under real-world conditions. We demonstrate this concept in three different application areas. First, we present a semi-algebraic framework to compute the Canonical Polyadic (CP) decompositions of multidimensional signals, a very fundamental tool in multilinear algebra with applications ranging from chemistry over communications to image compression. Compared to state-of-the art iterative solutions, our framework offers a flexible control of the complexity-accuracy trade-off and is less sensitive to badly conditioned data. The second application area is multidimensional subspace-based high-resolution parameter estimation with applications in RADAR, wave propagation modeling, or biomedical imaging. We demonstrate that multidimensional signals can be represented by tensors, providing a convenient description and allowing to exploit the multidimensional structure in a better way than using matrices only. Based on this idea, we introduce the tensor-based subspace estimate which can be applied to enhance existing matrix-based parameter estimation schemes significantly. We demonstrate the enhancements by choosing the family of ESPRIT-type algorithms as an example and introducing enhanced versions that exploit the multidimensional structure (Tensor-ESPRIT), non-circular source amplitudes (NC ESPRIT), and both jointly (NC Tensor-ESPRIT). To objectively judge the resulting estimation accuracy, we derive a framework for the analytical performance assessment of arbitrary ESPRIT-type algorithms by virtue of an asymptotical first order perturbation expansion. Our results are more general than existing analytical results since we do not need any assumptions about the distribution of the desired signal and the noise and we do not require the number of samples to be large. At the end, we obtain simplified expressions for the mean square estimation error that provide insights into efficiency of the methods under various conditions. The third application area is bidirectional relay-assisted communications. Due to its particularly low complexity and its efficient use of the radio resources we choose two-way relaying with a MIMO amplify and forward relay. We demonstrate that the required channel knowledge can be obtained by a simple algebraic tensor-based channel estimation scheme. We also discuss the design of the relay amplification matrix in such a setting. Existing approaches are either based on complicated numerical optimization procedures or on ad-hoc solutions that to not perform well in terms of the bit error rate or the sum-rate. Therefore, we propose algebraic solutions that are inspired by these performance metrics and therefore perform well while being easy to compute. For the MIMO case, we introduce the algebraic norm maximizing (ANOMAX) scheme, which achieves a very low bit error rate, and its extension Rank-Restored ANOMAX (RR-ANOMAX) that achieves a sum-rate close to an upper bound. Moreover, for the special case of single antenna terminals we derive the semi-algebraic RAGES scheme which finds the sum-rate optimal relay amplification matrix based on generalized eigenvectors. Numerical simulations evaluate the resulting system performance in terms of bit error rate and system sum rate which demonstrates the effectiveness of the proposed algebraic solutions

Singular Value Computation and Subspace Clustering

In this dissertation we discuss two problems. In the first part, we consider the problem of computing a few extreme eigenvalues of a symmetric definite generalized eigenvalue problem or a few extreme singular values of a large and sparse matrix. The standard method of choice of computing a few extreme eigenvalues of a large symmetric matrix is the Lanczos or the implicitly restarted Lanczos method. These methods usually employ a shift-and-invert transformation to accelerate the speed of convergence, which is not practical for truly large problems. With this in mind, Golub and Ye proposes an inverse-free preconditioned Krylov subspace method, which uses preconditioning instead of shift-and-invert to accelerate the convergence. To compute several eigenvalues, Wielandt is used in a straightforward manner. However, the Wielandt deflation alters the structure of the problem and may cause some difficulties in certain applications such as the singular value computations. So we first propose to consider a deflation by restriction method for the inverse-free Krylov subspace method. We generalize the original convergence theory for the inverse-free preconditioned Krylov subspace method to justify this deflation scheme. We next extend the inverse-free Krylov subspace method with deflation by restriction to the singular value problem. We consider preconditioning based on robust incomplete factorization to accelerate the convergence. Numerical examples are provided to demonstrate efficiency and robustness of the new algorithm. In the second part of this thesis, we consider the so-called subspace clustering problem, which aims for extracting a multi-subspace structure from a collection of points lying in a high-dimensional space. Recently, methods based on self expressiveness property (SEP) such as Sparse Subspace Clustering and Low Rank Representations have been shown to enjoy superior performances than other methods. However, methods with SEP may result in representations that are not amenable to clustering through graph partitioning. We propose a method where the points are expressed in terms of an orthonormal basis. The orthonormal basis is optimally chosen in the sense that the representation of all points is sparsest. Numerical results are given to illustrate the effectiveness and efficiency of this method

Constructing networks of quantum channels for state preparation

Entangled possibly mixed states are an essential resource for quantum computation, communication, metrology, and the simulation of many-body systems. It is important to develop and improve preparation protocols for such states. One possible way to prepare states of interest is to design an open system that evolves only towards the desired states. A Markovian evolution of a quantum system can be generally described by a Lindbladian. Tensor networks provide a framework to construct physically relevant entangled states. In particular, matrix product density operators (MPDOs) form an important variational class of states. MPDOs generalize matrix product states to mixed states, can represent thermal states of local one-dimensional Hamiltonians at sufficiently large temperatures, describe systems that satisfy the area law of entanglement, and form the basis of powerful numerical methods. In this work we develop an algorithm that determines for a given linear subspace of MPDOs whether this subspace can be the stable space of some frustration free k-local Lindbladian and, if so, outputs an appropriate Lindbladian. We proceed by using machine learning with networks of quantum channels, also known as quantum neural networks (QNNs), to train denoising post-processing devices for quantum sources. First, we show that QNNs can be trained on imperfect devices even when part of the training data is corrupted. Second, we show that QNNs can be trained to extrapolate quantum states to, e.g., lower temperatures. Third, we show how to denoise quantum states in an unsupervised manner. We develop a novel quantum autoencoder that successfully denoises Greenberger-Horne-Zeilinger, W, Dicke, and cluster states subject to spin-flip, dephasing errors, and random unitary noise. Finally, we develop recurrent QNNs (RQNNs) for denoising that requires memory, such as combating drifts. RQNNs can be thought of as matrix product quantum channels with a quantum algorithm for training and are closely related to MPDOs. The proposed preparation and denoising protocols can be beneficial for various emergent quantum technologies and are within reach of present-day experiments

Software for Exascale Computing - SPPEXA 2016-2019

This open access book summarizes the research done and results obtained in the second funding phase of the Priority Program 1648 "Software for Exascale Computing" (SPPEXA) of the German Research Foundation (DFG) presented at the SPPEXA Symposium in Dresden during October 21-23, 2019. In that respect, it both represents a continuation of Vol. 113 in Springer’s series Lecture Notes in Computational Science and Engineering, the corresponding report of SPPEXA’s first funding phase, and provides an overview of SPPEXA’s contributions towards exascale computing in today's sumpercomputer technology. The individual chapters address one or more of the research directions (1) computational algorithms, (2) system software, (3) application software, (4) data management and exploration, (5) programming, and (6) software tools. The book has an interdisciplinary appeal: scholars from computational sub-fields in computer science, mathematics, physics, or engineering will find it of particular interest