Implementing the retraction algorithm for Symmetric Indefinite Band matrices

1
Implementing the retraction algorithm for
Symmetric Indefinite Band matrices

• Linda Kaufman, and students Daniel Clark, Serafim
  Stamatis, Edwin Torres, Ory Diaz, Steve Bang,
  Philip Kang, Diamond Flores, Mark Fila, Philip
  Nelson and Yomi Idris
• William Paterson Universitysponsored by NSF
  grant

2
Story starts with Gene in 1973

• I am a graduate student and Gene is my advisor
• He asks me to referee a paper by Jim Bunch on
  solving symmetric systems where the matrix is
  indefinite
• But analysis does not fit algorithm
• Bunchs algorithm looks at 1 column when making
  decisions
• Analysis only works if it looks at 2 columns

3
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 magnitudes of
  a11 versus other elements

4
Why 2 by 2 and pivoting
5
2 x 2 vs 1 x 1
6
fast forward 25 years to designing Dispersion
Compensation fibers

• Signal degradation and Restoration

7
Fiber Optics Design and Solving symmetric Banded
systems

• Linda Kaufman
• William Paterson University

8
(No Transcript)
9
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
10
Typical refractive profile and parameterization
11
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
12
Marcuse model for fiber wrapped around spool

• For a radially symmetric fiber

if r is the radius from the center of the fiber,
R is the radius of the spool we are lead to a
matrix problem of the structure.
bs and cs are 0(r/R)
13
Bandwidth spread with Bunch-Kaufman on banded
matrix
14
Nonlinear Rayleigh quotient for nonlinear
eigenvalue problem

• Try to minimize

Solution is at
15
Safeguarded Nonlinear Rayleigh quotient
16
Calculating dispersion from eigenvalues
17
Calculating Dispersion from Eigenvalues
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
18
Objective Create a algorithm to factor a
symmetric indefinite banded matrix A
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 1 -5 is indefinite-
determinant is negative -5 2 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
19
Uses

• (1)Solve a system of equations Axb
• If LALT D , D is block diagonal, where D
  contains 1 by 1 and 2 by 2 blocks
• set z L b
• Solve D y z for y
• Set x LT y
• (2)Find the inertia of a system- the number of
  positive and negative eigenvalues of a matrix
• If LALT D , The inertia of A is the inertia of
  D.
• A 2 by 2 block corresponds to 1 positive
  eigenvalue and 1 negative eigenvalue

20
Competition

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

21
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

22
Difficulties-1

• Consider the linear system
• .0001x y 1.0
• x y 2.0
• True solution is about x 1.001, y .999

23
Comparison with interchanging rows and columns
24
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

25
But with band pivoting for stability can ruin
bandwidth
Worst case r m, 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 b x x x x x x
x x c x x x x x x x
d x x x x x x
x x x x x x
26
Partition A as
Reset B B Y D-1 YT YZ 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 x
dr x x x x x as bs cs ds x x x x
x x x x x x x x x x as x
x x x x x bp x x bs x x x x
x x cp cq x cs x x x x x x dp dq
dr ds x x x x x x
x x x x x x
27
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 x 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 bt x x x x x x
cq x ct x x x x x x dq dr dt
x x x x x x x x
x x x x
28
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 u m x x x
x x x x x x x x
x Now eliminate u to restore bandwidth
29
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.
30
Implicit 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

31
Storage configuration
Mathematics
In Fortran
In C
32
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
33
Alternative -3-data structure vs. Blas
2 rank 1 corrections that overlap- best strategy
depends on size of overlap
But no unsymmetric rank 1 correction that just
affects triangle. Could treat one column at a
time.
R
Symmetric- rank 1
Could use symmetric rank 1 updates for 2
triangles, general for region R. Does one break
up region P?
P
Rank 2
34
Using less expensive transformations for
retraction with elementary transformations
Assume we have 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 x dr x x x x x
as bs cs ds x x x x x x x x x
x x x x x as x x x x x x bp x
x bs x x x x x x cp cq x cs x
x x x x x dp dq dr ds x x x x
x x x x x x x x
If s is the larger than p, q, and r so no
pivoting is necessary, we can get rid of whole
unwanted lower triangle with one level 2
transformation for the columns
For row transformations in the instead of
row transformations with non unit increment one
can process a column at a time with unit increment
If unwanted stuff has k rows, only pivot about
log(k) times
35
Can we also replace row transformations with
column transformations when using orthogonal
transformations?
Gene and posse (Michael Saunders, Phil Gill, and
Walter Murray) to the rescue
In 1974 they showed in rank 1 updating paper that
the product of Givens transformations have
special structure which for us comes down to
Recursive definition to get ßs and ps but
divides by cs ugh!
About 2/3 time of drots
36
Properties

• Space- (2m1)n-in order to save transformations
  compared to (3m1)n for unsymmetric G.E.
• Operation count depends on types of
  transformations
• Elementary between m2n/2 and 5 m2n/4
• Compared with between m2n and 2m2n for G.E.

37
Problem 1
Sample Problems
Problems 2 and 3
Problem 4
38
Comparison using Atlas Blas on Unix
39
Comparison with Lapacks unsymmetric band codes,
n3000 on random matrices on 0,1
40
Data Locality on positive definite problems
Increments of one Less paging Block algorithms
For snap4 define a_ref(a_1,a_2) A(a_2)a_dim1
a_1 For snap4invdefine a_ref(a_1,a_2) A(a_2
a_1)a_dim1 disp-(a_1) Where
disp2subdiag_beg1
41
Random matrices and unusual results
In snap4- storage is columnwise in C In snap4in-
storage is row wise in C In retraction using
orthogonal transformations seems to be as fast as
using elementary transformations
42
Elementary vs. Orthogonal on matrices AsI, A
random on (-1,1), n2000, m200
Apparently orthogonal transformations are not
that much more expensive
43
Elementary vs. Orthogonal on matrices AsI, A
random on (-1,1), n2000, m200
Apparently orthogonal transformations do not
always buy stability
44
Pivoting for 1x1 or not
Only pivot for 2 x 2 s
time eig n 2 m
growth error 1
10010001000 0 100 0 1
1.1e-15 3.6e-15 2 170 5021000 154 100
66.6 1e01 2.8e-14 1.6e-13 3 380
5001000 500 100 89.7 1 1.6e-13
3.2e-13 4 140 3831000 92 100 65
2e02 1.5e-14 1.0e-11 5 260 4981000 393
100 57.8 6e01 7.0e-14 3.4e-12 6 150
5001000 199 100 24.5 5e02 5.4e-14
2.1e-11
Pivot also for 1 x 1
time n 2 1 growth
error 1 11010001000 0 0 100
0 1 1.1e-15 3.6e-15 2 190
5021000 167 79 100 44.11e01 7.7e-14
3.7e-13 3 390 5001000 500 0 100
89.7 1 1.6e-13 3.2e-13 4 210
3831000 52 203 100 69.14e01 1.7e-13
8.3e-11 5 310 4991000 336 134 100
63.36e01 6.6e-14 3.2e-12 6 190
5001000 77 394 100 25.51e02 4.2e-14
2.1e-12
45
Problem 3
46
Problem 4
47
Problem 5
48
Future

• Adaptations for small bandwidth as in optical
  fiber code and multiple right hand sides(Mark
  Fila).
• Use special structure for orthogonal(Phil Nelson)
• Condition estimation (Diamond Flores, math)
• Creating examples in Fiber Optics (Steve Bang and
  Orielsy Diaz)
• Block approach (me)

49
Reference

• The retraction algorithm for factoring symmetric
  banded matrices, J. of Numerical Linear Algebra
  and Applications, 2007
• Calculating Dispersion in Optical Fiber
  Design, IEEE J. on Lightwave Technology, March
  2007
• Eigenvalue Problems in Optical Fiber Design ,
  SIAM J. on Matrix Analysis and Applications, Jan.
  2006.