Fiber Optics Design and Solving Symmetric Banded Systems

1
Fiber Optics Design and Solving Symmetric Banded
Systems

• Linda Kaufman
• William Paterson University

2
Dispersion Compensation fibers

• Signal degradation and Restoration

3
Fiber Optics Design and Solving symmetric Banded
systems

• Linda Kaufman
• William Paterson University

4
(No Transcript)
5
Fiber Design

6
Model-Maxwells equation in cylindrical
coordinates

• For a radially symmetric fiber

r is the radius from the center of the fiber, m
and ? are known and ß and x are unknown. We have
a pde-eigenvalue problem. Want to find the
relative index profile ? such that the
dispersion
has particular values for specific values of
where
?
Note ? gives the width and dopant concentration
of the layers of the fiber and is usually defined
by about 20 parameters
7
Typical refractive profile and parameterization
8
Discretize the differential equation
9
Discretization
Leads to a nonlinear eigenvalue problem of the
form
The function s() takes into consideration the
boundary condition. If wrapped around a spool
Marcuse model leads to 7 diagonal problem
10
Each function evaluation of an optimizer to
determine the parameters of ? requires

• For several values of m
• For about 20 values of ?
• Determine the grid and make A and M
• Solve the nonlinear eigenvalue problem for the
  positive eigenvalues
• Might involve several solutions of a linear
  eigenvalue problem
• Determine the dispersion, dispersion slope and
  the effective index- an integral of the
  eigenvectors. The dispersion is given as

where ??2pc, c is the speed of light
11
Difficulties

• Numerical difficulties means that numerical
  differentiation does not work and must supply the
  gradient of the dispersion- yikes third
  derivatives!
• The eigenvalue problem is slightly nonlinear and
  gives real problems when ß is near zero
• Need accurate model so the dispersion has to be
  computed analytically- means differentiating an
  eigenvalue
• Want to move onto a fiber wrapped around a spool
  which leads to 7 diagonal problem and solving a
  symmetric system-where the outermost diagonal has
  the larger elements and for about 2000
  eigenvalues 20 are positive.

12
Definitions symmetric,indefinite, banded
An n x n matrix A is symmetric if ajk akj
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 A has band width 2m1 if ajk 0 for
k-j gtm
m2 x x x x x x x 0 x x x x x x x
x x x x x x x x 0 x x
x x x x x
13
Applications

• Finding intermediate eigenvalues of a symmetric
  problem that comes from a partial differential
  equation on a regular grid
• Optical fiber design of a fiber wrapped around a
  spool using Marcuse model
• Designing cavity of a linear accelerator
• Shift and invert Lanczos
• KKT equations-minimizing 1/2xtPx xtq such that
  Ax b gives the optimality conditions

14
Uses

• (1)Solve a system of equations Axb
• If A MDMT , D is block diagonal, one solves
• Mzb
• 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.

15
Competition

• Ignore symmetry- use Gaussian elimination- does
  not give inertia info- 0(nm2) time
• 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

16
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.001y .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
17
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 2 1 3
2 5 4 1 4 7
Unsymmetric pivoting yields
10 4 3 5 8 1 8 4 6 1 1 10 6 5 7
1 3 8 1 6 2 1 3 2 5
4 1 4 7
10 4 3 5 8 0 x x x f 0 x x f f
6 5 7 1 3 8 1 6 2 1
3 2 5 4 1 4 7
Eliminating first column yields
5 diagonal becomes 7 diagonal In general 2k1
diagonal becomes 3k1 diagonal
18
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 1
• 1000 0 10
• 1 10 0
• And use top 2 x 2 as a pivot
• But pivoting tends to upset the band structure!!

19
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)lt0

20
2 x 2 vs 1 x 1 for nonbanded symmetric system

• 4 2 1 1 0 0
• 2 3 3 1 x 1
  0 2 2.5
• 1 3 5 0 2.5 4.75
• 1 10 5 1 0 0
• 10 200 3 1 x 1 0 100 -47
• 5 3 8 0 -47 -17
• 1 10 5 1 10 0
• 10 50 3 2 x 2 10 50
  0
• 5 3 8 0 0 25.3

21
Banded algorithm based on B-K

• 1) Let c ar 1 max in abs. in col. 1
• 2) If a11 gt 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 cc lt 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

22
Fix up algorithm
Worst case r m1, 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 x x c x x x x x x
x x x 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
23
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 bp cp dp x x x x x x x
cq dq x x x x x x x cr dr x x
x x x as bs cs ds x x x x x x
x x x x x x x x x x as x x
x x x x x x bp x x bs x x x
x x x x x cp cq x cs x x x x
x x x x dp dq dr 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 If we dont remove elements
outside the band, the bandwidth could explode to
the full matrix.
24
Because of structure eliminating 1 element gets
rid of column
x x x x x x x x x x x x x
cq dq x x x x x x x cr dr x x
x x x x bt 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 bt x x x x x x x
x cq cr ct x x x x x x x x
dq dr 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
25
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 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 Now eliminate u to restore
bandwidth
26
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.
27
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 HQTY 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 mr-1

28
Alternative-2
Q is created such that G ZQ has the form G
where 0 has r-3 elements and G5 is a multiple
of G3 1)Explicitly use 0s 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-r2 rows) 0 1
1 (r-1 rows) 2) Store Q info in place of
0s and G3 to reduce space to (2m1)n 3) Reduce
solve time for each 2 x 2 from 4m6r mults to
4m2r- In worst case, 3mn vs. 5mn
29
Properties

• Space- (2m1)n-in order to save transformations
  compared to (3m1)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.

30
Elemental growth

• Let a max element of matrix
• For 1 x 1, Ajk Ajk A1k Aj1 /A11 taking
  norms element growth is a(11/w)
• For 2 x 2 if no retraction, it turns out that
  element growth after 1 step is a(1 (3w)/(1-w))
• For 2 x 2 if retraction, it turns out that
  element growth is (48/(1-w))a
• 2 steps of 1 x 1 1 step of retraction when
  w1/3, giving growth bound of 4n-1

31
Comparison using Atlas Blas on Sun
Problem 1 2 3 4
Time-retraction (ms) 98 161 341 244
Time-Snapback4 140 590 2200 1370
Time-Snap3 140 550 1640 1290
Time-Lapack-tf2/trf 195/136 395/235 395/235 395/235
Error-retraction 3e-15 3e-13 3e-13 2e-11
Error-snapback4 5e-15 1e-13 4e-12 2e-11
Error-snap3 5e-15 4e-11 2.2e-4 3e-5
Error Lapack-tf2 7e-15 1e-15 6e-15 1e-12
32
Comparison with snapback-reference Blas- n500,
m40
problem diagonal Off diagonal except outer Outer diagonal
1 1 .01 .01
2 1 .01 10
3 1 .01 1000
4 1 Increases by 10
problem Time retract Time Snap3 Time snap4 Error retract Error-snap3 Error- snap4
1 50 60 60 2e-15 2e-15 2e-15
2 70 140 160 4e-14 1e-12 4e-14
3 110 430 701 6e-13 1e-9 4e-12
4 90 300 260 8e-13 1e-11 1e-12
33
Comparison on random examples-Atlas blas
N2000,m100 N1000,m100 N1000,m200
No. of positive eigenvalues 994 497 497
2 x2s 432 223 164
Average r 48 45 91
Time-retract 380 180 532
error 3.8e-12 3e-12 5e-11
Time-tf2 750 358 1306
Error-tf2 8e-13 5e-13 1e-11
34
Future

• Could do pivoting with 1 by 1 at price of space
  and time- needs to be implemented
• Compare time/ element growth with Givens and
  stabilized elementary
• Comparison with Toledo/Irony code and Lapack on
  bigger machines
• Adaptations for small bandwidth as in optical
  fiber code.

35
The Nonlinear eigenvalue problem

• First attempt

Very slow when eigenvalues near 0.
36
Solve the linear eigenvalue problem using
safeguarded Rayleigh quotient
37
Second attempt Nonlinear Rayleigh quotient

• Try to minimize

Solution is at
38
Safeguarded Nonlinear Rayleigh quotient
39
Third Attempt Newton Approach
Find root of
Iterate
Where x is an eigenvector of
40
Easy Problem
z.5
41
Hard Problem

• z.081

42
Realistic problem
43
Dispersion
44
Calculating derivatives
45
(No Transcript)
46
Dispersion
When computing A , must also construct A and
A, but these are diagonal. For optimizer
really need the derivatives of the
dispersion with respect to the design parameters-
another layer of differentiation
47
Dispersion gradient
(17)
k indicates which of the parameters and there can
be about 20
48
(No Transcript)
49
Experimental results
5 positive eigenvalues, 15 design parameters of
widths of layers and concentrations, 600 grid
points