Transcript Fiber Optics Design and Solving symmetric Banded systems

Fiber Optics Design and
Solving Symmetric Banded
Systems
Linda Kaufman
William Paterson University
Dispersion Compensation fibers
Signal degradation and Restoration
High Speed Optical Communication Segment
EDFA
Transmission (~80km)
Amplification
Dispersion
Compensation
Fiber Optics Design and
Solving symmetric Banded
systems
Linda Kaufman
William Paterson University
Marcuse model for fiber wrapped
around spool
 For a radially symmetric fiber
where r is the radius from the center of the fiber, R is the
radius of the spool which leads to a matrix problem of the
structure.
b’s and
c’s are
0(r/R)
Objective: Create a algorithm to factor a symmetric
indefinite banded matrix A
An n x n matrix A is symmetric if
ajk = akj
 1 10 20 


 10 30 50 
 20 50  70


A matrix is indefinite if any of these hold
a. eigenvalues are not necessarily all
positive or all negative
b. One cannot factor A into LLT
1 -5
-5 2
is indefinite- determinant is negative
m=2
x x x
x x x x
0
x x x x x
A has band width 2m+1 if
ajk = 0 for |k-j| >m
x x x x x
x x x x x
0
x x x x
x x x
Uses
 (1)Solve a system of equations
Ax=b
If A = MDMT , D is block diagonal, one solves
Mz=b
Dy =z
MTx =y
(2)Find the inertia of a system- the number of
positive and negative eigenvalues of a matrix
If A = MDMT, The inertia of A is the inertia of D.
Given a matrix B, the inertia of a matrix A = B-cI
is the number of eigenvalues greater, less than and equal
to c.
Competition
 Ignore symmetry- use Gaussian elimination- does not
give inertia info- 0(nm2) time-(n is size of matrix, m is
bandwidth)
 Band reduction to tridiagonal (Schwarz,Bischof,
Lang, Sun, Kaufman) followed by Bunch for
tridiagonal-0(n2m)
 Snap Back-Irony and Toledo-Cerfacs-2004, 0(nm2)
time generally faster than GE but twice space
 Bunch – Kaufman for general matrices and hope
bandwidth does not grow as Jones and Patrick
noticed
Difficulties
Consider the linear system
.0001x + y = 1.0
x
+ y = 2.0
Gaussian elimination- 3 digit arithmetic
.0001x + y = 1.0
10000y = 10000
Giving y = 1, x =0
But true solution is about
x =1.001, y =.999
If interchange first and second rows and columns before
Gaussian elimination get
x + y = 2.0
y = 2.0 -1.0= 1.0
So x = 1.0- a bit better
Gaussian elimination with partial pivoting to prevent division
by zero and unbounded elemental growth
1 1 10
1 8 4 6
10 4 3 5 8
6 5 7 1 3
8 1 6 21
3 2 54
1 4 7
Eliminating first
column yields
Unsymmetric
pivoting yields
10 4
0 x
0 x
6
3
x
x
5
8
5 8
x f
f f
7 1 3
1 6 21
3 2 54
1 4 7
10 4 3 5 8
1 8 4 6
1 1 10
6 5 7 1 3
8 1 6 21
3 2 54
1 4 7
5 diagonal
becomes 7
diagonal
In general
2k+1 diagonal
becomes 3k+1
diagonal
Symmetric pivoting

But to maintain symmetry one must pivot rows and columns simultaneously
What if matrix is
0
1
?
1
0
Interchanging first row and column does not help
If matrix is
0
1
1000
1
0
10
1000
10
0
Pivot it to
0
1000
1000
0
1
10
1
10
0
And use top 2 x 2 as a “pivot”
But pivoting tends to upset the band structure!!
Bunch -Kaufman for symmetric
indefinite non banded
Partition A as
Where D is either 1 x 1 or 2 x 2
and B’ = B – Y D-1 YT
Choice of dimension of D depends on
magnitude of a11 versus other elements
If D is 2 x 2, det(D)<0
2 x 2 vs 1 x 1 for nonbanded symmetric system
 4 2 1
2 3 3
1 3 5
1 10 5
10 200 3
5 3 8
1 10 5
10 50 3
5 3 8
1x1
1 0 0
0 2 2.5
0 2.5 4.75
1x1
1 0 0
0 100 -47
0 -47 -17
2x2
1 10 0
10 50 0
0 0 25.3
Bandwidth spread with BunchKaufman on banded matrix
Banded algorithm based on B-K
 1) Let c = |ar 1 | = max in abs. in col. 1
 2) If |a11 | >= w c, use a 1 x 1 pivot. Here w is a scalar to balance
element growth, like 1/3
 Else
3)Let f= max element in abs. in column r
4) If w c*c <= |a11 | f, use a 1 x 1 pivot
Else
5)interchange the rth and second rows and
columns of A
6) perform a 2 by 2 pivot
7) fix it up if elements stick out beyond the
original band width
Never pivot with 1 x 1
Fix up algorithm
Worst case r =m+1, what happens in pivoting
x x x x x x
x x x x x x x a b c d
x x x x x x x x
x x x x x x x x x
x x x x x x x x x x
x x x x x x x
x x x x x x x x x x x
a x x x x x x x x x x
b
x x x x x x x x xx
c
x x x x x x x xx
d
x x x x x x xx
x x x x x x x x
x x x x x x x
x x x x x x
Partition A as
Reset B’ = B – Y D-1 YT
Let Z = D-1 YT =
x x x x x x x x x
p q r s x x x x x .
Then B’ looks like
x x x
x x x
x x x
x x x
x x x
x x x
bp x x
cp cq x
dp dq dr
x x
x x
x x
x x
x x
as x
bs x
cs x
ds x
x
x bp cp dp
x x cq dq
x x cr dr
as bs cs ds x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
If we don’t remove elements outside the band,
the bandwidth could explode to the full matrix.
Because of structure eliminating 1
element gets rid of column
x
x
x
x
x
x
x x x
x x x
x x x
x x x
x x x
x x x
x x bt
cq cr ct
dq dr dt
x
x
x
x
x
x
x
x
x
x
x
x x
x x
x bt
x x
x x
x x
x x
x x
x x
x x
x
cq dq
cr dr
ct dt
x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
Continue with eliminating another
element
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
x
u
x
x
x
x
x
x
x
x
m
x
x
x
x
x
x
x
x
x
x
x
x x
x x
x x
x x
x x
x x
x x
x x
x x
x x
x
x
x
x
x
x
x
x
x
x
x
u
m
x x
x x x x
x x x x
x x x x
x x x x
x x x x
x x x x
x x x x
Now eliminate u to restore bandwidth
In practice one would just use the
elements in the first part of Z to
determine the Givens/stabilized
elementary transformations and never
bother to actually form the bulge.
Thus there is never any need to
generate the elements outside the
original bandwidth.
Alternative Formulation
Partition A as
Where D is 2 x 2
Let Z = D-1 YT
Reset B’ = QT(B – Y Z)Q= QTB Q-HG
Where H=QTY and G =ZQ
Therefore do “retraction” followed by rank 2
correction.
Recall H looks like
h1 h2
0 h3
Where the “0” has length r-1 and the whole H
has length m+r-1
Alternative-2
Q is created such that G =ZQ has the form
 G1 G2 G3


G=
0
G
4
G
5


where 0 has r-3 elements and G5 is a multiple of G3
1)Explicitly use 0’s in rank-2 correction
Thus correction has form below where numbers indicate the rank
of the correction
1 1 0
(r-3 rows)
1 2 1
(m-r+2 rows)
0 1 1
(r-1 rows)
2) Store Q info in place of 0’s and G3 to reduce space to (2m+1)n
3) Reduce solve time for each 2 x 2 from 4m+6r mults to 4m+2rIn worst case, 3mn vs. 5mn
Properties
 Space- (2m+1)n-in order to save transformations





compared to (3m+1)n for unsymmetric G.E.
Never pivots for positive definite or 1 x 1
Decrease operations by not applying second column
of pivot when these will be undone
Operation count depends on types of transformations
Elementary –between m2n/2 and 5 m2n/4
Compared with between m2n and 2m2n for G.E.
Elemental growth
 Let a = max element of matrix
 For 1 x 1, Ajk = Ajk – A1k Aj1 /A11 taking norms
element growth is a(1+1/w)
 For 2 x 2 if no retraction, it turns out that
element growth after 1 step is a(1 + (3+w)/(1w))
 For 2 x 2 if retraction, it turns out that element
growth is (4+8/(1-w))a
 2 steps of 1 x 1 = 1 step of retraction when
w=1/3, giving growth bound of 4n-1
Comparison using Atlas Blas on Sun
n=1000, m=100
problem
diagonal
offdiagonal
outer diagonal
1
2
3
4
1
1
1
10*|i-j|
1
100
10000
10*|i-j|
100
10
10
1
pos. eigenvalues
/2x2/average r
1000/0/0/
502/161/66
500/500/89
493/397/57
Problem
1
2
3
4
Time-retraction
(ms)elem/orthog
98/98
161/240
296/520
220/220
Time-Snapback
elem/orthog
140/140
550/590
1640/
2200
1290/1370
Time-Lapack-tf2/trf 195/136
395/235
395/235
395/235
Error-retraction
3e-15
/3e-15
1.6e-13
/1.3e-13
3e-13/
9e-14
2e-11
/5e-12
Error-snapback
5e-15
/5e-15
4e-11
1e-13
2.2e-4
4e-12
3e-5
2e-11
Error Lapack-tf2
7e-15
1e-15
6e-15
1e-12
Retraction vs. Lapack
Lapack vs. retraction, n=2000, on random matrices,
with about 1000 + eigenvalues
30
25
time
20
dgbtrf
retraction
15
dgbtf2
10
5
0
0
100
200
300
400
half bandwidth
500
600
Elementary vs. orthogonal on 2 random matrices, n=1000,
m=100
Elementary
Orthogonal
Time
150
160
2 x2’s
204
179
error
3.04e-11
1.10e-11
growth
146.4
1292.1
Time
160
160
2 x2’s
221
178
error
1.00e-12
3e-12
growth
130.5
1762.8
I thought orthogonal would be more accurate, have
less elemental growth, and take much more time-but
these assumptions were wrong for random matrices
Future
 Could do pivoting with 1 by 1 at price of
space and time- needs to be implemented
 Adaptations for small bandwidth as in optical
fiber code.
 Condition estimator