Structuuraanbrengende diagonaalschalingen
Thomas Loots
14 juli 2017
Bachelorproject Begeleiding: Jan BrandtsKorteweg-de Vries Instituut voor Wiskunde
Samenvatting
We geven een bewijs voor de stelling van Sinkhorn-Knopp, die zegt dat iedere niet-negatieve matrix met dubbelstochastisch patroon diagonaal equivalent is aan een unieke dubbelstochas-tische matrix [10]. Voor matrices die bovendien symmetrisch zijn bewijzen we dat links en rechts dezelfde schaling kan worden gebruikt. We beschrijven lineair en kwadratisch conver-gente algoritmes om de schalingsmatrices te benaderen ([7], [11], [15]).
Verder bewijzen we constructief dat iedere positieve vierkante matrix diagonaal equivalent is aan een dergelijke matrix waarvan de rij- en kolomproducten gelijk zijn aan ´e´en. Voor niet-negatieve matrices nemen we producten over niet-nulentries en kennen we geen exacte formule voor de schalingsmatrices. Zodoende geven we een Sinkhorn-Knopp-achtig algoritme om deze te benaderen.
Er blijkt dat iedere (complexe) unitaire matrix naar quasistochastische unitaire vorm kan worden geschaald [9]. De schalingsmatrices kunnen vermoedelijk worden benaderd met een Sinkhorn-Knopp-achtig algoritme [1].
Tot slot besteden we aandacht aan de stelling van Birkhoff-von Neumann, omdat er re-centelijk een variant op deze stelling is bewezen voor quasistochastische unitaire matrices ([2]).
Titel: Structuuraanbrengende diagonaalschalingen
Auteur: Thomas Loots, thomas.loots@student.uva.nl, 10804579 Begeleiding: Jan Brandts
Einddatum: 14 juli 2017
Korteweg-de Vries Instituut voor Wiskunde Universiteit van Amsterdam
Inhoudsopgave
1. Inleiding 4
2. Opzet 6
2.1. Definities en proposities . . . 6
2.2. De afwezigheid van een exacte formule . . . 10
3. Dubbelstochastische schalingen 13 3.1. De stelling van Sinkhorn-Knopp . . . 13
3.2. DS-algoritmes . . . 18
3.2.1. Lineair convergent (SK, DEV, BAL en EQ) . . . 18
3.2.2. Kwadratisch convergent (FU) . . . 19
4. Andere schalingen 23 4.1. Dubbel-product-stochastische matrices . . . 23
4.1.1. Een gunstig verband tussen DS-algoritmes en productbalancing . . . . 24
4.1.2. BvN voor dubbel-product-stochastische matrices . . . 25
4.2. Quasistochastische unitaire matrices . . . 25
5. Conclusie 28 6. Populaire samenvatting 29 Bibliografie 31 A. BvN 32 B. FU 34 C. PB 35
1. Inleiding
In 1964 [15] bewees Richard Sinkhorn de volgende stelling.
Stelling 1.1. Voor iedere positieve vierkante matrix A bestaat er een unieke dubbelstochasti-sche matrix M die te schrijven is als M = RAC, waarbij R en C diagonaalmatrices zijn.
Het links en rechts vermenigvuldigen met diagonaalmatrices R en C noemen we het schalen van matrices en we verwijzen soms naar R en C als (diagonaal)schalingen.
Stelling 1.1 maakt deel uit van een groter geheel, waarvan in 2016 [8] een overzicht werd gegeven door Martin Idel. Dit overzicht vormt een terugblik op meer dan 70 jaar aan het schalen van matrices. Het uitgangspunt is het volgende vraagstuk. Zij A een niet-negatieve vierkante matrix. Zijn er diagonaalmatrices R en C te vinden zodat RAC dubbelstochastisch is? We nemen hetzelfde uitgangspunt en zijn ook ge¨ınteresseerd in hoe (snel) R en C gevonden kunnen worden. Richard Sinkhorn en Paul Knopp bewezen in 1967 [10] precies wanneer deze diagonaalschalingen bestaan. We proberen het originele bewijs moderner te presenteren en bewijzen een vergelijkbaar resultaat voor symmetrische matrices. In Figuur 1.1 [8] wordt een grafisch overzicht gegeven van (verbanden tussen) de vele bewijstechnieken voor de stelling van Sinkhorn-Knopp.
Deze bewijstechnieken brengen een verscheidenheid aan generalisaties en gerelateerde pro-blemen voort. Een voorbeeld is het schalen van (niet-negatieve) matrices naar matrices met voorgeschreven rij- en kolomsommen. Onder zekere voorwaarden kan zelfs worden aange-toond dat er een enkele diagonaalmatrix D bestaat zodat DAD−1 of DAD voorgeschreven rij- en kolomsommen heeft. Ook wordt regelmatig aangenomen dat A symmetrisch is.
Een gerelateerd vraagstuk is of er voor iedere (complexe) unitaire matrix U unitaire diago-naalmatrices R en C bestaan zodat de rij- en kolomsommen van RU C gelijk zijn aan ´e´en. In 2014 ([2], [3]) en 2015 [9] werd deze vraag bestudeerd en positief beantwoord. Naar aanleiding hiervan werd in 2016 ([1], [4]) door een aantal van dezelfde auteurs een variant gegeven op de stelling van Garrett Birkhoff en John von Neumann. Deze stelling zegt dat iedere dubbelsto-chastische matrix geschreven kan worden als convexe combinatie van permutatiematrices. De nieuwe variant geeft een vergelijkbare uitdrukking voor quasistochastische unitaire matrices. We zullen een kort overzicht van deze recente resultaten geven.
Garrett Birkhoff (1911-1996). John von Neumann (1903-1957). Tot slot kunnen niet-negatieve matrices onder zekere voorwaarden worden geschaald naar een unieke matrix van dezelfde afmetingen met voorgeschreven rij- en kolomproducten [14]. Voor positieve matrices zijn er exacte formules beschikbaar voor de diagonaalschalingen. Voor niet-negatieve matrices benaderen we de diagonaalschalingen iteratief.
2. Opzet
Sectie 2.1 vormt een wiskundige basis voor volgende hoofdstukken. In Sectie 2.2 verklaren we waarom R en C uit Stelling 1.1 in het algemeen benaderd moeten worden.
2.1. Definities en proposities
Schrijf {e1, . . . , en} voor de standaardbasis van Rn, e = (1, 1, . . . , 1)| voor de all-ones vector,
0 voor de nulmatrix en I voor de identiteitsmatrix. Laat Rm×n+ (resp. R m×n
≥0 ) de verzameling
van positieve (resp. niet-negatieve) m × n-matrices zijn. Schrijf Rm×n6=0 voor de verzameling
van m × n-matrices waarvan alle entries ongelijk zijn aan nul. Voor matrices A, B ∈ Rm×n schrijven we A > B als voor alle entries geldt dat aij > bij (analoog voor <, ≤ en ≥).
Definieer voor v ∈ Rn de diagonaalmatrix
diag(v) = diag(v1, . . . , vn) = v1 . .. vn en v−1 = (v−11 , . . . , vn−1)| voor v ∈ Rn
6=0. De commutatieve groep van diagonaalmatrices met
louter positieve getallen op de hoofddiagonaal duiden we aan met Dn+. Merk op dat de inverse van een willekeurig element diag(v) ∈ Dn+wordt gegeven door diag(v−1). Veel schrijfwijzen en definities voor Rm×nworden in Sectie 4.2 zonder verdere vermelding ook voor Cm×ngebruikt.
Definitie 2.1. Zij A ∈ Rm×n. Het getal e|iAe (resp. e|Aej) heet de i-de rijsom (resp. j-de
kolomsom) van A. Het i-de rijproduct (resp. j-de kolomproduct ) van A wordt gegeven door Y j:aij6=0 aij (resp. Y i:aij6=0 aij).
Hierbij is de conventie dat het lege product gelijk is aan ´e´en.
Definitie 2.2. Een matrix M ∈ Rn×n≥0 heet rijstochastisch (resp. kolomstochastisch) als
M e = e (resp. e|M = e|), oftewel als alle rijsommen (resp. kolomsommen) gelijk zijn aan ´e´en. Matrices die zowel rijstochastisch als kolomstochastisch zijn heten dubbelstochastisch. Algemener heet een matrix U ∈ Rn×n quasistochastisch als U e = e en e|U = e|. Een matrix P heet dubbel-product-stochastisch als alle rij- en kolomproducten van P gelijk zijn aan ´e´en.
De (hoofd)diagonaal van een matrix kan als volgt worden gegeneraliseerd.
Definitie 2.3. Zij A ∈ Rn×n en σ ∈ Sn. De entries a1,σ(1), . . . , an,σ(n) vormen een diagonaal
van A. Als al deze entries positief (resp. niet-nul) zijn, spreken we van een positieve diagonaal (resp. niet-nuldiagonaal ).
Definitie 2.4. De permanent van een matrix A ∈ Rn×n wordt gegeven door per(A) = X σ∈Sn Y 1≤i≤n ai,σ(i).
Het bewijs van de volgende propositie is direct.
Propositie 2.5. Voor elke matrix A ∈ Rn×n≥0 geldt dat A een positieve diagonaal heeft als en
alleen als per(A) > 0.
Definitie 2.6. Als voor matrices A, B ∈ Rm×ngeldt dat aij = 0 als en alleen als bij = 0, dan
hebben ze hetzelfde patroon. Voor A ∈ Rn×n spreken we van een dubbelstochastisch patroon
als er een dubbelstochastische matrix bestaat met hetzelfde patroon als A.
Definitie 2.6 is inzichtelijk, maar niet constructief. De volgende propositie biedt uitkomst (zie [13] voor een bewijs).
Propositie 2.7. Zij A ∈ Rn×n met A 6= 0. De volgende uitspraken zijn equivalent. 1. A heeft een dubbelstochastisch patroon.
2. Iedere niet-nulentry van A behoort tot een niet-nuldiagonaal. 3. Er bestaan geen permutatiematrices Π1 en Π2 zodat
Π1AΠ2= A1 0 A2 A3 , waarbij A1∈ Rk×k (0 < k < n) en A2 6= 0.
De volgende definitie vormt een verscherping van (irreducibiliteit en) de derde uitspraak in de vorige propositie.
Definitie 2.8. Een matrix A ∈ Rn×n heet partly decomposable als er permutatiematrices Π1
en Π2 bestaan zodat Π1AΠ2 = A1 0 A2 A3 ,
waarbij A1 ∈ Rk×k (0 < k < n). We noemen A fully indecomposable als deze niet partly
decomposable is.
De definities en proposities worden verduidelijkt in het volgende voorbeeld. Voorbeeld 2.9. Wegens Propositie 2.7 heeft de matrix
A = 1 0 0 0 2 8 0 1 4
met per(A) = (1 · 2 · 4) + (1 · 8 · 1) = 16 (Def. 2.4) een dubbelstochastisch patroon (Def. 2.6), aangezien elke niet-nulentry bevat is in (minstens) ´e´en van de niet-nuldiagonalen (Def. 2.3)
1 0 0 0 2 8 0 1 4 en 1 0 0 0 2 8 0 1 4 .
Zodoende bestaat er een dubbelstochastische matrix (Def. 2.2) met hetzelfde patroon als A, bijvoorbeeld M = 1 0 0 0 1/2 1/2 0 1/2 1/2 .
Merk op dat A en M partly decomposable zijn (door Π1 = I = Π2 te nemen in Definitie 2.8).
Definitie 2.10. Matrices A, M ∈ Rn×nheten diagonaal equivalent als er matrices R, C ∈ D+n bestaan zodat M = RAC.
Met gebruik van de groepseigenschappen van Dn+ ziet men eenvoudig dat diagonale equi-valentie een equiequi-valentierelatie is. De bewijzen van de volgende twee resultaten zijn direct. Propositie 2.11. Diagonaal equivalente matrices hebben hetzelfde patroon.
Gevolg 2.12. Voor diagonaal equivalente matrices A en M geldt dat A fully indecomposable is als en alleen als M fully indecomposable is.
Definitie 2.13. De directe som van vierkante matrices A1, . . . , Amis de blokdiagonale matrix
A1⊕ · · · ⊕ Am = A1 . .. Am
Propositie 2.14. Als A ∈ Rn×n een dubbelstochastisch patroon heeft, dan bestaan er per-mutatiematrices Π1 en Π2 zodat Π1AΠ2 een directe som is van fully indecomposable matrices
(met dubbelstochastisch patroon).
Bewijs. We bewijzen met inductie naar n.
• Basisstap: Voor n = 1 is de propositie triviaal (met Π1= 1 = Π2).
• Inductiehypothese: Stel dat de propositie geldt voor 1, . . . , n − 1.
• Inductiestap: Zij A ∈ Rn×n zoals in de propositie. Als A fully indecomposable is, neem
dan Π1= I = Π2. Als A partly decomposable is, bestaan er permutatiematrices Π1 en
Π2 zodat Π1AΠ2= A1 0 A2 A3 ,
waarbij A1 ∈ Rk×k (0 < k < n). Wegens Propositie 2.7 geldt dat A2 = 0, dus A =
A1⊕ A3. De matrices A1 en A3 hebben een dubbelstochastisch patroon, anders kon A
wegens Propositie 2.7 via vermenigvuldiging met permutatiematrices op de vorm
A1 0 0 A3 = A11 0 0 0 A12 A13 0 0 0 0 A31 0 0 0 A32 A33
en dus op de vorm A11 0 0 0 0 A31 0 0 A12 0 A13 0 0 A32 0 A33
met A12 6= 0 en A32 6= 0 worden gebracht (door de tweede en derde rij en kolom
van blokmatrices te verwisselen), maar dan heeft A wegens Propositie 2.7 geen dub-belstochastisch patroon. Wegens de inductiehypothese bestaan er permutatiematrices Π11, Π12, Π31, en Π32 zodat Π11A1Π12 en Π31A3Π32 directe sommen zijn van fully
indecomposable matrices (met dubbelstochastisch patroon). Hieruit volgt direct dat Π1
en Π2 in Π1AΠ2 = Π11 0 0 Π31 A1 0 0 A3 Π12 0 0 Π32 =Π11A1Π12 0 0 Π31A3Π32
permutatiematrices zijn zodat Π1AΠ2 een directe som is van fully indecomposable
ma-trices (met dubbelstochastisch patroon).
Gevolg 2.15. Als M ∈ Rn×n≥0 dubbelstochastisch is, dan bestaan er permutatiematrices Π1 en
Π2 zodat Π1AΠ2 een directe som is van fully indecomposable dubbelstochastische matrices.
Bewijs. Wegens Propositie 2.14 bestaan er permutatiematrices Π1 en Π2 zodat Π1M Π2 een
directe som is van fully indecomposable matrices M1, . . . , Mk. Omdat Π1M Π2
dubbelsto-chastisch is, zijn M1, . . . , Mk ook dubbelstochastisch.
Definitie 2.16. Zij A ∈ Rm×n, E ⊂ {1, . . . , m} en F ⊂ {1, . . . , n}. We schrijven A[E, F ] = (Aij)i∈E,j∈F en A(E, F ) = (Aij)i /∈E,j /∈F voor deze submatrices van A. Verder schrijven we
A(E, F ] = (Aij)i /∈E,j∈F en A[E, F ) = (Aij)i∈E,j /∈F.
Opmerking 2.17. Gevolg 2.15 (resp. Propositie 2.14) toont aan dat iedere dubbelstochasti-sche matrix (resp. matrix met dubbelstochastisch patroon) volledig is opgebouwd uit nullen en fully indecomposable dubbelstochastische submatrices (resp. fully indecomposable subma-trices met dubbelstochastisch patroon) die geen gemeenschappelijke entries hebben.
Voorbeeld 2.18. De matrix M = 0.2 0 0 0.6 0 0.2 0 0.8 0 0 0.2 0 0.4 0 0 0.4 0 0.2 0.4 0 0 0 0 0.6 0 0 1 0 0 0 0 0.2 0 0 0.8 0 = 0.2 0 0 0 0 0 0 0.2 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.1 0 0 0 0 0 0 1 0 0 0 0 0 0 0.1 5 0 0 15 0 5 0 40 0 0 1 0 4 0 0 4 0 2 20 0 0 0 0 30 0 0 1 0 0 0 0 20 0 0 8 0 0.2 0 0 0 0 0 0 0.1 0 0 0 0 0 0 1 0 0 0 0 0 0 0.2 0 0 0 0 0 0 1 0 0 0 0 0 0 0.2 = RAC
is dubbelstochastisch. De matrices A en M zijn diagonaal equivalent (Def. 2.10) en hebben hetzelfde patroon (Prop. 2.11). Gebruikmakend van permutatiematrices Π1 en Π2
corres-ponderend met permutaties (24)(56) en (24)(36) ziet men dat A en M partly decomposable zijn (Gevolg 2.12). Bovendien is Π1AΠ2 de directe som (Def 2.13) van fully indecomposable
matrices met dubbelstochastisch patroon (Prop. 2.14). Oftewel A is opgebouwd uit nullen en de fully indecomposable submatrices (Def. 2.16)
A1= A[{1, 3, 4}, {1, 4, 6}] (rood), A2 = A[{2, 6}, {2, 5}] (blauw), A3= A[{5}, {3}] (groen)
met dubbelstochastisch patroon (Opmerking 2.17). De matrix Π1M Π2 is de directe som van de
bijbehorende fully indecomposable dubbelstochastische matrices M1, M2 en M3 (Gevolg 2.15).
2.2. De afwezigheid van een exacte formule
Het lijkt de bewijstechnieken van Stelling 1.1 en generalisaties ([8], [9], [10], [14], [15]) weleens te ontbreken aan elegantie. Niet-constructieve bewijzen treden ver buiten het gebied van lineaire algebra, terwijl constructieve bewijzen omslachtig lijken. In deze sectie proberen we aannemelijk te maken dat R en C uit Stelling 1.1 in het algemeen niet op een voor de hand liggende wijze gevonden kunnen worden.
Soms kunnen R en C wel eenvoudig gevonden worden. Zo is het triviaal dat voor 1 × 1-matrices A > 0 geldt dat RAC met R = 1 en C = A−1 dubbelstochastisch is. Voor rang-1-matrices kan de all-ones matrix ee| handig worden ingezet om R en C te vinden.
Voorbeeld 2.19. Zij A ∈ Rn×n+ een rang-1-matrix. Dan zijn er vectoren v, w ∈ Rn+ zodat
A = vw|= diag(v)ee|diag(w).
Er volgt dat RAC = n−1diag(v−1) A diag(w−1) = n−1ee| dubbelstochastisch is. Voor 2 × 2-matrices helpt het om de entries te schrijven als kwadraten.
Voorbeeld 2.20. Zij A ∈ R2×2+ . Omdat alle entries van A positief zijn, kunnen ze worden
geschreven als kwadraten:
A =a
2 b2
c2 d2
. Elke rij- en kolomsom van
ˆ RAC =cd 0 0 ab a2 b2 c2 d2 bd 0 0 ac =cda 2bd cdb2ac abc2bd abd2ac = abcdad bc bc ad .
is gelijk aan µ = abcd(ad + bc). Er volgt direct dat RAC = µ−1RAC dubbelstochastisch is.ˆ De vraag is of men R en C kan vinden als A ∈ Rn×n+ willekeurig is. Een idee is recursieve
Voorbeeld 2.21. Zij
A =A1 A2 A3 A4
met Ai ∈ R2×2+ . In Voorbeeld 2.20 zijn exacte formules voor Ri en Ci gegeven zodat Mi =
RiAiCi dubbelstochastisch is. De matrix
M = 1 2
M1 M2
M3 M4
is dubbelstochastisch. Er is echter geen voor de hand liggende manier om bijbehorende diago-naalmatrices R en C met M = RAC te construeren uit de matrices Ri en Ci. Dit probleem
lijkt typerend te zijn voor recursieve en lokale strategie¨en.
In Stelling 4.2 zullen voor elke A ∈ Rn×n+ matrices R0, C0 ∈ Dn+ worden gegeven zodat
P = R0AC0 dubbel-product-stochastisch is. Een voor de hand liggend idee is om de logaritme te gebruiken om de dubbel-product-stochastische matrix P te tranformeren tot een dubbels-tochastische matrix M . Als men deze transformatie kan schrijven als vermenigvuldiging met matrices uit Dn
+, dan geeft dit aanleiding tot exacte formules voor matrices R, C ∈ D+n met
M = RAC. Voorbeeld 2.22. Zij A = 5 10 3 1 5 3 3 30 1 ∈ R3×3+ .
Met Stelling 4.2 vinden we (op afrondfouten na) de dubbel-product-stochastische matrix
P = R0AC0 ≈ 0.73 0 0 0 1.58 0 0 0 0.87 5 10 3 1 5 3 3 30 1 0.41 0 0 0 0.09 0 0 0 0.48 ≈ 1.48 0.64 1.06 0.64 0.69 2.27 1.06 2.27 0.42 . De matrix M ≈ 0.73 −0.11 0.39 −0.11 −0.04 1.15 0.39 1.15 −0.54
gegeven door mij = n−1+ log(pij) is (op afrondfouten na) quasistochastisch, maar niet
dub-belstochastisch en kan niet worden geschreven als het product van P en matrices uit D3 +.
De reden dat we moeite hebben om R en C te vinden wordt geschetst in het volgende voorbeeld.
Voorbeeld 2.23. Zij A ∈ R3×3+ . Om te garanderen dat RAC = r1 0 0 0 r2 0 0 0 r3 a11 a12 a13 a21 a22 a23 a31 a32 a33 c1 0 0 0 c2 0 0 0 c3 = r1a11c1 r1a12c2 r1a13c3 r2a21c1 r2a22c2 r2a23c3 r3a31c1 r3a32c2 r3a33c3
dubbelstochastisch is, moet men het stelsel r1a11c1 + r1a12c2 + r1a13c3 = 1 r2a21c1 + r2a22c2 + r2a23c3 = 1 r3a31c1 + r3a32c2 + r3a33c3 = 1 r1a11c1 + r2a21c1 + r3a31c1 = 1 r1a12c2 + r2a22c2 + r3a32c2 = 1 r1a13c3 + r2a23c3 + r3a33c3 = 1
met rij- en kolomsommen van RAC oplossen naar ri en cj. De vergelijkingen zijn lineair in
de producten ricj, maar niet-lineair (namelijk bilineair) in ri en cj. Zodoende is een exacte
formule voor R en C in het algemeen niet beschikbaar.
3. Dubbelstochastische schalingen
In Sectie 3.1 bewijzen we de stelling van Sinkhorn-Knopp en een vergelijkbaar resultaat voor symmetrische matrices. Vervolgens beschouwen we in Sectie 3.2 algoritmes voor het benaderen van R en C.
3.1. De stelling van Sinkhorn-Knopp
Stelling 1.1 kan worden gegeneraliseerd naar niet-negatieve matrices A [10]. Uniciteit van de dubbelstochastische matrix M wordt bewezen in Stelling 3.1. We bewijzen in Stelling 3.5 (met behulp van [16]) dat de rij in Definitie 3.3 convergeert naar een dubbelstochastische matrix. Hiermee wordt in Stelling 3.6 bewezen dat M bestaat als en alleen als A een dubbelstochastisch patroon heeft.
Merk voor het bewijs van de volgende stelling op dat Π−1= Π|voor iedere
permutatiema-trix Π.
Stelling 3.1 (Uniciteit). Zij A ∈ Rn×n≥0 . Als er een dubbelstochastische matrix M = RAC
met R, C ∈ Dn
+ bestaat, dan is M uniek. De matrices R en C zijn op een constante na uniek
als en alleen als A fully indecomposable is.
Bewijs. Stel dat M = RAC en M0 = R0AC0 met R, C, R0, C0 ∈ Dn
+ dubbelstochastisch zijn en schrijf P = R0R−1 en Q = C−1C0. Zij E = {i | pii= min i pii= p} en F = {j | qjj = maxj qjj = q}. Merk op dat e = M0e = R0AC0e = P RACQe = P M Qe ≥ pM Qe. Voor entries i ∈ E treedt gelijkheid op:
1 = e|ie = e|ipM Qe,
oftewel
p−1 = e|iM Qe ≤ e|iM qe = qe|ie = q, (3.1) waarbij gelijkheid optreedt als en alleen als j /∈ F =⇒ aij = 0. Zo ook
e|= e|M0 = e|P M Q ≤ e|P M q. Voor entries j ∈ F treedt gelijkheid op:
1 = e|ej = e|P T qej,
oftewel
waarbij gelijkheid optreedt als en alleen als i /∈ E =⇒ aij = 0. Er geldt p−1 ≤ q en p ≤ q−1,
oftewel pq = 1. Zodoende treedt er gelijkheid op in Formules 3.1 en 3.2, dus
A[E|F ) = 0 en A(E|F ] = 0. (3.3)
als E, F 6= {1, . . . , n}. Er bestaan permutatiematrices Π1 en Π2 zodat
Π1AΠ2 =
A[E, F ] A[E, F ) A(E, F ] A(E, F )
. (3.4)
Als A fully indecomposable is, dan geldt E = F = {1, . . . , n}. Anders gold A[E|F ) = 0 en A(E|F ] = 0 wegens Formule 3.3 en was A wegens Formule 3.4 partly decomposable. Hieruit volgt direct dat R0R−1= P = pI en C−1C0 = Q = qI, oftewel R0 = pR en C0 = p−1C. Dus R en C zijn op een constante na uniek. Bovendien volgt M = RAC = pRAp−1C = R0AC0= M0, dus M is uniek.
Als A partly decomposable is, bestaan er wegens Propositie 2.14 permutatiematrices Π1 en
Π2 zodat
Π1AΠ2 = A1⊕ · · · ⊕ Am,
waarbij elke Ai fully indecomposable is. Wegens de voorgaande alinea geldt voor elke i dat
een dubbelstochastische matrix Mi= RiAiCi met Ri, Ci∈ Dn+i uniek is als deze bestaat. Er
volgt direct dat als de matrix
ˆ M = M1 . .. Mm = R1 . .. Rm A1 . .. Am C1 . .. Cm = ˆRΠ1AΠ2Cˆ
bestaat, deze de unieke dubbelstochastische matrix is die diagonaal equivalent is aan Π1AΠ2.
Zodoende is de matrix Π|1M Πˆ |2 = RAC met R = Π|1RΠˆ 1, C = Π2CΠˆ |2 ∈ D+n een
dubbelsto-chastische matrix die diagonaal equivalent is aan A en is uniek wegens uniciteit van ˆM . Stel namelijk dat ˜M = ˜RA ˜C dubbelstochastisch is met ˜R, ˜C ∈ Dn+, dan geldt dat de dubbels-tochastische matrix Π1M Π˜ 2 = Π1RΠ˜ |1Π1AΠ2Π|2CΠ˜ 2 diagonaal equivalent is aan Π1AΠ2 en
dus wegens uniciteit van ˆM gelijk is aan ˆM , oftewel ˜M = Π|1M Πˆ |2. De matrices R en C zijn niet op een constante na uniek (vervang bijvoorbeeld R1 en C1 door 2R1 en 12C1).
We leiden Definitie 3.3 in voor 2 × 2-matrices. Voorbeeld 3.2. Zij A =a b
c d
∈ R2×2≥0 met per(A) > 0. De matrix
A1= AC0 = a b c d (a + c)−1 0 0 (b + d)−1 =(a + c) −1a (b + d)−1b (a + c)−1c (b + d)−1d =a 0 b0 c0 d0 is kolomstochastisch en A2= R1A1 = (a0+ b0)−1 0 0 (c0+ d0)−1 a0 b0 c0 d0 =(a 0+ b0)−1a0 (a0+ b0)−1b0 (c0+ d0)−1c0 (c0+ d0)−1d0
is rijstochastisch. Zo krijgen we door alternerend rechts en links te vermenigvuldigen met matrices uit D+2 een rij (Ak)k∈N waarvan de elementen om en om kolom- en rijstochastisch
zijn. Merk op dat deze constructie niet werkt als per(A) = 0, omdat er dan door nul wordt gedeeld.
Een generalisatie hiervan wordt gegeven in de volgende definitie.
Definitie 3.3. Zij A ∈ Rn×n≥0 met per(A) > 0. Schrijf R0 = I en C0 = diag(e|A)−1. Definieer
recursief Rk = diag(Rk−1ACk−1e)−1Rk−1 en Ck = diag(e|RkACk−1)−1Ck−1 voor k ∈ N. De
rij (A, R0AC0, R1AC0, R1AC1, . . .) = (A0, A1, A2, A3, . . . ) heet de SK-rij van A (naar Sinkhorn
en Knopp). De matrices A2k−1(resp. A2k) zijn kolomstochastisch (resp. rijstochastisch) voor
k ∈ N. Net als in Voorbeeld 3.2 is de SK-rij van A niet goedgedefinieerd als per(A) = 0. Lemma 3.4. Zij A ∈ Rn×n≥0 een matrix met een dubbelstochastisch patroon en aij > 0. Als
(xi,k) en (yj,k) positieve rijen zijn zodat (xi,kyj,k) convergeert naar een positieve limiet, dan
bestaan er convergente positieve rijen (x0i,k) en (yj,k0 ) met positieve limieten zodat x0i,kyj,k0 = xi,kyj,k voor alle k ∈ N.
We bewijzen dit lemma voor positieve matrices. Het bewijs voor niet-negatieve matrices wordt gegeven in [10].
Bewijs. Definieer xi,k0 = x−11,kxi,k en yj,k0 = x1,kyj,k. Dan geldt x0i,ky 0
j,k = xi,kyj,k. Er geldt
x01,k = 1 → 1. Hieruit volgt dat yj,k0 = x1,k0 y0j,k= x1,kyj,k een positieve limiet heeft. Als gevolg
geldt ook voor i < 1 ≤ n dat x0i,k = x0i,ky0j,k(yj,k0 )−1 = xi,kyj,k(yj,k0 )−1 een positieve limiet
heeft.
Stelling 3.5 (Convergentie). Zij A ∈ Rn×n≥0 met per(A) > 0. De SK-rij (A0, A1, A2, . . . ) van
A convergeert naar een dubbelstochastische matrix M .
Bewijs. Merk op dat voor elke rijstochastische n × n-matrix B geldt dat per(B) = X σ∈Sn Y i bi,σ(i)≤X i1 · · ·X in b1,i1· · · bn,in = Y i (bi1+ · · · + bin) = Y i (Be)i= Y i 1 = 1.
Een analoog argument kan worden toegepast op kolomstochastische matrices. Zodoende is de rij (per(Ak))k∈N begrensd door 1, aangezien iedere Ak(k ∈ N) rij- of kolomstochastisch is.
Zonder verlies van algemeenheid nemen we aan dat e|Ae = n. Als µ = e|Ae 6= n, beschou-wen we de matrix A0= µ−1A. Er geldt (Ck) = (µCk0), maar de rijen (Ak)k∈N en (A0k)k∈N zijn
identiek en convergeren dus naar dezelfde limiet. Wegens de AM-GM ongelijkheid geldt voor alle k ∈ N dat per(A2k−1) = per(A2k−2) Y i (e|A2k−2)−1i ≥ per(A2k−2)(n−1e|A2k−2e)−n= per(A2k−2),
met gelijkheid als en alleen als e|A2k−2 = e|, oftewel als A2k−2 dubbelstochastisch is. Op
analoge wijze geldt dat per(A2k) ≥ per(A2k−1), met gelijkheid als en alleen als A2k−1
dub-belstochastisch is. Zodoende is (per(Ak)) een stijgende rij die convergeert naar een zekere
x ∈ (0, 1]. Er volgt direct dat (Ak) dubbelstochastische limietpunten heeft.
Er rest aangetoond te worden dat er slechts ´e´en limietpunt is. Stel dat M en M0 limietpun-ten van de SK-rij zijn en Alk → M en Al0k → M
0. Per constructie zijn er diagonaalmatrices
R0k, Ck0 ∈ Dn +zodat Alk = R 0 kAl0 kC 0
k, dus Alk en Alk0 hebben hetzelfde patroon voor alle k ∈ N.
We tonen aan dat dit ook geldt voor de limietpunten. Uit per(Ak) → x volgt dat per(RkCk)
en per(Rk+1Ck) convergeren naar y = xper(A)−1. Er geldt voor elke permutatie σ dat
Y
i
mi,σ(i)=Y
i
ai,σ(i) lim
k→∞(Rzk)ii(Clk)σ(i),σ(i) = limk→∞per(RzkClk)
Y
i
ai,σ(i)= yY
i
ai,σ(i),
waarbij zk ∈ {lk, lk + 1} voor alle k ∈ N. Hierbij wordt gebruikt dat limieten en eindige
producten verwisselbaar zijn. Een analoog argument geeft Y i m0i,σ(i)= yY i ai,σ(i)= Y i mi,σ(i). (3.5)
voor elke permutatie σ. Stel dat mij > 0. Omdat M dubbelstochastisch is, maakt mij
wegens Propositie 2.7 deel uit van een positieve diagonaal. Het product van de elementen in deze positieve diagonaal is positief, dus wegens Formule 3.5 geldt ook dat m0ij > 0. Net zo impliceert m0ij > 0 dat mij > 0, dus M en M0 hebben hetzelfde patroon.
Zij mij > 0. Dan (R0k)ii(Ck0)jj → mijm0−1ij > 0. Wegens Lemma 3.4 bestaan convergente
positieve rijen (ri,k) → ri > 0 en (cj,k) → cj > 0 met ri,kcj,k = (R0k)ii(Ck0)jj voor alle k ∈ N.
Dan geldt M = diag(r1, . . . , rn) M0 diag(c1, . . . , cn), dus M en M0 zijn diagonaal equivalent.
Wegens Stelling 3.1 geldt M = M0.
Dankzij het voorwerk resteert er een bondig existentiebewijs.
Stelling 3.6 (Existentie). Zij A ∈ Rn×n≥0 . Er bestaat een dubbelstochastische matrix M die
diagonaal equivalent is aan A als en alleen als A een dubbelstochastisch patroon heeft. Bewijs. =⇒ Wegens Propositie 2.11 heeft A net als M een dubbelstochastisch patroon. ⇐= Omdat A een dubbelstochastisch patroon heeft, heeft A wegens Propositie 2.7 een positieve diagonaal, dus wegens Propositie 2.5 geldt dat per(A) > 0. Wegens Stelling 3.5 convergeert de SK-rij van A naar een dubbelstochastische limiet M . Uit formule 3.5 volgt dat A en M hetzelfde patroon hebben. Zij aij > 0. Omdat RkACk → M , geldt
dat (Rk)ii(Ck)jj → mija−1ij > 0. Wegens Lemma 3.4 bestaan convergente positieve rijen
De volgende stelling vormt een combinatie van Stelling 3.1 en Stelling 3.6 voor symmetrische matrices. Voor positieve matrices wordt het bewijs van de volgende stelling gegeven in [8]. We geven een bewijs voor niet-negatieve matrices.
Stelling 3.7. Als A ∈ Rn×n≥0 een dubbelstochastisch patroon heeft en symmetrisch is, dan
bestaat er een matrix X ∈ D+n zodat XAX de unieke dubbelstochastisch matrix is die diagonaal equivalent is aan A. De matrix X is uniek als en alleen als A fully indecomposable of een directe som van fully indecomposable matrices is.
Bewijs. Zij M = RAC met R, C ∈ D+n de unieke dubbelstochastische matrix die diagonaal equivalent is aan A. Er geldt dat
M|= (RAC)|= C|A|R|= CAR = CR−1M C−1R.
Omdat M|een dubbelstochastische matrix is die diagonaal equivalent is aan A, geldt wegens uniciteit dat M = M|. Als M fully indecomposable is, geldt wegens M = IM I en Stelling 3.1 dat CR−1 = λI, oftewel C = λR voor een zekere constante λ > 0. Dan moet X voldoen aan X = αR en X = α−1λR voor een zekere constante α > 0. Oftewel α = α−1λ, dus α =√λ. We concluderen dat X =√λR de unieke matrix is die voldoet.
Als M partly decomposable is, dan bestaat M = RAC wegens Opmerking 2.17 uit (nullen en) fully indecomposable dubbelstochastische submatrices die geen gemeenschappelijke entries hebben. Laat M [E, F ] een dergelijke submatrix zijn en merk op dat M [E, F ]|= M [F, E] ook een fully indecomposable dubbelstochastische submatrix is. De matrices M [E, F ] en M [F, E] zijn daarom dezelfde submatrix of hebben geen gemeenschappelijke entries. Zodoende maken het volgende gevalsonderscheid.
• E = F . Wegens het voorgaande bestaat er een unieke matrix Y ∈ D+#Ezodat M [E, F ] = Y A[E, F ]Y , dus er moet gelden dat X[E, E] = X[F, F ] = Y .
• E ∩ F = ∅. Er geldt
M [F, E] = M [E, F ]|= (R[E, E]A[E, F ]C[F, F ])|= C[F, F ]A[F, E]R[E, E], oftewel
X[E, E] = λR[E, E] en X[F, F ] = λ−1C[F, F ] voldoet voor alle λ > 0, dus X is niet uniek.
We illustreren de stellingen in een voorbeeld. Voorbeeld 3.8. Zij A = 0 0 0 5 0 12 184 0 0 184 184 0 5 0 0 0 ∈ R 4×4 ≥0 .
Deze matrix heeft een dubbelstochastisch patroon, dus er bestaan R, C ∈ D4+zodat M = RAC de unieke (Stelling 3.1) dubbelstochastische matrix is behorend bij A (Stelling 3.6), bijvoorbeeld
M = 0 0 0 1 0 14 34 0 0 34 14 0 1 0 0 0 = 1 0 0 0 0 1 0 0 0 0 13 0 0 0 0 15 0 0 0 5 0 12 184 0 0 184 184 0 5 0 0 0 1 0 0 0 0 12 0 0 0 0 16 0 0 0 0 15 = RAC.
Omdat A symmetrisch is, geldt Stelling 3.7 ook. Als we de constructie in het bewijs van deze stelling volgen, vinden we dat
X = λ 0 0 0 0 12√2 0 0 0 0 16√2 0 0 0 0 5λ1 ∈ D4 +
met λ > 0 de oplossingen zijn van M = XAX. De matrix X is niet uniek omdat A partly decomposable is (neem rijpermutatie (14)) en geen directe som is van fully indecomposable matrices.
Nu we hebben bewezen dat iedere matrix A ∈ Rn×n≥0 met dubbelstochastisch patroon op
dubbelstochastische vorm kan worden gebracht, is de volgende stelling des te interessanter (zie [5] voor een bewijs). Recentelijk is een variant op deze stelling voor quasistochastische unitaire matrices bewezen (zie Stelling 4.14). Convexe combinaties van permutatiematrices zijn duidelijk dubbelstochastisch. De stelling van Birkhoff-von Neumann (BvN) toont aan dat het omgekeerde ook geldt.
Stelling 3.9 (Birkhoff-von Neumann). Iedere dubbelstochastische matrix M kan worden ge-schreven als convexe combinatie van eindig veel permutatiematrices:
M = λ1Π1+ · · · + λmΠm,
waarbij Π1, . . . , Πm permutatiematrices en λ1, . . . , λm∈ [0, 1] scalairen zijn met Pmi=1λi = 1.
Als M een n×n matrix is, dan geldt uiteraard dat m ≤ n!. De bovengrens voor m kan echter aanzienlijk worden aangescherpt tot n2− 2n + 2. In Hoofdstuk A geven we een (afleiding van een) algoritme dat de scalairen en permutatiematrices vindt voor integer matrices waarvan de rij- en kolomsommen gelijk zijn.
3.2. DS-algoritmes
Definitie 3.3 beschrijft een algoritme SK voor het benaderen van R en C. Drie andere algo-ritmes (DEV, BAL, en EQ) met hetzelfde doel worden uitvoerig behandeld in [11]. In Subsectie 3.2.1 wordt lineaire convergentie van deze vier algoritmes experimenteel aangetoond. Kwadra-tische convergentie van algoritme FU (gebaseerd op een algoritme uit [7]) wordt in Subsectie 3.2.2 experimenteel aangetoond. In deze (en latere) sectie(s) worden de testmatrices A, B, C, D, H1, . . . , H5 uit [11] zonder verdere vermelding gebruikt. Voor het gemak noemen we
de algoritmes voor het benaderen van R en C DS-algoritmes (naar dubbelstochastisch).
3.2.1. Lineair convergent (SK, DEV, BAL en EQ)
Elk algoritme schaalt per iteratie een niet-negatieve matrix A (met dubbelstochastisch pa-troon) naar een dergelijke matrix A0. We geven de strategie van elk algoritme.
In SK worden iteratief de rij- en kolomsommen zoals in Definitie 3.3 alternerend op ´e´en geschaald.
gelijk is aan n−1e|A0e. Als bijvoorbeeld van alle rijen en kolommen de som e|i
0Ae van rij i0
het meest afwijkt van n−1e|Ae, dan vermenigvuldigt DEV de elementen in rij i0 van A met
(e|i
0Ae)
−1(n − 1)−1(e − e i0)
|Ae
Hieruit volgt dat e|i0A0e = n−1e|A0e. Verder schaalt BAL het rij-kolompaar
(i0, j0) = argmax (i,j)
|e|iAe − e|Aej|
waarvan de sommen maximaal verschillen, zodat voor hun sommen in A0 geldt dat e|iA0e = e|A0ej
en a0i0j0 = ai0j0.
Tot slot is EQ een versie van DEV waarin een stap van BAL wordt toegepast zodra DEV de laatst geschaalde kolom (resp. rij) opnieuw zou gaan schalen zonder eerst een andere kolom (resp. rij) geschaald te hebben.
Convergentie naar een diagonaal equivalente dubbelstochastische limiet treedt in DEV, BAL en EQ (net als voor SK) op als A een dubbelstochastisch patroon heeft, zo wordt bewezen in [11]. Om het convergentiegedrag te bestuderen, moet men eerst kunnen meten ’hoe dubbels-tochastisch’ een matrix is. Zij λ(X) = n−1e|Xe de gemiddelde rij/kolomsom van een matrix X en laat λ(X)−1X de bijbehorende genormaliseerde matrix zijn. In [11] wordt
µ(X) = max
i max{|e |
iXe − 1|, |e|Xej− 1|}
als maat genomen. Een algoritme termineert wanneer µ(X) < tol · λ(X) voor een gekozen waarde tol. In de genormaliseerde matrix wijken de rij- en kolomsommen dan maximaal tol af van 1. Hoe kleiner µ(X) is, des te ‘dubbelstochastischer’ is λ(X)−1X.
Voor elk algoritme kan per iteratie de convergentiemaat op logaritmische schaal worden geplot. Overigens gebruikt niet elk algoritme hetzelfde aantal floating point operations per iteratie. Zodoende geven dergelijke plots geen informatie over verhoudingen in tijdsduur. Figuur 3.1 (resp. Figuur 3.2) toont dat SK voor sommige matrices relatief weinig (resp. veel) iteraties vereist. Ook is duidelijk te zien dat voor alle algoritmes sprake is van lineaire convergentie.
3.2.2. Kwadratisch convergent (FU)
In deze subsectie leiden we af dat men het zoeken van R en C altijd kan reduceren tot het zoeken van X voor een symmetrische matrix A (zie Stelling 3.7). Vervolgens passen we de methode van Newton-Raphson toe als poging tot het verkrijgen van kwadratische convergentie naar R en C. We vergelijken deze methode met een algoritme FU dat wordt beschreven in [7]. Zij A ∈ Rn×n≥0 een matrix met een dubbelstochastisch patroon. We zoeken R, C ∈ D+n
zodat RAC de unieke dubbelstochastische matrix is die diagonaal equivalent is aan A (zie Stelling 3.1 en Stelling 3.6). We stellen zonder verlies van algemeenheid dat A symmetrisch is. Anders bekijken we de symmetrische matrix
0 A A| 0
Figuur 3.1.: Convergentieplot van vier DS-algoritmes op testmatrix H4.
Figuur 3.2.: Convergentieplot van zes DS-algoritmes op testmatrix D.
Men gaat eenvoudig na dat deze matrix een dubbelstochastisch patroon heeft. Als R 0 0 C 0 A A| 0 R 0 0 C = 0 RAC CA|R 0 dubbelstochastisch is met R, C ∈ Dn
+, dan is RAC dubbelstochastisch. Zodoende nemen we
aan dat A symmetrisch is en zoeken we wegens Stelling 3.7 een matrix X ∈ Dn+ zodat XAX dubbelstochastisch is. Schrijf
Xs= es1 . .. esn
en As= XsAXs voor s ∈ Rn. In feite zoeken we X ∈ D+n zodat XAX rijstochastisch is, want
XAX is dan wegens symmetrie ook kolomstochastisch. Zodoende geeft een nulpunt t van de functie
f : Rn→ Rn: s 7→ Ase − e.
aanleiding tot een oplossing X = Xt. Stelling 3.7 garandeert dat een dergelijk nulpunt bestaat.
We passen de methode van Newton-Raphson toe om t te benaderen.
Definitie 3.10 (Newton-Raphson, [6]). Beschouw g : B ⊂ Rn→ Rnmet g(t) = 0 voor zekere
t ∈ B en veronderstel dat Dg(s) inverteerbaar is voor alle s ∈ B. De dekpuntiteratie s(k+1) = s(k)− [Dg(s(k))]−1g(s(k))
met gegeven startwaarde s(1) ∈ B is de methode van Newton-Raphson. Stelling 3.11. Als A ∈ Rn×n en s ∈ Rn, dan geldt Df (s) = diag(A
se) + As.
Bewijs. We bepalen Df (s) via de formule
We schrijven ≈ voor stappen waarin we termen O(khk2) = o(khk) weglaten. Merk overigens op dat we ondanks het weglaten van termen de totale afgeleide exact berekenen. Schrijf H = diag(h). Wegens ex= 1 + x + O(x2) voor x ∈ R geldt
Xh = eh1 . .. ehn = 1 + h1+ O(h21) . .. 1 + hn+ O(h2n) ≈ I + H.
Hieruit volgt de lineaire benadering
f (s + h) − f (s) = As+he − e − (Ase − e) = (XhAsXh− As)e ≈ ((I + H)As(I + H) − As)e
= (HAs+ AsH + HAsH)e ≈ (HAs+ AsH)e = (diag(Ase) + As)h,
oftewel Df (s) = diag(Ase) + As. Hierbij gebruiken we in de laatste gelijkheid dat
HAse = h1 . .. hn e|1Ase .. . e|nAse = e|1Aseh1 .. . e|nAsehn = e|1Ase . .. e|nAse h1 .. . hn = diag(Ase)h.
Zodoende kunnen we de methode van Newton-Raphson toepassen (zolang er geen proble-men optreden met inverteerbaarheid van de totale afgeleide). Als beginwaarde neproble-men we s(1) = 0 en we gebruiken het eerder beschreven stopcriterium µ(A
s(k)) < tol · λ(As(k)) voor
een gekozen waarde tol.
In [7] wordt een variant op deze methode beschreven. In iedere stap wordt de update u(k)= −[Df (s(k))]−1f (s(k)) berekend. Vervolgens wordt u(k) iteratief vervangen door u(k)/2
zolang max i |e | if (s (k)+ u(k)/2)| < max i |e | if (s (k)+ u(k))|. (3.6)
Pas daarna berekent men s(k+1)= s(k)+ u(k) als correctie op de vorige benadering s(k).
We geven een heuristische verklaring voor deze aanpak. Het is typerend voor de methode van Newton-Raphson dat (s(k)) lineair convergeert naar een nulpunt t, totdat (s(k)) in een nabije omgeving Y van t komt te liggen en er kwadratische convergentie optreedt. Zodoende is het de kunst om (snel) in het gebied Y te geraken. Als de afstand tussen opvolgende benaderingen s(k) en s(k+1) groot is, dan loopt men het risico ’over het gebied Y te stappen’. Zodoende worden de updates u(k) iteratief gehalveerd om het verschil tussen opvolgende
benaderingen klein te houden. Formule 3.6 garandeert dat f (s(k+1)) hierdoor niet verder van nul komt te liggen ten opzichte van k·k∞. In onze experimenten treedt kwadratische
convergentie overigens snel op en veroorzaken de halveringsiteraties geen merkbaar verschil. Deze variant op de methode van Newton-Raphson noemen we FU. In hoofstuk B geven we een Matlab-implementatie van dit algoritme. Figuur 3.3 toont de kwadratische convergentie van FU. Voor testmatrix D is de convergentiemaat na twaalf iteraties kleiner dan 10−15, terwijl de lineair convergente algoritmes duizend(en) iteraties vereisen om een convergentiemaat klei-ner dan 10−5te behalen. Voor alle testmatrices vereist FU minder dan twintig iteraties. Zoals eerder is opgemerkt zegt dit overigens weinig over de tijdsduur van het algoritme. Daarom
Figuur 3.3.: Convergentieplot van DS-algoritme FU op testmatrix D. FU DEV BAL EQ SK A 0.556 0.089 0.725 0.093 0.042 B 0.476 7.461 1.510 1.164 2.302 C 0.602 62.31 7.704 0.988 28.79 D 0.771 135.9 1.069 0.626 36.26 H1 0.757 12.90 12.05 13.30 0.787 H2 0.861 13.90 14.52 11.83 1.116 H3 0.963 14.45 14.94 12.92 0.920 H4 1.043 14.86 15.23 15.03 0.909 H5 1.003 226.1 272.9 15.72 12.78
Tabel 3.1.: Tijdsduur (sec.) van vijf DS-algoritmes voor verschillende testmatrices.
voeren we het volgende experiment uit. Elk algoritme wordt honderd keer toegepast op iedere testmatrix. De nauwkeurigheid is tol = 10−5 en de tijdsduur wordt gemeten in seconden. De resultaten worden gepresenteerd in Tabel 3.1.
Algoritme FU is altijd relatief snel. De testmatrices A, B, C en D (resp. Hi) zijn overigens
geconstrueerd zodat SK (resp. DEV, BAL, EQ) ineffici¨ent zijn. Het is niet uitgesloten dat er matrices bestaan waarvoor FU relatief ineffici¨ent is.
4. Andere schalingen
Tot dusver hebben we alleen geschaald naar dubbelstochastische matrices. In Sectie 4.1 schalen we naar dubbel-product-stochastische matrices en in Sectie 4.2 naar quasistochastische unitaire matrices.
4.1. Dubbel-product-stochastische matrices
We weten nu dat diagonaalschalingen kunnen worden gebruikt om rij- en kolomsommen gelijk te maken aan ´e´en. In deze sectie wordt hetzelfde idee voor rij- en kolomproducten bekeken. Dat wil zeggen: bij iedere matrix A ∈ Rn×n≥0 worden diagonaalmatrices R en C gezocht zodat
P = RAC dubbel-product-stochastisch is. Dit proces noemen we productbalancing. Voorbeeld 4.1. Gegeven is de matrix
A = 1 b3 1 1 1 1 1 1 1 ∈ R3×3+ met b > 0. Zij R = 1 0 0 0 b 0 0 0 b en C = b 0 0 0 1 0 0 0 b ,
dan geldt dat
b−53RAC = b− 5 3 b b3 b b2 b b2 b2 b b2 dubbel-product-stochastisch is.
Door op deze manier entries te ‘verspreiden’ kan iedere matrix A ∈ Rn×n+ op
dubbel-product-stochastische vorm worden gebracht. Voor 2 × 2-matrices is dit precies de constructie uit Voorbeeld 2.20 die leidt tot de matrix ˆRAC. De matrix (abcd)−32RAC is dubbel-product-ˆ
stochastisch. Een andere manier om voor A ∈ Rn×n+ diagonaalmatrices R en C te vinden
wordt beschreven in het bewijs van de volgende stelling.
Stelling 4.2. Voor iedere matrix A ∈ Rn×n+ bestaat er een dubbel-product-stochastische matrix
P die diagonaal equivalent is aan A. Bewijs. Zij C = diag((Q
iai1) −1/n
, . . . , (Q
iain) −1/n
). Elk kolomproduct van AC is gelijk aan ´e´en: Y i cjjaij = Y i (Y i0 ai0j)−1/naij = ( Y i0 ai0j)−1 Y i aij = 1.
Zij R = diag((Q
jcjja1j)
−1/n, . . . , (Q
jcjjanj)
−1/n). Elk rijproduct van P = RAC is gelijk
aan ´e´en: Y j riicjjaij = Y j (Y j0 cj0j0aij0)−1/ncjjaij = ( Y j0 cj0j0aij0)−1 Y j cjjaij = 1.
De kolomproducten van P zijn ook gelijk aan ´e´en: Y i riicjjaij = Y i rii Y i cjjaij = Y i (Y j cjjaij)−1/n· 1 = Y j (Y i cjjaij)−1/n = Y j 1−1/n = 1.
Voor willekeurige matrices A ∈ Rn×n≥0 werkt deze constructie niet langer: in het algemeen
blijven de kolomproducten in de overgang van AC naar RAC niet gelijk aan ´e´en. Ook het verspreiden van entries werkt niet meer. In Hoofdstuk C geven we daarom een Matlab-implementatie van een algoritme PB dat R en C benadert voor matrices A ∈ Rn×n≥0 . In iedere
iteratie wordt een variant op Stelling 4.2 gebruikt. Het algoritme termineert zodra ieder rij-en kolomproduct maximaal tol afwijkt van 1. Naar eigrij-en vermoedrij-en treedt in PB voor alle A ∈ Rn×n≥0 convergentie op naar een dubbel-product-stochastische matrix RAC. Aangezien PB
net als SK alternerend rijen en kolommen schaalt, verwachten we lineaire convergentie. Deze verwachting wordt bevestigd door Figuur 4.1.
Figuur 4.1.: Convergentieplot van PB op testmatrix B.
Figuur 4.2.: Convergentieplot van SK op testmatrix D met herhaalde-lijke toepassing van PB.
4.1.1. Een gunstig verband tussen DS-algoritmes en productbalancing
Voor DS-algoritmes kunnen matrices worden ‘voorbehandeld’ door eerst productbalancing toe te passen. Deze voorbehandeling duiden we aan met het voorvoegsel ‘prod’. Figuur 3.2 toont dat convergentie hierdoor sterk verbeterd kan worden. Een mogelijke verklaring is dat door productbalancing entries meer in de buurt van het interval (0, 1) komen te liggen (zie Voorbeeld 2.22). Figuur 4.2 illustreert in ieder geval dat productbalancing ongeschikt is voor herhaaldelijk gebruik: in SK wordt de convergentie naar een dubbelstochastische matrix iedere
Opmerking 4.3. In DS-algoritmes levert men nauwelijks effici¨entie in met de voorbehande-ling als A ∈ Rn×n+ , omdat productbalancing zoals in Stelling 4.2 nauwelijks werk vereist. Voor
niet-positieve matrices A ∈ Rn×n≥0 moeten we PB toepassen en is effici¨entie niet gegarandeerd.
4.1.2. BvN voor dubbel-product-stochastische matrices
Iedere combinatie
λ1Π1+ · · · + λmΠm (4.1)
van permutatiematrices Π1, . . . , Πm met λ1· · · λm = 1 is dubbel-product-stochastisch.
Zo-doende kan men zich afvragen of er een variant op de stelling van Birkhoff-von Neumann bestaat voor dubbel-product-stochastische matrices: is iedere dubbel-product-stochastische matrix P te schrijven zoals in Formule 4.1?
Voorbeeld 4.4. De nulmatrix of bijvoorbeeld de matrix 1 1
0 1
zijn dubbel-product-stochastisch maar kunnen duidelijk niet worden geschreven zoals in For-mule 4.1. Het baat niet om aan te nemen dat P een dubbelstochastisch patroon heeft. Neem bijvoorbeeld de dubbel-product-stochastische matrix
P = 1 1 1 1/2 1/2 4 2 2 1/4
Iedere entry ligt op precies twee positieve diagonalen. Als we P schrijven als combinatie van de zes 3 × 3-permutatiematrices, forceert a11 = 1 (resp. a22 = 12 en a23 = 4) de eerste twee
matrices (resp. de derde en vierde matrix) in
P = α 1 0 0 0 1 0 0 0 1 + (1 − α) 1 0 0 0 0 1 0 1 0 + ( 1 2− α) 0 0 1 0 1 0 1 0 0 + (4 + α − 1) 0 1 0 0 0 1 1 0 0 + · · · .
Maar nu ligt de waarde 72 6= p31vast op positie (3, 1), dus P kan niet worden geschreven zoals in Formule 4.1.
4.2. Quasistochastische unitaire matrices
Notatie 4.5. SchrijfU (n) = {U ∈ Cn×n | U is unitair},
W U (n) = {W ∈ U (n) | W is quasistochastisch},
ZU (n) = {Z ∈ U (n) | Z is een diagonaalmatrix en Z11= 1}.
.
In artikel [3] uit 2014 wordt het volgende vermoeden uitgesproken.
Vermoeden 4.6. Voor iedere U ∈ U (n) bestaan er diagonaalmatrices R, C ∈ U (n) zodat W = RU C ∈ W U (n).
Eerder in 2014 werd door dezelfde auteurs de volgende stelling bewezen. Stelling 4.7 ([2]). Elke matrix U ∈ U (n) kan worden geschreven als
U = eiφZ1W1· · · Zp−1Wp−1Zp,
met p ≤ n, φ ∈ [0, 2π), Zi ∈ ZU (n) en Wi ∈ W U (n).
Een herformulering van Vermoeden 4.6 is dat p ≤ 2 in Stelling 4.7. Vermoeden 4.8. Elke matrix U ∈ U (n) kan worden geschreven als
U = eiφZ1W Z2, (4.2)
met φ ∈ [0, 2π), Z1, Z2∈ ZU (n) en W ∈ W U (n).
Voorbeeld 4.9. Elk element U ∈ U (2) kan worden geschreven zoals in Vermoeden 4.8: U = eiφ
cos(θ)eiψ sin(θ)eiξ − sin(θ)e−iξ cos(θ)e−iψ
= ei(φ+θ+ψ)1 0 0 ie−i(ψ+ξ)
cos(θ)e−iθ i sin(θ)e−iθ
i sin(θ)e−iθ cos(θ)e−iθ
1 0
0 −iei(−ψ+ξ)
= ei(φ+θ+ψ)Z1W Z2
Zie [3] voor meer details.
Opmerking 4.10. We tonen aan dat Vermoeden 4.6 en Vermoeden 4.8 daadwerkelijk equi-valent zijn. Zij U ∈ U (n). Als er diagonaalmatrices R, C ∈ U (n) bestaan zodat
W = RU C ∈ W U (n),
dan volgt dat U = R−1W C−1 met R−1, C−1 ∈ U (n). Er geldt dat r−111 = eiψ en c−111 = eiξ voor zekere ψ, ξ ∈ [0, 2π). Zij Z1= e−iψR−1, Z2 = e−iξC−1∈ ZU (n), dan geldt dat
U = eiφZ1W Z2,
waarbij eiφ = ei(ψ+ξ) en φ ∈ [0, 2π). Andersom, als U geschreven kan worden zoals in Formule 4.2, dan volgt dat
W = RAC = e−iφZ1−1U Z2−1, waarbij R, C ∈ U (n) diagonaalmatrices zijn.
Vergelijk de volgende definitie met Definitie 3.3. Definitie 4.11. Zij U ∈ U (n) en
φ(y) = (
|y|−1y als y 6= 0, 1 als y = 0,
voor y ∈ C. Schrijf φ(v) = (φ(v1), . . . , φ(vn))| voor v ∈ Cn. Zij R0 = I en C0 =
diag(φ(e|U ))−1. Definieer recursief
Rk= diag(φ(Rk−1U Ck−1e))−1Rk−1 en Ck= diag(φ(e|RkU Ck−1))−1Ck−1.
We noemen
Omdat dit een SK-achtig algoritme is verwachten we lineaire convergentie. Dit wordt ex-perimenteel bevestigd in Figuur 4.3, waarin een convergentieplot wordt getoond van een Matlab-implementatie QS toegepast op een willekeurige matrix U ∈ U (100).
Figuur 4.3.: Convergentieplot van QS op een willekeurige matrix U ∈ U (100). Vermoeden 4.6 is gebaseerd op het volgende vermoeden uit [3].
Vermoeden 4.12. Voor elke U ∈ U (n) convergeert de QS-rij naar een matrix W = RU C ∈ W U (n) met R, C ∈ U (n).
Vermoeden 4.12 wordt gebaseerd op
1. het feit dat de verzamelingen U (n) en {eiφZ1W Z2} beide n2 vrijheidsgraden hebben [3].
2. het bewijs dat W bestaat voor een 2n-dimensionale deelverzameling van U (n) [2]. 3. het succes van duizend numerieke experimenten met n = 3, 4, 5 [3].
4. exacte formules voor Z1 en Z2 in U (2) (zie Voorbeeld 4.9) [3].
Artikel [9] uit 2015 bevestigt niet de convergentie van de QS-rij, maar wel existentie van R en C uit Vermoeden 4.6.
Stelling 4.13. Voor iedere U ∈ U (n) bestaan er diagonaalmatrices R, C ∈ U (n) zodat W = RU C ∈ W U (n).
Naar aanleiding van Stelling 4.13 wordt in [1] een unitaire variant op de stelling van Birkhoff-von Neumann bewezen. Er worden zelfs formules (met keuzevrijheid) voor de co¨effici¨enten zi
gegeven.
Stelling 4.14. Iedere matrix W ∈ W U (n) kan worden geschreven als combinatie van eindig veel permutatiematrices:
W = z1Π1+ · · · + zmΠm,
waarbij Π1, . . . , Πm permutatiematrices en z1, . . . , zm ∈ C scalairen zijn met Pmi=1zi = 1 en
Pm
i=1|zi|2= 1 .
5. Conclusie
We hebben bewezen dat er voor iedere niet-negatieve matrix A met dubbelstochastisch pa-troon diagonaalmatrices R, C ∈ Dn+ bestaan zodat M = RAC dubbelstochastisch is. De matrix M is uniek. De diagonaalmatrices R en C zijn op een constante na uniek als en alleen als A fully indecomposable is. Als A ook symmetrisch is, dan bestaat er zelfs een diagonaalmatrix X ∈ Dn+ zodat M = XAX. De matrix X is uniek als en alleen als A fully indecomposable of een directe som van fully indecomposable matrices is. Omdat er in het algemeen geen directe formule voor R en C beschikbaar is, bekeken we vijf zogenoemde DS-algoritmes voor het iteratief benaderen van R en C. We toonden experimenteel aan dat vier van deze algoritmes (SK, DEV, BAL en EQ) lineair convergent zijn. We bewezen dat het vijfde algoritme FU een variant op de methode van Newton-Raphson is en toonden experimenteel kwadratische convergentie aan. Het zou interessant zijn om de convergentie van FU nader te bestuderen.
Verder bewezen we constructief dat er voor iedere positieve vierkante matrix A diago-naalmatrices R, C ∈ D+n bestaan zodat P = RAC dubbel-product-stochastisch is. Voor niet-negatieve matrices gaven we een ogenschijnlijk lineair convergent SK-achtig algoritme PB om R en C te benaderen en zagen dat er een gunstig verband lijkt te bestaan tussen pro-ductbalancing en DS-algoritmes. In nader onderzoek kan de convergentiefactor van PB worden bestudeerd en kan er ook voor niet-negatieve matrices een exacte formule voor R en C worden gezocht.
Vervolgens bleek dat er voor iedere (complexe) unitaire matrix U unitaire diagonaalma-trices R en C bestaan zodat W = RU C quasistochastisch is. De diagonaalmadiagonaalma-trices R en C kunnen vermoedelijk worden benaderd met een SK-achtig algoritme, dat zoals verwacht line-air convergent lijkt te zijn. Tot slot toonden we een variant op de stelling van Birkhoff-von Neumann, die quasistochastische unitaire matrices uitdrukt als combinatie van permutatie-matrices. In vervolgonderzoek kan het interessant zijn om meer aandacht te besteden aan unitaire matrices.
6. Populaire samenvatting
Gegeven is een matrix A met niet-negatieve entries, bijvoorbeeld
A = 1 0 0 0 2 8 0 1 4 .
De rijsommen (resp. kolomsommen) van A worden gegeven door de vector
rA= 1 10 5 (resp. cA= 1 3 12). De matrix M = 1 0 0 0 1/2 1/2 0 1/2 1/2
is in dit opzicht bijzonder; de rij- en kolomsommen zijn allemaal gelijk aan ´e´en. Oftewel
rM = 1 1 1 en cM = 1 1 1 .
Niet-negatieve matrices waarvan de rij- en kolomsommen allemaal gelijk zijn aan ´e´en noemen we dubbelstochastisch. Verder zien we dat
M = 1 0 0 0 1/2 1/2 0 1/2 1/2 = 1 0 0 0 12 0 0 0 1 1 0 0 0 2 8 0 1 4 1 0 0 0 12 0 0 0 18 = RAC,
oftewel A kan door links- en rechtsvermenigvuldiging met diagonaalmatrices (lees: het schalen van rijen en kolommen) R en C worden getransformeerd naar een dubbelstochastische matrix. Dat dit mogelijk is voor A blijkt te danken te zijn aan de posities van de nullen in deze matrix. We bewijzen precies voor welke nulpatronen dergelijke transformaties mogelijk zijn. Tevens bewijzen we dat de resulterende dubbelstochastische matrix uniek is.
In het vorige voorbeeld vonden we eenvoudig diagonaalmatrices R en C zodat M = RAC dubbelstochastisch is, maar voor
A = 1 1 5 5 3 3 2 1 5 5 1 2 1 3 3 2 2 4 4 5 4 4 2 2 4 met rA= 15 16 10 17 16 en cA= 11 11 13 19 20
liggen deze niet voor de hand. Merk echter op dat we de rijsommen als volgt gelijk kunnen maken aan ´e´en:
R1A = 1 15 0 0 0 0 0 161 0 0 0 0 0 101 0 0 0 0 0 171 0 0 0 0 0 161 1 1 5 5 3 3 2 1 5 5 1 2 1 3 3 2 2 4 4 5 4 4 2 2 4 = 0.07 0.07 0.33 0.33 0.20 0.19 0.13 0.06 0.31 0.31 0.10 0.20 0.10 0.30 0.30 0.12 0.12 0.24 0.24 0.29 0.25 0.25 0.13 0.13 0.25 met rR1A= 1 1 1 1 1 en cR1A= 0.72 0.76 0.86 1.31 1.36 ,
op afrondfouten na. Vervolgens kunnen we op analoge wijze door rechtsvermenigvuldiging met een diagonaalmatrix C1 de kolomsommen gelijk maken aan ´e´en:
rR1AC1 = 0.97 0.97 0.97 0.99 1.1 en cR1AC1 = 1 1 1 1 1 .
Merk echter op dat de rijsommen van deze matrix niet gelijk zijn aan ´e´en. Wel zien we dat ze dichtbij ´e´en liggen. We bewijzen dat door dit proces vaak genoeg te herhalen alle rij- en kolomsommen willekeurig dichtbij ´e´en komen. In dit geval hebben we na honderd iteraties (een eitje voor de computer) al vijftien decimalen precisie. Dat wil zeggen dat de de rij- en kolomsommen van R100· · · R1AC1· · · C100 gelijk zijn aan 1.000000000000000. Merk op dat
R100· · · R1 en C1· · · C100 nog steeds diagonaalmatrices zijn.
We implementeren ook andere algoritmes en bekijken hoe snel deze R en C vinden. Daar-naast bestuderen we hetzelfde idee voor unitaire matrices in plaats van niet-negatieve matri-ces. Tot slot komt ook hetzelfde idee maar dan met rij- en kolomproducten aan bod.
Bibliografie
[1] S. de Baerdemacker, L. Chen, A. de Vos en L. Yu (2016). The Birkhoff theorem for unitary matrices of arbitrary dimensions. Linear Algebra Appl., 514:151–164.
[2] S. de Baerdemacker en A. de Vos (2014) The decomposition of U (n) into XU (n) and ZU (n). The 44th International Symposium on Multiple-Valued Logic, Bremen, 173–177.
[3] S. de Baerdemacker en A. de Vos (2014). Scaling a unitary matrix. Open Syst. Inf. Dyn., 21:1450013. [4] S. de Baerdemacker en A. de Vos (2016). The Birkhoff theorem for unitary matrices of prime dimension.
Linear Algebra Appl., 493:455–468.
[5] R.B. Bapat en T.E.S. Raghavan (1997). Nonnegative matrices and applications. Cambridge University Press, New York.
[6] J.H. Brandts (2017). Inleiding numerieke wiskunde. Collegedictaat. Universiteit van Amsterdam. [7] M. F¨urer (2004). Quadratic convergence for scaling of matrices. Conference: Proceedings of the sixth
workshop on algorithm engineering and experiments and the first workshop on analytic algorithmics and combinatorics, New Orleans.
[8] M. Idel (2016). A review of matrix scaling and Sinkhorn’s normal form for matrices and positive maps. arXiv:1609.06349 [math.RA].
[9] M. Idel en M.M. Wolf (2015). Sinkhorn normal form for unitary matrices. Linear Algebra Appl., 471:76–84. [10] P. Knopp en R. Sinkhorn (1967). Concerning nonnegative matrices and doubly stochastic matrices. Pacific
J. Math., 21:343–348.
[11] T.L. Landis en B.N. Parlett (1982). Methods for scaling to doubly stochastic form. Linear Algebra Appl., 48:53–79.
[12] N. Linial, A. Samorodnitsky en A. Widgerson (2000). A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents. Combinatorica, 20:545–568.
[13] L. Mirsky en H. Perfect (1965). The distribution of positive elements in doubly stochastic matrices. The Journal of the London Mathematical Society, 40:689–698.
[14] U.G. Rothblum en S.A. Zenios (1992). Scalings of matrices satisfying line-product constraints and gene-ralizations. Linear Algebra Appl., 175:159-175.
[15] R. Sinkhorn (1963). A relationship between arbitrary positive matrices and doubly stochastic matrices. Ann. Math. Statist., 35:876–879.
A. BvN
Het bewijs van Stelling 3.9 in [5] is constructief en kan worden omgezet in een recursief algoritme. Elke recursie van dit algoritme reduceert een gegeven dubbelstochastische matrix tot een dubbelstochastische matrix met minstens ´e´en extra nulentry. Zo wordt uiteindelijk een permutatiematrix verkregen en termineert het algoritme. Iedere recursie vereist een positieve diagonaal. De lastigheid van het algoritme is het vinden van deze positieve diagonalen. Het begin van dit proces wordt in [5] onvolledig beschreven. Vandaar volgt hier een poging tot een heldere beschrijving.
Gegeven is een dubbelstochastische matrix M . Als M een permutatiematrix is, vormen de entries gelijk aan ´e´en een positieve diagonaal. Als M geen permutatiematrix is, bestaat er een entry ai1j1 ∈ (0, 1). Omdat M dubbelstochastisch is, bestaat er een j2 6= j1 zodat
ai1j2 ∈ (0, 1). Om dezelfde reden bestaat er een i2 6= i1 zodat ai2j2 ∈ (0, 1). Zo verkrijgt
men een rij x = (i1, i1, i2, i2, i3. . .) met rij-indices en een rij y = (j1, j2, j2, j3, j3, . . .) met
corresponderende kolomindices. We stoppen het vergroten van beide rijen simultaan direct nadat een rij- of kolomindex voor de derde keer voorkomt (in zijn eigen rij x of y). Schrijf ik
voor de laatste rij-index in x en jlvoor de laatste kolomindex in y (merk op dat l ∈ {k, k +1}).
De rijen x en y worden als volgt aangepast.
• Stel jl= j1. Verwijder uit zowel x als y de laatste entry.
• Stel jl 6= j1. Als ik drie keer voorkomt in x, dan geldt ik = ir voor een zekere r < k.
Verwijder de eerste 2r − 1 entries uit zowel x als y. Als jl drie keer voorkomt in y, dan
geldt jl= jr voor een zekere r < k. Verwijder de eerste 2r − 2 entries uit zowel x als y.
Laat t het aantal elementen van de nieuwe rij(en) x (en y) zijn. Zij P = ((x1, y1), . . . , (xt, yt))
en δ = mxsys = min{mij : (i, j) ∈ P}. Zonder verlies van algemeenheid mag men aannemen
dat s = 1. Beschouw anders de rij P = ((xs, ys), . . . , (xt, yt), (x1, y1), . . . , (xs−1, ys−1)).
De-finieer C zodat cij = 0 als (i, j) /∈ P en cxiyi = (−1)
i+1 voor i ∈ {1, . . . , t}. Vanaf hier kan
men hervatten met het volgen van de beschrijving in [5] (“Now the matrix A1 . . . ”), waar het
vinden van een positieve diagonaal in zijn geheel kort maar verhelderend wordt ge¨ıllustreerd in een voorbeeld.
Wanneer we het gehele algoritme implementeren in Matlab, doen zich problemen voor wegens eindige precisie. Om deze problemen te omzeilen kan men het probleem versimpelen en als input een niet-negatieve matrix met gehele entries nemen waarvan iedere rij en kolom optelt tot hetzelfde getal k ∈ N. De output is een combinatie van permutatiematrices waarvan de niet-negatieve scalairen optellen tot k. We geven een dergelijke implementatie.
function [scalars,perm] = BvN(T) n = size(T,1); I = eye(n); if nnz(T) == n scalars = T(find(T,1)); perm = sum((T~=0)*((1:n)’),2)’; return end A = T; while nnz(A) ~= n i = [1 1]; while nnz(A(i(1),:)) == 1 i = i + 1; end j = find(A(i(1),:),2);
while sum(i==i(end)) < 3 && sum(j==j(end)) < 3
i = [i, find([A(1:i(end)-1,j(end));0;A(i(end)+1:end,j(end))],1)]; j = [j, j(end)]; if sum(i==i(end)) == 3 || sum(j==j(end)) == 3 break end i = [i, i(end)]; j = [j, find([A(i(end),1:j(end)-1),0,A(i(end),j(end)+1:end)],1)]; end M = [i;j]; if j(end) == j(1) M = M(:,1:end-1); else if sum(i==i(end)) == 3 M = M(:,find(i==i(end),1)+1:end); else M = M(:,find(j==j(end),1)+1:end); end end delta = min(diag(A(M(1,:),M(2,:)))); s = find(diag(A(M(1,:),M(2,:)))==delta,1); M = [M(:,s:end),M(:,1:s-1)]; for k = 1:length(M)
A(M(1,k),M(2,k)) = A(M(1,k),M(2,k)) + ((-1)^k)*delta; end end pos_diag = sum((A~=0)*((1:n)’),2)’; labda = min(diag(T(1:n,pos_diag))); T = T - labda*I(pos_diag,:); [scalars1,perm1] = BvN(T); scalars = [labda;scalars1]; perm = [pos_diag;perm1];
B. FU
function [M,R,C] = FU(A,tol) n = size(A,2);
blok = zeros(n); B = [blok A; A’ blok]; [e,X] = deal(ones(2*n,1)); conv = max(abs(sum(B)-1)); while conv > tol
s = pinv(B+diag(sum(B,2)))*(e-sum(B,2)); while max(abs(exp(s).*(B*exp(s))-1))>max(abs(exp(s/2).*(B*exp(s/2))-1)) s = s/2; end X = X.*exp(s); B = diag(exp(s))*B*diag(exp(s)); conv = max(abs(sum(B)-1)); end R = diag(X(1:n)); C = diag(X(n+1:2*n)); M = B(1:n,n+1:2*n); R = R*C(1); C = C/C(1);
C. PB
function [R,C] = PB(A,tol) n = size(A,1); [R,C] = deal(eye(n)); [powers1,powers2] = deal([]); for k = 1:n m1 = nnz(A(:,k)); if m1 == 0 powers1 = [powers1 1]; else powers1 = [powers1 1/m1]; end m2 = nnz(A(k,:)); if m2 == 0 powers2 = [powers2 1]; else powers2 = [powers2 1/m2]; end end R1 = diag((1./prod((A+(A==0))’)).^powers2); A = R1*A; R = R1*R;while max(abs(prod(A+(A==0))- 1)) >= tol C1 = diag((1./prod(A+(A==0))).^powers1); A = A*C1; C = C*C1; R1 = diag((1./prod((A+(A==0))’)).^powers2); A = R1*A; R = R1*R; end
