4. Operationen auf strukturierten Matrizen BCCB-Matrizen

4. Operationen auf strukturierten Matrizen

BCCB-Matrizen

Zurück zum (3 × 3)-Beispiel: Sind die Randwerte periodisch, so ergibt sich etwa

b 2,1 = p 3,3 · x 1,3 + p 3,2 x 1,1 + p 3,1 x 1,2

+ p _2,3 · x _2,3 + p _2,2 x _2,1 + p _2,1 x _2,2 + p 1,3 · x 3,3 + p 1,2 x 3,1 + p 1,1 x 3,2 . Insgesamt:



























 b 1,1

b 2,1

b 3,1

b 1,2

b 2,2

b _3,2 b _1,3 b _2,3 b _3,3





























=





























p _2,2 p _1,2 p _3,2 p _2,1 p _1,1 p _3,1 p _2,3 p _1,3 p _3,3 p _3,2 p _2,2 p _1,2 p _3,1 p _2,1 p _1,1 p _3,3 p _2,3 p _1,3 p _1,2 p _3,2 p _2,2 p _1,1 p _3,1 p _2,1 p _1,3 p _3,3 p _2,3 p 2,3 p 1,3 p 3,3 p 2,2 p 1,2 p 3,2 p 2,1 p 1,1 p 3,1

p 3,3 p 2,3 p 1,3 p 3,2 p 2,2 p 1,2 p 3,1 p 2,1 p 1,1

p _1,3 p _3,3 p _2,3 p _1,2 p _3,2 p _2,2 p _1,1 p _3,1 p _2,1 p 2,1 p 1,1 p 3,1 p 2,3 p 1,3 p 3,3 p 2,2 p 1,2 p 3,2

p 3,1 p 2,1 p 1,1 p 3,3 p 2,3 p 1,3 p 3,2 p 2,2 p 1,2

p 1,1 p 3,1 p 2,1 p 1,3 p 3,3 p 2,3 p 1,2 p 3,2 p 2,2























































 x 1,1

x 2,1

x 3,1

x 1,2

x 2,2

x _3,2 x _1,3 x _2,3 x _3,3





























mit einer BCCB-Matrix.

Zirkulante Matrizen

Ist A ∈ C ^n×n zirkulant, d.h.

A =











a 0 a 1 · · · a n−1

a n−1 a 0 · · · a n−2

.. . .. . . .. .. . a ₁ a ₂ · · · a ₀











= p(S n ) :=

n−1

X

j=0

a j S _n ^j

(mit dem zirkulanten Shift S n ∈ R ^n×n ), so gilt

A = F _n ^H DF n , D = diag(λ ₀ , . . . , λ _n−1 ), mit der Matrix der diskreten Fourier-Transformation (DFT) (kurz:

Fourier-Matrix) (Jean Baptiste Joseph Fourier, 1768–1830) F _n = 1

√ n h ω ^−kj _n i

0≤k,j≤n−1 ∈ C ^n×n , ω _n := exp ^2πi _n

= cos ^2π _n

+ i sin ^2π _n

(i ² = −1). Die Eigenwerte von A sind gegeben durch

λ _j = p(ω ^j _n ), p(z) = a 0 + a 1 z + · · · + a n−1 z ⁿ⁻¹ , j = 0, . . . , n − 1.

Fourier-Matrix

Dies folgt aus der Spektralzerlegung

S n = F _n ^H D S F n , D S = diag(ω _n ⁰ , ω _n ¹ , . . . , ω _n ⁿ⁻¹ ) des zirkulanten Shifts.

Wichtige Eigenschaften der Fourier-Matrix:

F _n ^T = F n (aber F _n ^H 6= F n für n > 2).

F _n ^H F _n = I _n , d.h. F _n ist unitär.

F _n ⁻¹ = F _n ^H = F n . Bezeichnungen.

• Diskrete Fourier-Transformation (diskrete Fourier-Analyse):

x ˆ = F n x , d.h. x ˆ k = ^{√ 1}

n

P n

j =1 exp(−2πi(j − 1)(k − 1)/n) x j .

• Inverse diskrete Fourier-Transformation (diskrete Fourier-Synthese):

x = F _n ⁻¹ x ˆ , d.h. x _k = ^{√ 1}

n

P n

j=1 exp(2πi (j − 1)(k − 1)/n) ˆ x _j .

Fourier-Transformationen

Die „klassische“ Berechnung einer (inversen) diskreten Fourier-Transforma- tion (ein Matrix-Vektor-Produkt mit F n bzw. F _n ⁻¹ ) erfordert offenbar O (n ² ) komplexe Multiplikationen. Bei Anwendung der schnellen Fourier-Transforma- tion (FFT) reduziert sich dieser Aufwand auf O (n log ₂ n) (James William Cooley (∗1926) and John Wilder Tuckey (1915–2000): An algorithm for the machine calculation of complex Fourier series, Math. Comp. 19, 297–301 (1965)).

Diese Verbesserung kann nicht überbewertet werden: „It [the FFT] has

changed the face of science and engineering so much that it is not an

exaggeration to say that life as we know it would be very different without the

FFT.“ [Charles Van Loan, Computational Frameworks for the Fast Fourier

Transform, SIAM, Philadelphia 1992, p. ix].

Die Idee der FFT

Wir setzen (aus schreibtechnischen Gründen) im Folgenden

ζ _n = ¯ ω _n = ω ⁻¹ _n = exp

− 2πi n

= cos 2π

n

− i sin 2π

n

,

sodass F _n = [ζ _n ^kj ] _{0≤k,j≤n−1} .

Außerdem sei n gerade, z.B. n = 8: Mit ζ := ζ ₈ ist

F 8 =

























1 1 1 1 1 1 1 1

1 ζ ζ ² ζ ³ ζ ⁴ ζ ⁵ ζ ⁶ ζ ⁷ 1 ζ ² ζ ⁴ ζ ⁶ ζ ⁸ ζ ¹⁰ ζ ¹² ζ ¹⁴ 1 ζ ³ ζ ⁶ ζ ⁹ ζ ¹² ζ ¹⁵ ζ ¹⁸ ζ ²¹ 1 ζ ⁴ ζ ⁸ ζ ¹² ζ ¹⁶ ζ ²⁰ ζ ²⁴ ζ ²⁸ 1 ζ ⁵ ζ ¹⁰ ζ ¹⁵ ζ ²⁰ ζ ²⁵ ζ ³⁰ ζ ³⁵ 1 ζ ⁶ ζ ¹² ζ ¹⁸ ζ ²⁴ ζ ³⁰ ζ ³⁶ ζ ⁴² 1 ζ ⁷ ζ ¹⁴ ζ ²¹ ζ ²⁸ ζ ³⁵ ζ ⁴² ζ ⁴⁹

























(wir lassen den Vorfaktor 1/ √

n weg).

Die Idee der FFT

Wegen ζ ⁸ = 1, d.h. ζ ^j = ζ ^k , wenn j − k (ohne Rest) durch 8 teilbar ist, folgt

F 8 =

























1 1 1 1 1 1 1 1

1 ζ ζ ² ζ ³ ζ ⁴ ζ ⁵ ζ ⁶ ζ ⁷ 1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ³ ζ ⁶ ζ ζ ⁴ ζ ⁷ ζ ² ζ ⁵ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁵ ζ ² ζ ⁷ ζ ⁴ ζ ζ ⁶ ζ ³ 1 ζ ⁶ ζ ⁴ ζ ² 1 ζ ⁶ ζ ⁴ ζ ² 1 ζ ⁷ ζ ⁶ ζ ⁵ ζ ⁴ ζ ³ ζ ² ζ ¹























 .

Jetzt nummerieren wir die Zeilen von F ₈ um: zuerst werden die mit geradem

(0,2,4,6), danach die mit ungeradem Index (1,3,5,7) gezählt. Die zugehörige

Pemutationsmatrix wird mit P bezeichnet.

Die Idee der FFT

PF 8 =

























1 1 1 1 1 1 1 1

1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁶ ζ ⁴ ζ ² 1 ζ ⁶ ζ ⁴ ζ ² 1 ζ ζ ² ζ ³ ζ ⁴ ζ ⁵ ζ ⁶ ζ ⁷ 1 ζ ³ ζ ⁶ ζ ζ ⁴ ζ ⁷ ζ ² ζ ⁵ 1 ζ ⁵ ζ ² ζ ⁷ ζ ⁴ ζ ζ ⁶ ζ ³ 1 ζ ⁷ ζ ⁶ ζ ⁵ ζ ⁴ ζ ³ ζ ² ζ ¹

























=:

B 1,1 B 1,2

B 2,1 B 2,2

Wir untersuchen die einzelnen Blöcke: Wegen ζ = ζ ₈ ist ζ ² = ζ ₄ , d.h.

B _1,1 = B _1,2 =









1 1 1 1

1 ζ ₄ ζ ₄ ² ζ ₄ ³ 1 ζ ₄ ² ζ ₄ ⁴ ζ ₄ ⁶ 1 ζ ₄ ³ ζ ₄ ⁶ ζ ₄ ⁹









= F ₄ .

Die Idee der FFT

Aus den Spalten 0,1,2 bzw. 3 von B 2,1 „klammern“ wir ζ ⁰ , ζ ¹ , ζ ² bzw. ζ ³ „aus“:

B 2,1 =









1 1 1 1

1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁶ ζ ⁴ ζ ²

















1 0 0 0

0 ζ 0 0

0 0 ζ ² 0 0 0 0 ζ ³









= F 4 D 4 .

Analog:

B _2,2 =









1 1 1 1

1 ζ ² ζ ⁴ ζ ⁶ 1 ζ ⁴ 1 ζ ⁴ 1 ζ ⁶ ζ ⁴ ζ ²

















ζ ⁴ 0 0 0

0 ζ ⁵ 0 0 0 0 ζ ⁶ 0

0 0 0 ζ ⁷









= F ₄ (ζ ⁴ D ₄ ) = −F ₄ D ₄ .

Insgesamt erhalten wir

PF ₈ =

F ₄ F ₄ F ₄ D ₄ −F ₄ D ₄

=

F ₄ O O F ₄

I ₄ I ₄ D ₄ −D ₄

.

FFT

Fazit:

Seien n gerade, σ die folgende (even/odd) Permutation σ = [0, 2 . . . , n − 2, 1, 3, . . . , n − 1]

und P = P σ die zugehörige Permutationsmatrix.

Dann besitzt die zeilenpermutierte Fourier-Matrix F n die Zerlegung

PF _n =

F _n/2 F _n/2 F _n/2 D _n/2 −F _n/2 D _n/2

=

F _n/2 O O F _n/2

I _n/2 I _n/2 D _n/2 −D _n/2

.

Dabei bezeichnet D n/2 die Diagonalmatrix D _n/2 = diag

ζ _n ⁰ , ζ _n ¹ , . . . , ζ _n ^n/2−1

∈ C (n/2)×(n/2)

mit ζ _n = ¯ ω _n = exp(−2πi/n).

FFT

Berechne jetzt x ˆ = F n x für ein x ∈ C ⁿ (n gerade).

Wegen der obigen Zerlegung von F n unterteilen wir dies in zwei Schritte:

1. Reduktionsschritt: Berechne z =

I n/2 I n/2

D n/2 −D _n/2

x .

Für m = 8:

z ₀ = x ₀ + x ₄ , z ₁ = x ₁ + x ₅ , z ₂ = x ₂ + x ₆ , z ₃ = x ₃ + x ₇ z 4 = (x 0 − x 4 ), z 5 = (x 1 − x 5 )ζ _n , z 6 = (x 2 − x 6 )ζ _n ² , z 7 = (x 3 − x 7 )ζ _n ³

(n/2 komplexe Multiplikationen und n komplexe Additionen).

2. Berechne

F _n/2 z (0 : n/2 − 1) und F _n/2 z (n/2 : n − 1)

(zwei Fourier-Transformationen der Dimension n/2).

FFT

Ist n = 2 ^p eine Potenz von 2, so ist n/2 ebenfalls gerade und die beiden Fourier-Transformationen der Dimension n/2 können auf vier

Fourier-Transformationen der Dimension n/4 reduziert werden. Der Aufwand zur Reduktion beträgt 2 · n/4 = n/2 komplexe Multiplikationen (und

2 · n/2 = n komplexe Additionen). Diese Prozess wird solange fortgesetzt bis man eine Multiplikation mit F _n auf n Multiplikationen mit F ₁ = [1] reduziert hat (eine Multiplikation mit F ₁ erfordert offenbar keinen Aufwand). Dieses

Reduktionsverfahren heißt schnelle Fourier-Transformation (FFT = Fast Fourier Transform).

Damit ist gezeigt:

Zur Durchführung einer schnellen Fourier-Transformation der Ordnung n = 2 ^p

sind n

2 p = n

2 log ₂ (n) komplexe Multiplikationen

und n log ₂ (n) komplexe Additionen erforderlich.

Bit-reversal

Die konventionelle Berechnung einer Fourier-Transformation der Länge n = 2 ^p durch F _n x erfordert also 2 ^p+1 /p-mal mehr Multiplikationen als ihre Berechnung durch FFT. Wenn z.B. für p = 20 die FFT-Version eine Sekunde benötigt, so benötigt F _n x etwa 29 Stunden.

Verbleibendes Problem: Bestimmt man x ˆ = F n x durch FFT, so erhält man nicht x, sondern einen permutierte Version ˆ y = Q x ˆ mit einer

Permutationsmatrix Q ∈ R ^n×n .

Es gilt: Besitzt für n = 2 ^p der Index j ∈ {0, 1, . . . , n − 1} die Binärdarstellung j = b p−1 2 ^p−1 + · · · + b 2 2 ² + b 1 2 + b 0 =: [b p−1 . . . b 2 b 1 b 0 ] ₂ , und ist

r (j) := [b ₀ b ₁ b ₂ . . . b p−1 ] ₂ = b ₀ 2 ^p−1 + b ₁ 2 ^p−2 + b ₂ 2 ^p−3 + · · · + b p−1

(bit-reversal), dann gelten

x ˆ _j = y _{r (j)} und y _j = ˆ x _{r (i} ) .

FFT und Faltung

Bei periodischen Randbedingungen entspricht die eindimensionale Faltung der Multiplikation mit einer zirkulanten Matrix b = Ax .

Mit der Spektralzerlegung A = F _n ^H DF _n bietet sich folgendes Verfahren an:

Berechne y = F n x durch eine FFT.

Berechne z = Dy durch n Multiplikationen z _j = d _j,j y _j , j = 1, 2, . . . , n.

Berechne b = F _n ^H z durch eine FFT.

Die Eigenwerte von A (= die Diagonalelemente d j,j von D) ergeben sich aus F n A = DF n . Betrachet man die erste Spalte dieser Identität, so folgt

√ nF n a 1 = [d j,j ] _1≤j≤n ,

was eine weitere FFT erfordert.

Damit kann eine eindimensionale Faltung mit ( ³ ₂ log ₂ (n) + 1)n Multiplikationen

realisiert werden.

BCCB-Matrix: 3 × 3 Beispiel

Für unser lineares Modell mit periodischen Randbedingungen erhalten wir im 2D-Fall für P, X ∈ R ^3×3



























 b _1,1 b _2,1 b _3,1 b 1,2

b 2,2

b 3,2

b 1,3

b 2,3

b 3,3





























=





























p 2,2 p 1,2 p 3,2 p 2,1 p 1,1 p 3,1 p 2,3 p 1,3 p 3,3

p _3,2 p _2,2 p _1,2 p _3,1 p _2,1 p _1,1 p _3,3 p _2,3 p _1,3 p _1,2 p _3,2 p _2,2 p _1,1 p _3,1 p _2,1 p _1,3 p _3,3 p _2,3 p 2,3 p 1,3 p 3,3 p 2,2 p 1,2 p 3,2 p 2,1 p 1,1 p 3,1

p 3,3 p 2,3 p 1,3 p 3,2 p 2,2 p 1,2 p 3,1 p 2,1 p 1,1

p 1,3 p 3,3 p 2,3 p 1,2 p 3,2 p 2,2 p 1,1 p 3,1 p 2,1

p _2,1 p _1,1 p _3,1 p _2,3 p _1,3 p _3,3 p _2,2 p _1,2 p _3,2 p _3,1 p _2,1 p _1,1 p _3,3 p _2,3 p _1,3 p _3,2 p _2,2 p _1,2 p _1,1 p _3,1 p _2,1 p _1,3 p _3,3 p _2,3 p _1,2 p _3,2 p _2,2























































 x _1,1 x _2,1 x _3,1 x 1,2

x 2,2

x 3,2

x 1,3

x 2,3

x 3,3





























(3)

mit einer BCCB-Matrix.

BCCB-Matrizen: Spektralzerlegung

Eine BCCB-Matrix

A =











A ₁ A ₂ · · · A _m A _m A ₁ A m−1

.. . . .. .. . A 2 A 3 · · · A m











∈ R ^mn×mn ,

die aus zirkulanten Blöcken A j ∈ R ^n×n mit Spektralzerlegung A j = F _n ^H D ˜ j F n , D ˜ j = diag

˜ λ ^(j) ₁ , ˜ λ ^(j) ₂ , . . . , λ ˜ ^(j _n ⁾

, (j = 1, 2, . . . , m) zusammengesetzt ist, besitzt die Faktorisierung

A = (I _m ⊗ F _n ^H )











D ˜ ₁ D ˜ ₂ · · · D ˜ _m D ˜ m D ˜ 1 D ˜ m−1

.. . . .. .. . D ˜ ₂ D ˜ ₃ · · · D ˜ ₁











(I _m ⊗ F _n ). (4)

Permutationsmatrizen

Um die Spektralzerlegung von A zu bestimmen, ist eine Permutation der blockzirkulanten Matrix aus Diagonalblöcken erforderlich.

Ist

π =

π(1) π(2) . . . π(n)

eine Permutation der Indexmenge {1, 2, . . . , n} (eine bijektive Abbildung von {1, 2, . . . , n} auf {1, 2, . . . , n}), dann ist die zugehörige Permutationsmatrix P = P π = [p i,j ] ∈ R ^n×n durch

p i,j =

1, falls j = π(i ), 0, falls j 6= π(i ), definiert. Alternativ

P π =











 e ^T _π(1) e ^T _π(2)

.. . e ^T _π(n)













(e _j = j -ter Einheitsvektor).

Permutationsmatrizen: Eigenschaften

Beispiel.

π =

4 1 3 2 , P π =









0 0 0 1

1 0 0 0

0 0 1 0

0 1 0 0







 .

P π x =









0 0 0 1

1 0 0 0

0 0 1 0

0 1 0 0















 x 1

x 2

x 3

x 4









=







 x 4

x 1

x 3

x 2









=







 x _π(1) x _π(2) x π(3)

x π(4)







 .

Es gelten:

P _π ⁻¹ = P _π

−1

= P _π ^T .

Ist A = [a i ,j ] ∈ C ^n×n mit den Zeilen z ^T ₁ , z ^T ₂ , . . . , z ^T _n und Spalten

s 1 , s 2 , . . . , s n , so besitzt PA die Zeilen z ^T _π(1) , z _π(2) , . . . , z ^T _π(n) und AP ^T die

Spalten s _π(1) , s _π(2) , . . . , s _π(n) . Insbesondere ist PAP ^T = [a π(i),π(j ) ].

Transponierungspermutation

Wir betrachten speziell die sog. Transponierungspermutation π _× der Menge {1, 2, . . . , mn} bzw. die zugehörige Permutationsmatrix P × ∈ R ^mn×mn :

π × =

[ 1 1 + n 1 + 2n · · · 1 + (m − 1)n 2 2 + n 2 + 2n · · · 2 + (m − 1)n

· · · · · ·

n − 1 2n + 1 3n − 1 · · · mn − 1

n 2n 3n · · · mn ]

Zuerst kommen die Indizes j = 1 mod n, dann die Indizes j = 2 mod n, . . . , dann die Indizes j = n − 1 mod n, schließlich die Indizes j = n mod n (j = 0 mod n) (jeweils in in ihrer natürlichen Reihenfolge).

Soll die Abhänigkeit von m und n betont werden, so schreiben wir π _n,m bzw.

P _n,m statt π _× bzw. P × .

Transponierungspermutation: Beispiel

Für m = 4 und n = 3:

π _× =

1 4 7 10 2 5 8 11 3 6 9 12

und

P × = P _3,4 =









































1 0 0 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 0 0 1









































Transponierungspermutation: Eigenschaften

Ihren Namen verdankt die Transpositionspermutation der Eigenschaft vec(X ^T ) = P × vec(X ) für X ∈ C ^n×m .

Beispiel.

X =





11 12 13 14 21 22 23 24 31 32 33 34



 , X ^T =









11 21 31 12 22 32 13 23 33 14 24 34







 ,

P × vec(X ) =

11 21 31 12 22 32 13 23 33 14 24 34 T

=

11 12 13 14 21 22 23 24 31 32 33 34 T

,

also P × vec(X ) = vec(X ^T ).

Transponierungspermutation: Eigenschaften

P m,n = P _n,m ^T .

Seien A ∈ R ^n×n , B ∈ R ^m×m sowie P n,m ∈ R ^mn×mn , dann B ⊗ A = P _n,m ^T (A ⊗ B)P n,m .

Beispiel. Mit

F 4 = 1 2









1 1 1 1

1 −ı −1 ı

1 −1 1 −1

1 ı −1 −ı







 :

P _3,4 ^T (I ₃ ⊗ F ₄ )P _3,4 = P _3,4 ^T





F ₄ O O O F ₄ O O O F 4



 P _3,4 = 1 2









I ₃ I ₃ I ₃ I ₃ I ₃ −ıI ₃ −I ₃ ıI ₃ I 3 −I ₃ I 3 −I ₃ I 3 ıI 3 −I ₃ −ıI ₃









= F 4 ⊗ I 3 .

Transponierungspermutation: Anwendung auf Blockmatrizen

P _n,m vertauscht die Blockstruktur einer (m × m)-Blockmatrix mit (n × n)-Blöcken in die einer (n × n)-Blockmatrix mit (m × m)-Blöcken:

P

n,m















































˜ λ

⁽¹⁾₁

0 0 λ ˜

⁽²⁾₁

0 0 ˜ λ

⁽³⁾₁

0 0 λ ˜

⁽⁴⁾₁

0 0 0 ˜ λ

⁽¹⁾₂

0 0 λ ˜

⁽²⁾₂

0 0 ˜ λ

⁽³⁾₂

0 0 λ ˜

⁽⁴⁾₂

0 0 0 λ ˜

⁽¹⁾₃

0 0 λ ˜

⁽²⁾₃

0 0 ˜ λ

⁽³⁾₃

0 0 λ ˜

⁽⁴⁾₃

˜ λ

⁽⁴⁾₁

0 0 λ ˜

⁽¹⁾₁

0 0 ˜ λ

⁽²⁾₁

0 0 λ ˜

⁽³⁾₁

0 0 0 ˜ λ

⁽⁴⁾₂

0 0 λ ˜

⁽¹⁾₂

0 0 ˜ λ

⁽²⁾₂

0 0 λ ˜

⁽³⁾₂

0 0 0 λ ˜

⁽⁴⁾₃

0 0 λ ˜

⁽¹⁾₃

0 0 ˜ λ

⁽²⁾₃

0 0 λ ˜

⁽³⁾₃

˜ λ

⁽³⁾₁

0 0 λ ˜

⁽⁴⁾₁

0 0 ˜ λ

⁽¹⁾₁

0 0 λ ˜

⁽²⁾₁

0 0 0 ˜ λ

⁽³⁾₂

0 0 λ ˜

⁽⁴⁾₂

0 0 ˜ λ

⁽¹⁾₂

0 0 λ ˜

⁽²⁾₂

0 0 0 λ ˜

⁽³⁾₃

0 0 λ ˜

⁽⁴⁾₃

0 0 ˜ λ

⁽¹⁾₃

0 0 λ ˜

⁽¹⁾₃

˜ λ

⁽²⁾₁

0 0 λ ˜

⁽³⁾₁

0 0 ˜ λ

⁽⁴⁾₁

0 0 λ ˜

⁽¹⁾₁

0 0 0 ˜ λ

⁽²⁾₂

0 0 λ ˜

⁽³⁾₂

0 0 ˜ λ

⁽⁴⁾₂

0 0 λ ˜

⁽¹⁾₂

0 0 0 λ ˜

⁽²⁾₃

0 0 λ ˜

⁽³⁾₃

0 0 ˜ λ

⁽⁴⁾₃

0 0 λ ˜

⁽¹⁾₃















































P

^T_n,m

Transponierungspermutation: Anwendung auf Blockmatrizen

=

















































λ ˜ ⁽¹⁾ ₁ ˜ λ ⁽²⁾ ₁ ˜ λ ⁽³⁾ ₁ ˜ λ ⁽⁴⁾ ₁ 0 0 0 0 0 0 0 0 λ ˜ ⁽⁴⁾ ₁ ˜ λ ⁽¹⁾ ₁ ˜ λ ⁽²⁾ ₁ ˜ λ ⁽³⁾ ₁ 0 0 0 0 0 0 0 0 λ ˜ ⁽³⁾ ₁ ˜ λ ⁽⁴⁾ ₁ ˜ λ ⁽¹⁾ ₁ ˜ λ ⁽²⁾ ₁ 0 0 0 0 0 0 0 0 λ ˜ ⁽²⁾ ₁ ˜ λ ⁽³⁾ ₁ ˜ λ ⁽⁴⁾ ₁ ˜ λ ⁽¹⁾ ₁ 0 0 0 0 0 0 0 0 0 0 0 0 λ ˜ ⁽¹⁾ ₂ λ ˜ ⁽²⁾ ₂ λ ˜ ⁽³⁾ ₂ λ ˜ ⁽⁴⁾ ₂ 0 0 0 0 0 0 0 0 λ ˜ ⁽⁴⁾ ₂ λ ˜ ⁽¹⁾ ₂ λ ˜ ⁽²⁾ ₂ λ ˜ ⁽³⁾ ₂ 0 0 0 0 0 0 0 0 λ ˜ ⁽³⁾ ₂ λ ˜ ⁽⁴⁾ ₂ λ ˜ ⁽¹⁾ ₂ λ ˜ ⁽²⁾ ₂ 0 0 0 0 0 0 0 0 λ ˜ ⁽²⁾ ₂ λ ˜ ⁽³⁾ ₂ λ ˜ ⁽⁴⁾ ₂ λ ˜ ⁽¹⁾ ₂ 0 0 0 0 0 0 0 0 0 0 0 0 λ ˜ ⁽¹⁾ ₃ ˜ λ ⁽²⁾ ₃ ˜ λ ⁽³⁾ ₃ ˜ λ ⁽⁴⁾ ₃ 0 0 0 0 0 0 0 0 λ ˜ ⁽⁴⁾ ₃ ˜ λ ⁽¹⁾ ₃ ˜ λ ⁽²⁾ ₃ ˜ λ ⁽³⁾ ₃ 0 0 0 0 0 0 0 0 λ ˜ ⁽³⁾ ₃ ˜ λ ⁽⁴⁾ ₃ ˜ λ ⁽¹⁾ ₃ ˜ λ ⁽²⁾ ₃ 0 0 0 0 0 0 0 0 λ ˜ ⁽²⁾ ₃ ˜ λ ⁽³⁾ ₃ ˜ λ ⁽⁴⁾ ₃ ˜ λ ⁽¹⁾ ₃















































 .

Die BCCB-Struktur bleibt dabei erhalten.

BCCB-Matrizen: Spektralzerlegung

Mit der Transponierungspermutation P × geht der mittlere Faktor in (4) über in

P ×











D ˜ 1 D ˜ 2 · · · D ˜ m

D ˜ _n D ˜ ₁ D ˜ m−1

.. . . .. .. . D ˜ ₂ D ˜ ₃ · · · D ˜ ₁











P _× ^T = diag(˜ A ₁ , A ˜ ₂ , . . . , A ˜ _n ),

mit den zirkulanten Diagonalblöcken

A ˜ k =













λ ˜ ⁽¹⁾ _k ˜ λ ⁽²⁾ _k · · · λ ˜ ^(m) _k

˜ λ ^(m) _k ˜ λ ⁽¹⁾ _k ˜ λ ^(m−1) _k .. . . .. .. . λ ˜ ⁽²⁾ _k ˜ λ ⁽³⁾ _k · · · ˜ λ ⁽¹⁾ _k













∈ C ^m×m , k = 1, 2, . . . , n.

Diese werden ihrerseits durch die Fourier-Matrix F _m diagonalisiert:

A ˜ k = F _m ^H D ˆ k F m , D ˆ k = diag

λ ˆ ^(k) ₁ , λ ˆ ^(k) ₂ , . . . , ˆ λ ^(k) _m

, k = 1, 2, . . . , n.

BCCB-Matrizen: Spektralzerlegung

Zusammenfassung: Für eine BCCB-Matrix A gilt











A ₁ A ₂ · · · A _m A _m A ₁ A m−1

.. . . .. .. . A ₂ A ₃ · · · A ₁











= (I _m ⊗ F _n ^H )











D ˜ 1 D ˜ 2 · · · D ˜ m

D ˜ m D ˜ ₁ D ˜ m−1

.. . . .. .. . D ˜ ₂ D ˜ ₃ · · · D ˜ ₁











(I _m ⊗ F _n )

= (I m ⊗ F _n ^H )P _× ^T P ×











D ˜ 1 D ˜ 2 · · · D ˜ m

D ˜ _m D ˜ ₁ D ˜ m−1

.. . . .. .. . D ˜ 2 D ˜ 3 · · · D ˜ 1











P _× ^T P × (I m ⊗ F n )

= (I m ⊗ F _n ^H )P _× ^T











A ˜ ₁ O · · · O O A ˜ 2 O .. . . .. .. . O O · · · A ˜ _n











P × (I m ⊗ F n )

BCCB-Matrizen: Spektralzerlegung

= (I _m ⊗ F _n ^H )P _× ^T (I _n ⊗ F _m ^H )











D ˆ ₁ O · · · O O D ˆ ₂ O .. . . .. .. . O O · · · D ˆ n











(I _n ⊗ F _m )P × (I _m ⊗ F _n ).

Wir schreiben diese Spektralzerlegung als

(I m ⊗ F _n ^H )P _× ^T (I n ⊗ F _m ^H )P × P _× ^T











D ˆ 1 O · · · O O D ˆ 2 O .. . . .. .. . O O · · · D ˆ _n











P × P _× ^T (I n ⊗ F m )P × (I m ⊗ F n ),

weil P _× ^T (I n ⊗ F m )P × = F m ⊗ I n und folglich P _× ^T (I n ⊗ F m )P × (I m ⊗ F n ) = F m ⊗ F n

gelten.

P _× ^T diag

D ˆ 1 , D ˆ 2 , . . . , D ˆ n

P × =: D ist natürlich wieder eine Diagonalmatrix.

BCCB-Matrizen: Spektralzerlegung

Wir erhalten insgesamt die Spektralzerlegung:

A = U ^H DU (5)

mit

U = F _m ⊗ F _n

U ist die Matrixdarstellung der zweidimensionalen DFT.

Die Eigenwerte einer BCCB-Matrix erhält man durch Vergleich der ersten Spalten in der Identität UA = DU, d.h. Ua ₁ = Du ₁ , unter Beachtung von

u 1 = (F m ⊗ F n )e 1 = 1

√ mn e, e ^T = [1, 1, . . . , 1].

Also liefert die Anwendung der zweidimensionalen DFT auf die erste Spalte von A einen Vektor, dessen Komponenten die Eigenwerte von A sind (modulo des Faktors 1/ √

mn).

BCCB-Matrizen: Spezialfall PSF-Matrix

Im Spezialfall der zu einer PSF-Matrix bei periodischen Randbedingungen gehörigen BCCB-Matrix A ergibt sich die erste Spalte von A durch zirkulantes Shiften der PSF-Matrix.

Im Beispiel (3) ist die 2D-DFT anzuwenden auf

vec ⁻¹ (Ae 1 ) =





p 2,2 p 2,3 p 2,1

p 3,2 p 3,3 p 3,1

p 1,2 p 1,3 p 1,1



 .

Diese Matrix erhält man durch zirkulantes Shiften der Zeilen und Spalten der PSF-Matrix um einen Index in negativer Zeilen- und Spaltenrichtung:





p 1,1 p 1,2 p 1,3

p 2,1 p 2,2 p 2,3

p _3,1 p _3,2 p _3,3





circshift(−1,−1)

−−−−−−−−−−−−→





p 2,2 p 2,3 p 2,1

p 3,2 p 3,3 p 3,1

p _1,2 p _1,3 p _1,1



 .

Allgemein: Liegt die Mitte der PSF-Matrix an der Position (j _c , k _c ), so ist ein

zirkulanter Shift der Länge (1 − j _c , 1 − k _c ) erforderlich.

BCCB-Matrizen: Grundoperationen

Mit der schnellen Fourier-Transformation (FFT) lassen sich die DFT und damit Grundoperationen mit aus PSF-Matrizen aufgebauten BCCB-Matrizen effizient ausführen:

In M AT L AB bewirkt Y=fft2(X) bei einer m × n Matrix X die Anwendung der 2D-DFT auf X , d.h. die Multiplikation von vec(X) mit der (unnormierten) Matrix U aus (5). Die Umkehroperation bewirkt X=ifft2(Y).

Enthält der zweidimensionale Vektor c die Indizes des Zentrums der PSF-Matrix P ∈ R ^m×n , so berechnet der folgende M AT L AB -Befehl sämtliche Eigenwerte der zugehörigen BCCB-Matrix A:

D = fft2(circshift(P, 1-c));

BCCB-Matrizen: Grundoperationen

Das Produkt b = Ax eines Vektors x und einer BCCB-Matrix A mit Hilfe der FFT berechnet man in M AT L AB durch:

D = fft2(circshift(P, 1-c));

B = ifft2(D .* fft2(X));

B = real(B);

Bei der Multiplikation x = A ⁻¹ b mit der Inversen muss man wegen A ⁻¹ = U ^H D ⁻¹ U die Eigenwerte von A durch ihre Kehrwerte ersetzen:

D = fft2(circshift(P, 1-c));

X = ifft2(fft2(B) ./ D);

X = real(X);

Die Befehle B = real(B); bzw. X = real(X); unterdrücken durch

Rundungsfehler generierte Imaginärteile, die in exakter Arithmetik identisch

verschwinden.

BTTB+BTHB+BHTB+BHHB-Matrizen

Schreibt man in unserem (3 × 3)-Beispiel Spiegelungsrandbedingungen vor, so ergibt sich

A = A ₀ + A ₁ + A ₂ + A ₃ mit

A 0 =





























p _2,2 p _1,2 0 p _2,1 p _1,1 0 0 0 0 p _3,2 p _2,2 p _1,2 p _3,1 p _2,1 p _1,1 0 0 0 0 p _3,2 p _2,2 0 p _3,1 p _2,1 0 0 0 p 2,3 p 1,3 0 p 2,2 p 1,2 0 p 2,1 p 1,1 0 p 3,3 p 2,3 p 1,3 p 3,2 p 2,2 p 1,2 p 3,1 p 2,1 p 1,1

0 p _3,3 p _2,3 0 p _3,2 p _2,2 0 p _3,1 p _2,1 0 0 0 p _2,3 p _1,3 0 p _2,2 p _1,2 0 0 0 0 p 3,3 p 2,3 p 1,3 p 3,2 p 2,2 p 1,2

0 0 0 0 p 3,3 p 2,3 0 p 3,2 p 2,2





























(BTTB, Anteile der Pixel im ursprünglichen Bild),

BTTB+BTHB+BHTB+BHHB-Matrizen

A 1 =





























p _3,2 0 0 p _3,1 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 p _1,2 0 0 p _1,1 0 0 0

p 3,3 0 0 p 3,2 0 0 p 3,1 0 0

0 0 0 0 0 0 0 0 0

0 0 p 1,3 0 0 p 1,2 0 0 p 11

0 0 0 p _3,3 0 0 p _3,2 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 p _1,3 0 0 p _1,2





























(BTHB, Randbedingungen oben und unten),

BTTB+BTHB+BHTB+BHHB-Matrizen

A 2 =





























p _2,3 p _1,3 0 0 0 0 0 0 0 p _3,3 p _2,3 p _1,3 0 0 0 0 0 0

0 p _3,3 p _2,3 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 p _2,1 p _1,1 0

0 0 0 0 0 0 p _3,1 p _2,1 p _1,1

0 0 0 0 0 0 0 p _3,1 p _2,1





























(BHTB, Randbedingungen links und rechts)

BTTB+BTHB+BHTB+BHHB-Matrizen

und

A ₃ =





























p 3,3 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 p _1,3 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 p _3,1 0 0

0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 p _1,1





























(BHHB, Randbedingungen in NO-, SO-, SW- und NW-Richtung).

BTTB+BTHB+BHTB+BHHB-Matrizen

Diese etwas kompliziertere Matrixstruktur vereinfacht sich, wenn der nichtverschwindende Teil P ˜ der PSF-Matrix

P =





O O O O P ˜ O O O O



 , P ˜ ∈ R (2k−1)×(2k−1) mit Zentrum in (k , k ),

die folgende zweifache Symmetriebedingung (ZSB) erfüllt: (M AT L AB -Notation) P ˜ = flipud( ˜ P) = fliplr( ˜ P),

was wir im ab jetzt immer annehmen.

Besipiel: Die folgenden Matrizen erfüllen die ZSB





1 2 1

3 4 3

1 2 1



 ,









1 2 1 0

3 4 3 0

1 2 1 0

0 0 0 0







 ,













0 0 0 0

0 1 2 1

0 3 4 3

0 1 2 1

0 0 0 0













.

Ist etwa in unserem (3 × 3)-Beispiel

P ˜ =





α β α

γ δ γ

α β α



 ,

so folgt

A =





M ₁ + M ₂ M ₁ O M 1 M 2 M 1

O M ₁ M ₁ + M ₂





mit

M 1 :=





α + γ α 0

α γ α

0 α α + γ



 und M 2 :=





β + δ β 0

β δ β

0 β β + δ



 .

A ist blocksymmetrisch mit symmetrischen Blöcken.

BTTB+BTHB+BHTB+BHHB-Matrizen:

Spektralzerlegung

Im ZSB-Fall ist die Matrix A also symmetrisch, besitzt daher eine Spektralzerlegung

A = C ^T DC.

Die orthogonale Matrix C (deren Spalten die Eigenvektoren von A sind) repräsentiert die zweidimensionale diskrete Kosinustransformation (DCT).

Die 2D DCT kann (ähnlich wie die 2D DFT) durch schnelle Algorithmen realisiert werden; in der M AT L AB -Image Processing Toolbox stehen dazu die Befehle dct2 und idct2 (inverse Kosinustransformation) zur Verfügung.

Analog zum BCCB-Fall lassen sich die Eigenwerte von A bestimmen durch

λ _j = [Ca 1 ] _j c _j,1

mit der ersten Spalte a ₁ von A und dem j-ten Eintrag c _j,1 der ersten Spalte

von C.

BTTB+BTHB+BHTB+BHHB-Matrizen: Erste Spalte

Die erste Spalte von A lässt sich wieder durch einfache Manipulationen der PSF-Matrix P bestimmen. Sei

P =













p _1,1 p _1,2 p _1,3 p _1,4 p _1,5 p _2,1 p _2,2 p _2,3 p _2,4 p _2,5 p _3,1 p _3,2 p _3,3 p _3,4 p _3,5 p _4,1 p _4,2 p _4,3 p _4,4 p _4,5 p _5,1 p _5,2 p _5,3 p _5,4 p _5,5













.

BTTB+BTHB+BHTB+BHHB-Matrizen: Erste Spalte

Bei Spiegelungs-Randbedingungen setzt sich (unter wiederholter Ausnutzung der ZSB) die erste Spalte der BTTP+BTHB+BHTB+BHHB-Matrix A ∈ R ^25×25 zusammen aus den Anteilen













p 3,3 p 3,4 p 3,5 0 0 p 4,3 p 4,4 p 4,5 0 0 p _5,3 p _5,4 p _5,5 0 0

0 0 0 0 0

0 0 0 0 0











 +













p 3,4 p 3,5 0 0 0 p 4,4 p 4,5 0 0 0 p _5,4 p _5,5 0 0 0

0 0 0 0 0

0 0 0 0 0













+













p _4,3 p _4,4 p _4,5 0 0 p _5,3 p _5,4 p _5,5 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0











 +













p _4,4 p _4,5 0 0 0 p _5,4 p _5,5 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0











 .

Hierbei entsprechen die Anteile Pixeln im ursprünglichen Bildgebiet, jenseits

des linken Randes, jenseits des oberen Randes bzw. links oberhalb der linken

oberen Bildecke.

BTTB+BTHB+BHTB+BHHB-Matrizen: Erste Spalte

Allgemein konstruiert man eine Matrix aus den Einträgen der ersten Spalte von A wie folgt:

Sei P ∈ R (2k−1)×(2k−1) mit Zentrum an Position (k , k ) (hier k = 3).

Definiere die Shiftmatrizen Z ₁ und Z ₂ durch die M AT L AB -Befehle

Z ₁ = diag(ones(k,1), k-1) =













0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

0 0 0 0 0

0 0 0 0 0













;

Z 2 = diag(ones(k-1,1), k) =













0 0 0 1 0

0 0 0 0 1

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0











 .

Dann enthält die Matrix Z ₁ PZ ₁ ^T + Z ₁ PZ ₂ ^T + Z ₂ PZ ₁ ^T + Z ₂ PZ ₂ ^T die erste

Spalte von A.

BTTB+BTHB+BHTB+BHHB-Matrizen:

Spektralzerlegung

Diese Konstruktion lässt sich leicht auf PSF-Matrizen mit beliebiger Größe und beliebigem Zentrum verallgemeinern. Die HNO M AT L AB -Routine dctshift, bewerkstelligt diese Manipulation (analog zu circshift).

Die Eigenwerte lassen sich somit berechnen durch:

e1 = zeros(size(P)); e1(1,1) = 1;

D = dct2(dctshift(P,c)) ./ dct2(e1);

Multiplikation eines Vektors x = vec(X ) mit A bzw. A ⁻¹ erfolgt entsprechend durch

B = idct2( D.* dct2(X));

bzw. durch

B = idct2( dct2(X) ./ D);

BTTB-Matrizen: Toeplitz-Matrizen

Die Spektralzerlegung einer BTTB-Matrix ist im Gegensatz zu

BCCB-Matrizen oder BTTB+BTHB+BHTB+BHHB-Matrizen nicht auf einfache Weise mittels FFT möglich. Das Matrix-Vektor-Produkt lässt sich jedoch nach einer Einbettung in eine zirkulante Matrix mittels FFT effizient durchführen.

Zunächst betrachten wir Toeplitz-Matrizen

T = toeplitz(t 1−n , . . . , t −1 , t 0 , t 1 , . . . , t n−1 ) :=











t 0 t −1 . . . t 1−n

t 1 t 0 . . . t 2−n

.. . . .. ... .. . t n−1 t n−2 . . . t 0









 .

Toeplitz-Matrizen gehören zur (größeren) Klasse der persymmetrischen Matrizen, das sind Matrizen B = [b _i,j ] ∈ C ^n×n , die die Symmetrieeigenschaft b _i,j = b n−j+1,n−i+1 bzw. B = EB ^T E besitzen mit der

„Umkehrpermutationsmatrix“

E = [e n , e n−1 , . . . , e 2 , e 1 ].

(Beachte E = E ^T = E ⁻¹ .)

BTTB-Matrizen: Zirkulante Einbettung von Toeplitz-Matrizen

Bei gegebener Toeplitz-Matrix T = toeplitz(t 1−n , . . . , t −1 , t 0 , t 1 , . . . , t n−1 ) ∈ C ^n×n wird durch

T _ext :=

T S S T

∈ C ^2n×2n

mit S := toeplitz(t ₁ , . . . , t n−1 , 0, t 1−n , . . . , t −1 ) ∈ C ^n×n eine zirkulante Matrix definiert.

Beispiel: (n = 3)

T =





t

0

t

−1

t

−2

t

1

t

0

t

−1

t

2

t

1

t

0



 , S =





0 t

2

t

1

t

−2

0 t

2

t

−1

t

−2

0



 ,

T

ext

=

















t

0

t

−1

t

−2

0 t

2

t

1

t

1

t

0

t

−1

t

−2

0 t

2

t

2

t

1

t

0

t

−1

t

−2

0

0 t

2

t

1

t

0

t

−1

t

−2

t

−2

0 t

2

t

1

t

0

t

−1

t

−1

t

−2

0 t

2

t

1

t

0

















.

BTTB-Matrizen

Zirkulante Einbettung von Toeplitz-Matrizen

Zur Multiplikation von T mit einem Vektor v ∈ C ⁿ berechnet man das Produkt

T _ext v

0

= T S

S T v 0

= T v

Sv

mittles FFT und erhält T v in den ersten n Komponenten des Ergebnisvektors.

Im Block-Fall verfährt man entsprechend.

Separable PSF-Matrizen

Wir betrachten nochmals den separablen Fall

A = A r ⊗ A c ∈ R ^mn×mn , A r ∈ R ^n×n , A c ∈ R ^m×m .

Dieser Fall liegt vor, wenn P = cr ^T gilt für zwei Vektoren c ∈ R ^m und r ∈ R ⁿ . Beispiel:

P =



 c 1

c 2

c 3





r 1 r 2 r 3

=





c 1 r 1 c 1 r 2 c 1 r 3

c 2 r 1 c 2 r 2 c 2 r 3

c 3 r 1 c 3 r 2 c 3 r 3





führt bei Null-Randbedingungen auf A = A r ⊗ A c mit

A r =



 r 2 r 1

r 3 r 2 r 1

r 3 r 2



 , A c =



 c 2 c 1

c 3 c 2 c 1

c 3 c 2



 .

Separable PSF-Matrizen: Konstruktion der Kronecker-Faktoren

Ist eine PSF-Matrix separabel (oft explizit an der PSF erkennbar), so können die Vektoren r und c durch Berechnung der dominanten Singulärvektoren bestimmt werden.

In M AT L AB :

[u, s, v] = svds(P, 1);

c = sqrt(s) * u;

r = sqrt(s) * v;

Zur Sicherheit sollte geprüft werden, ob P überhaupt separabel ist, d.h., ob

der Rang von P numerisch gleich Eins ist. (Wie?)

Separable PSF-Matrizen: Grundoperationen

Matrix-Vektor-Multiplikation:

Ist A = A r ⊗ A c , b = vec(B), x = vec(X ), so gilt vec(B) = b = Ax = vec(A c XA ^T _r ).

In M AT L AB : B = Ac * X * Ar’.

Lösung linearer Gleichungssysteme mit A:

Wegen (A _r ⊗ A _c ) ⁻¹ = A ⁻¹ _r ⊗ A ⁻¹ _c gilt

A ⁻¹ b = (A ⁻¹ _r ⊗ A ⁻¹ _c )b = vec(A ⁻¹ _c BA ^−T _r ).

In M AT L AB : X = Ac \ B / Ar’.

Separable PSF-Matrizen: Grundoperationen

Singulärwertzerlegung (SVD):

Gelten A _r = U _r Σ _r V _r ^T und A _c = U _c Σ _c V _c ^T , so folgt A = A _r ⊗ A _c = (U _r Σ _r V _r ^T ) ⊗ (U _c Σ _c V _c ^T )

= (U _r ⊗ U _c )(Σ _r ⊗ Σ _c )(V _r ⊗ V _c ) ^T .

Dies ist bis auf die Reihenfolge der Singulärwerte und -vektoren die SVD von A.

Die (naive) Lösung V Σ ⁻¹ U ^T b des LGS Ax = b mit b = vec(B) ist daher

x = vec(X ) mit X = A ⁻¹ _c BA ^−T _r = (V c Σ ⁻¹ _c U _c ^T )B(U r Σ ⁻¹ _r V _r ^T ).

Separable PSF-Matrizen: Grundoperationen

Singulärwertzerlegung (Fortsetzung):

Der mittlere Ausdruck Σ ⁻¹ _c U _c ^T BU r Σ ⁻¹ _r lässt sich einfach durch die komponentenweise Division von M AT L AB berechnen gemäß

U _c ^T BU r ./ S.

Hier enthält S ∈ R ^m×n die Singulärwerte von A, und zwar so angeordnet, dass vec(S) die Diagonalelemente von Σ _r ⊗ Σ _c enthält.

In M AT L AB :

[Uc, Sc, Vc] = svd(Ac);

[Ur, Sr, Vr] = svd(Ar);

S = diag(Sc) * diag(Sr)’;

X = Vc * ( ( Uc’ * B * Ur) ./ S) * Vr’;