Nonlinear Eigenvalue ProblemsbyRoss Adams LippertSubmitted to the Department of Mathematicsin partial fulfillment of the requirements for the degree ofDoctor of Philosophyat theMASSACHUSETTS INSTITUTE OF TECHNOLOGYJune 1998@Ross A. Lippert, 1998.The author hereby grants to MIT permission to reproduce and distributepublicly paper and electronic copies of this thesis document in whole or inpart.Author.... ...............Department of MathematicsMay 1, 1998Certified by..Accepted by..Ch "rman,Alan EdelmanAssistant ProfessorThesis SupervisorHung ChengApplied a ematics CommitteeAccepted by..Chairman, Departmental CommitteeRichard Melroseon Graduate StudentsJUN 01198LBRF IESNonlinear Eigenvalue ProblemsbyRoss Adams LippertSubmitted to the Department of Mathematicson May 1, 1998, in partial fulfillment of therequirements for the degree ofDoctor of PhilosophyAbstractThis thesis begins by exploring the problem left open by Wilkinson of computing the distanceof a given diagonalizable matrix to the nearest matrix with degenerate structure. Algorithmsfor finding the nearest nondiagonalizable matrix to a given matrix are presented and ana-lyzed, and an algorithm based on the work of Alexander Malyshev is presented for findingthe nearest nondiagonalizable matrix with a specific nontrivial canonical form.This thesis will then recast the Wilkinson problem in terms on a more general nonlineareigenvalue problem. Software written by the author to solve general nonlinear eigenvalueproblems is presented, as well as a variety of examples of problems which the software cansolve, including Local Density Approximation (LDA) problems arising in solid state physics.Lastly, this thesis presents several algorithms for computing operators commonly arisingin LDA problems (Laplacian, pointwise multiplication, etc.) using a truncated basis ofinterpolating scaling functions (the Deslauriers-Dubuc functions). The algorithms suggest ameans of improving the performance of common linear operations in any truncated compactlysupported wavelet basis.Thesis Supervisor: Alan EdelmanTitle: Assistant ProfessorAcknowledgmentsI would like to thank both Alan Edelman and Tomas Arias. Both of whom argued againstme at times when I thought I should stop because I could not derive some result or other.Needless to say, they usually ended up being right. I think that this is one of the mostvaluable lessons about research that I have learned.I would also like to thank Alan for discussing all aspects of this thesis with me. Whilepushing and prodding me to stay focused, he nonetheless gave me room to indulge my owncuriosities and form my own style. I would also like to thank my academic siblings, my fellowadvisees, Yanyuan Ma, Parry Husbands, and Peter McCorquodale, for all their help in mygraduate existence.There are many people who have stimulated my thoughts and provided me with a greatdeal of material with which to work. I'd like to thank Gil Strang for his wavelets class, SteveSmith for his valuable work on Stiefel-Grassmann optimization, Alexander Malyshev for hisinspirations on general multiple eigenvalue problems, and Erik Elmroth for his insights intothe double eigenvalue problem.I would like to thank the many friends I have made at MIT whose company provided awelcome distraction from academic work, notably Elizabeth Shapere and Nancy Nangeroni,as well as sundry others who have shared drinks with me at the Phoenix Landing.I finally would like to thank my parents for their encouragement and for providing mewith the means to pursue my own interests wherever those interests took me.Contents1 Foreword 172 The Computation and Sensitivity of Double Eigenvalues 192.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2 The geometry of eigenvalues ................... ........ 222.2.1 The level sets of eig ........ ......... ....... 222.2.2 The double eigenvalue variety ................... ... 252.2.3 An illustration by characteristic polynomials . ............. 262.2.4 A pseudospectra "movie" ................... ...... 292.3 Finding the nearest double eigenvalue ................... ... 332.3.1 A perturbative look at double eigenvalue averaging . ......... 332.3.2 An iterative fixed-point algorithm . .................. . 362.3.3 The singularities of the spectral portrait . ............... 372.4 Nearest A's are basically J2's . . . . . . . . . . . . . . . . . . . . . . . . . . 412.5 Conditioning and stability ............................ 462.5.1 Sensitivity analysis of A = A + auvH . ................ . 472.5.2 Ill-posed problems and principal curvatures . .............. 532.5.3 Stability of the fixed point algorithm . ................. 592.5.4 Sensitivity of the critical points of the spectral portrait ........ 602.6 Nearest A's come from simple amin . ................ . . . . . . 622.7 Triple eigenvalues .............. .. . .. ... .. .... 663 The Computation of Nearest k-fold Eigenvalues 693.13.23.33.43.53.63.7Introduction . ..... ............Malyshev's formula for dist(A, J2) . . . . . .Formula for dist(A, Jk) . . ...........Finding A....... ............3.4.1 Generic success: tr(UHV) 0 . . . .3.4.2 Nongeneric success: tr(UHV) = 0 . .The 2-norm "normals" . . .... .......A possible rationale for the 2-norm normalsConcluding remarks ............... . . . . . . . . . 69.. . . . 70. . . . . . . . 70... .. . . 73. . . . . . 76. . . . . . . . 77. . . . . . . . . . 78. . . . . . . . . . 82. . . . . . . . 8485... . .. .... . 85. . . . . . . 86. . . . . . . . . . . 90. . . . . . . . 91. . . . . . . 92. . . . . . . . 93. . . . . . . . . . . . . 93. . . . . . . . . . . . 94. . .. .. . 95Solving General Nonlinear Eigenvalue Problems4.1 Introduction ........ .............4.2 Calling sgmin . . . ....................4.3 Sample problems and their differentials . . . . . . .4.3.1 The Procrustes problem . ..........4.3.2 Nearest Jordan structure . . . . . . . . ...4.3.3 Trace minimization ..............4.3.4 Trace minimization with a nonlinear term4.3.5 Simultaneous Schur decomposition problem4.3.6 INDSCAL ...... . ......4.4 Running our examples . ...... .........4.4.1 A sample Procrustes problem . . . . . . . .4.4.2 A sample Jordan block problem . . . . . . .4.4.3 A sample trace minimization problem . . . .4.4.4 A sample LDA toy problem (trace minimization with nonlinear4.4.5 A sample simultaneous Schur decomposition problem . . . . .4.4.6 A sample INDSCAL problem (simultaneous diagonalization) .4.5 Making modifications (the structure of the code) . . . . . . . . .4.5.14.5.2. . . . 96. . . . 97. . . . 98. . . . 99term) 101. . . . 102. . . . 103... . 104Subroutine dependencies . . . . . . . . . . . . . . . . . . . . . .Under the hood ..............................1051074. . . .. . . . .. . . ... . . .. . . ... . . .. . . .. . . . .4.5.3 What to modify .. . . . . . . . . . . . . .4.6 Geometric technicalities . . . . . . . . . . . . .109. . . . . . . . . 111Manifolds .. ..... .......... .......The difference between a tangent and a differential . .Inner products, gradients, and differentials . . . . . . .Getting around Stief(n, k) ................Covariant differentiation . . . . . . . . . . . . . . . . .Inverting the covariant Hessian (technical considerations1111131141151161185 Multiscale Computation with Interpolating Wavelets5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . .5.2 Introduction to interpolets ...... .........5.2.1 Interpolet Multiresolution Analysis . . . . . . .5.2.2 Interpolet transforms . . . . . . . . . . . .5.3 Accuracy of interpolet approximation . . . . . . . . . .5.4 Truncated bases ............. .........5.5 Multilevel algorithms for V 2 and other operators . . .5.5.1 V 2 in 1-level decomposition . . . . . . . . . . .5.5.2 Computing other operators ...........5.5.3 V 2 in 2-level decomposition . . . . . . . . . . .5.6 Efficient implementation .................5.6.1 Data structures . . .................5.6.2 Implementation overview . . . . . . . . . . . ..5.6.3 Block-to-block subroutines . . . . . . . . . . ..5.7 Results and conclusions ........ . .......A Matlab source for Stiefel-Grassmann optimization121.. ... ...... 121.. ... .... .. 126. . . . . . . . . . . 126. . . . . . . . . 131. . . . . . . . . . . 133.. ... ... ... 134. . . . . . . . . . 142. . . . . . . . . . . 144.. ... ... ... 149. . . . . . . . . . . 151.. ... ... ... 154. . . . . . . 155. . . . . . . 155. . . . . . . . . 156.. ... ... ... 1561674.6.14.6.24.6.34.6.44.6.54.6.610List of Figures2-1 I]etA 12 (solid) and I etAl 2 (dashed) vs. t for A defined by (2.1) and its nearestdefective matrix ............. ... ................ 212-2 (Left) Mx defined by the constraint eig(A) - A = 0 for fixed A. The gradientof the constraint function points in the direction of increasing A and hasmagnitude inversely proportional to the distance to MA+6A. (Right) the MAdefined by the constraint amin(A - AI) = 0 for fixed A, giving a gradient ofunit length in the direction of increasing A, the distance to next MA beingproportional to s = vHu............. ...... .......... .. 242-3 The envelope can be identified by the critical points of the distance from Ato M Ax. .. ......... . . . .. . .. ... . . .. . .. . . . . ...... .. 262-4 A graph in the a-b plane of the lines of polynomials of the form p(x) =3 + ax + b satisfying p(A) = 0. .................. ....... 272-5 A graph in a-b-c space of the double root variety, or swallowtail. ....... . 282-6 At a = 0 we have a single eigenvalue at -0.3 and a parabolic valley for thedouble eigenvalue at 0.15. The conservation law indicates that there is also asaddle in the portrait (it sits between the valleys). . .............. 302-7 As we begin to move off D, the double eigenvalue splits into two complexconjugates which move off from their origin very quickly. The valley of thedouble eigenvalue (charge +1) has split into two valleys (+1 each) with asaddle (-1) between them. . .................. . ...... 302-8 We then see a critical change to the left of A = 0.15 as the saddle whichinitially sat between the two valleys at a = 0 decays into two saddles and amaximum. .................. ................. 312-9 The maximum which was spun off from the decay collides with the saddleat A = 0.15. At the focal distance a = 0.45 the maximum and saddle haveannihilated each other, creating a degenerate shoulder. . .......... . 312-10 After the collision, the critical points have exchanged charges. The A = 0.15critical point is now a maximum while a saddle heads from it to the right. 322-11 The fixed point algorithm converging. .......... ....... . . . . . 382-12 The fixed point algorithm approaching a limit-cycle . ............. . 392-13 A possible problem finding a normal from A to A illustrated in characteristicpolynomial space. ................... ............. .. 422-14 The normals of the cubic double root variety forming their parabola-like enve-lope. The concave down arc in the plot is another part of the envelope whichoccurs over cubics with complex coefficients. . ................ . 472-15 The normals of the swallowtail. A faint hood-like shape begins to emerge... 482-16 The swallowtail and its hood-like evolute of ill-conditioned problems .... 482-17 In the cubic polynomial space, a particle moves along the double root variety,with velocity v and centripetal acceleration a(v, v). . .............. 542-18 The centripetal acceleration of the particle moving along the double root va-riety can be associated with a circle centered at the focal point of the doubleroot variety. ................. ... ............... 552-19 Below the focal point of a hemispherical surface, the distance minimizing pon the hemisphere to a given q is at the bottom of the surface. . ...... . 642-20 Above the focal point of the hemispherical surface, the bottom point of thehemisphere is now a maximum of distance (the nearest p on the hemisphereto q has become either of the end points of the hemisphere). . ....... . 644-1 Procrustes Problem ................... ............ 984-2 Jordan Problem ................... ........... .. 994-3 Trace Minimization Problem ........... ...... ......... 1014-4 LDA Toy Problem .............. ................. 1024-5 Simultaneous Schur Problem ........... ............... 1034-6 (INDSCAL) Simultaneous Diagonalization Problem . ............. 1054-7 Our code tree: dependencies of various modules . ............. . 1064-8 The unconstrained differential of F(Y) can be projected to the tangent spaceto obtain the covariant gradient, G, of F. . .................. . 1154-9 In a flat space, comparing vectors at nearby points is not problematic..... . 1174-10 In a curved manifold, comparing vectors at nearby points can result in vectorswhich do not lie in the tangent space. . .................. ... 1175-1 (a) shows 0(x/2) = -'(x - 3) + -2(x - 1) + O(x) + -Q(x + 1) + -'1(x + 3).(b) shows the functions =1q(x - 3), 9-(x - 1), q(x), 9¢(x + 1), - ¢(x + 3). 1275-2 (a) the smooth function f(x) = ( + 3)3e-( )4 . (b) the approximation to f(x)in 1 0(Z). (c) component of the approximation in 11 (C1). (d) component ofthe approximation in 1(D 1). .......................... 1355-3 A visual summary of the 1-level touching condition. Solid plots representfunctions centered at points in S. Dotted plots represent functions centeredat points in Z n - S. Tick marks delimit the overlap of the two functions. 1375-4 Our example of the 1-level touching good basis in one dimension. Note thatthe two functions plotted do not touch. Tick marks denote the set S C Z. 1385-5 A generic example of a truncation which meets our definitions of good and1-level touching. In black are the points of S C Z. . .............. 1395-6 An example of a 2-level touching good basis which can be used for a diatomicmolecule (two atomic cores). The points in S C Z are where residual valuescould be significant ............................... 1405-7 The previously used implementation is on the left, and the implementationemploying a 1-level touching algorithm is on the right. (Note the differencein scale on the vertical axes.) ................... ....... 1585-8 The previously used implementation is on the left, and the implementationemploying a 2-level touching algorithm is on the right. (Note the differencein scale on the vertical axes.) . ....... ........... ....... 1595-9 The previously used implementation is on the left, and the implementationemploying a multilevel algorithm is on the right. (Note the difference in scaleon the vertical axes.) . . . ....... ....... ......... . 1605-10 The condition number of the Laplacian operator represented in truncated mul-tiresolution interpolet basis as a function of the number of refinement levels kwith and without simple diagonal preconditioning and compared with the con-dition number in an orthogonal basis with the same resolution. Lines indicateresults for bases with a single atomic center of refinement, and points representresults for two atomic centers corresponding to the nitrogen molecule. . . 1625-11 Convergence of the solution to Poisson's equation for the nuclear potentialin a nitrogen molecule in an interpolet basis with k = 8 levels of refinementabout each nucleus. ................... ............ 1635-12 Effective condition number of Poisson's equation for the nuclear potential ina nitrogen molecule with simple diagonal preconditioning as a function of k,the number of levels of refinement about each nucleus. . ......... . . 163List of Tables2.1 phenomena observed along the normal .................. .. . 324.1 A short list of the optional arguments for sgmin . ............... 885.1 Notation for multiscale computations. . .................. ... 125^II.~OI~I~C'SC" L~-- ir~lllll*IIT--~ L-X-I_-LII.1-~--~- 1~1 _-I I------II~---)YI----- s 11 1111 --- -16Chapter 1ForewordWhat is a nonlinear eigenvalue problem? For us, it has become a catch-all phrase for a flavorof eigenvalue problem that one does not currently expect to see in an LAPACK style softwarelibrary or a Matrix Computations style text. One way to "nonlinearize" the usual Ax = Axproblem is to have the elements of A depend on x and A. Another example arises when onewishes to regularize the eigenstructure of A according to some prescribed canonical form(i.e., given A what is the nearest A such that A has some prescribed eigenstructure). Char-acteristic of such nonlinear eigenvalue problems is a variational or optimizational procedureat the forefront of the computation. Polynomial eigenproblems are sometimes referred to asnonlinear (e.g. Demmel [23] uses this nomenclature, but we are not in favor of this use ofthe term).Some nonlinear eigenvalue problems have structure from which specialized algorithmsfor solution can be found. There is a wealth of geometric theory which can lead one to thenearest defective A to a given A. This geometric theory is elaborated in Chapter 2. Therealso appears to be some geometric structure to finding the nearest A to A with an arbitrarycanonical structure, though this is a much more difficult problem than defectiveness, andthe theory is still incomplete (Chapter 3).More general nonlinear eigenvalue problems may not have as rich a geometric structure aswas found in Chapters 2 or 3 but can still be formulated as optimizations of some functionF(Y) over the set of k-dimensional subspaces of ER or Cn, much the same way that theeigenvalue problem can be cast as an optimization of the Rayleigh quotient:min xHAx.This leads one to consider optimizations over the Stiefel or Grassmann manifolds (or a setof manifolds in between which we call Stiefel-Grassmann manifolds a.k.a. flag manifolds).Fortunately, these manifolds have enough structure that their differential geometry is com-putable in terms of elementary functions and matrix exponentiation. Chapter 4 providesthe computational details for abstracting the optimization algorithms of flat spaces to theircurved analogues on Stiefel-Grassmann manifolds.One particular Stiefel-Grassmann optimization (actually, it is a pure Grassmann opti-mization) is the computation of electronic ground states in the Local Density Approximation(LDA). In an LDA problem the eigenvalue problem isA(Y)Y = YA,where Y is n x k and A is diagonal. This is an eigenvalue problem where the operator Adepends on the eigenspace Y, of the sort mentioned earlier.Work on this problem led to development of the practical aspects of wavelet computations,as the Y were being represented in a wavelet basis of interpolating scaling functions. Thisled to a digression from general nonlinear eigenvalue problems to the improvement of theperformance of common wavelet operators. Thus Chapter 5 may appear to be somewhat outof place in this thesis, but I hope that this explanation suffices. Chapter 5 appears to havegeneralizations to wavelet bases beyond those of interpolating scaling functions. Whetherthis is possible or not remains to be seen.Chapter 2The Computation and Sensitivity ofDouble Eigenvalues2.1 IntroductionA fundamental problem in numerical analysis is to understand how eigenstructure affects thenumerical behavior of algorithms. Classical insights come from studying the conditioning ofthe eigenvalue problem. From this point of view, ill-conditioning is related to nearness to adefective matrix. More recent insights are based on pseudospectra where interesting behavioris associated with non-normality. In this chapter, we revisit the classical insights with theadded technology of differential geometry. Pseudospectra are windows into n2 dimensionalspace. Singularities in the spectral portrait are relics of a defective matrix. "Singularities"of these singularities (to be defined later) indicate a center of curvature of the defectivematrices and play a role in conditioning.If A is an n x n, diagonalizable matrix, we can find a neighborhood about A in matrixspace where the eigenvalues are all analytic functions of the elements of A. Thus we haven analytic functions A = eigj(A). An eigenvalue Aj is said to be ill-conditioned if firstorder variations A result in very large first order variations in Aj. A textbook example is the12 x 12 Frank matrix12 11 10 ... 3 2 111 11 10 ... 3 2 10 10 10 ... 3 2 1F 12 = 0 0 9 ... 3 2 10 0 0 ... 2 2 10 0 0 ... 0 11An c perturbation of the upper right element results in a change in the smallest two eigen-values of order 10'6 (a relative change of order 2 x 108).Wilkinson [67, 68] showed that if the conditioning of some of the eigenvalues of a matrixis poor, then that matrix must lie near some defective matrix. He further presented severalbounds on this distance. In particular, the Frobenius norm distance from F12 to the nearestmatrix with a double eigenvalue is only 1.822 x 10-10.Nearness to a defective matrix can also sometimes explain transient behaviors in thematrix exponenital. Such transient behavior is explored by Trefethen [63], who exhibitedthe transient behavior ofA = , (2.1)0 -2a matrix with distinct negative eigenvalues. His message from this example is that defectiveeigenvalues are neither necessary nor sufficient for transient growth. It is curious, however,that in his example, there is a defective matrix lurking. Let U be defined by the singularvalue decomposition UEVT = A + I. ThenUTAU = -1.5 5.04950.049509 -1.5and thus A is a distance 0.049509 (= 9x 10-31 IAl F) from a matrix A with a double eigenvalue,-1.5. A plot of the norms of the exponentiations of A and A in Figure 2-1 shows that thedefective matrix does play a role.1.20.8 -0.6 -0.4 -0.200 0.5 1 1.5 2 2.5 3 3.5 4Figure 2-1: IletA 12 (solid) and IletA 2 (dashed) vs. t for A defined by (2.1) and its nearestdefective matrixGiven the significance of a lurking nearby defective matrix in these examples, this chapterstudies the computation of the shortest (least squares) distance of any matrix, A, with dis-tinct eigenvalues, to the nearest matrix, A, with at least one repeated eigenvalue. Wilkinsonreviewed and critiqued several bounds on this distance, although no procedure was given forthis computation. Recently, Malyshev [49] has studied the 2-norm version of this problem.There has also been recent work on the computation of the distance to more general Jordanstructures in [22, 33, 29, 30], and also in our sg_min template (see Chapter 4).This chapter closely studies the n2 - 1 dimensional surface of defective matrices andpresents two methods (the sg_min approach is a third, though it will not be discussed inthis chapter) which will identify the A for a given A and give appropriate condition numberson A and its repeated eigenvalue. One method proposed will relate the nearby A to thecritical points of the spectral portrait. We shall find that the conditioning of A is directlyrelated to the degeneracy of the associated critical point. A second method proposed relies onsuccessive improvements to an approximation of the repeated eigenvalue. This method willhave a stability which is also dependent upon conditioning of the critical points of the spectralportrait. Although we feel that there is more theoretical value in the first method, we presentthe second because of its easy implementation in terms of common Matlab operations.Throughout this chapter, unless otherwise stated, we will use the Frobenius matrix normand inner product. We also assume passing familiarity with concepts from differential ge-ometry and singularity theory. (We recommend [13] as an excellent reference for singularitytheory.)2.2 The geometry of eigenvaluesIn this section, we will examine the level sets of the eig function and their envelope, thedouble eigenvalue variety. We will consider a matrix A with eigenvalue A to belong to asurface MA in matrix space. This is a surface of co-dimension 1 with a continuous normalvector field almost everywhere. The envelope of the MA, D, the double eigenvalue varietyis comprised of all matrices with at least one repeated eigenvalue. This surface is also ofco-dimension 1 and also has a continuous normal field almost everywhere.2.2.1 The level sets of eigConsider the function eig(A) returning some eigenvalue of A. If A is diagonal, then eig maybe presumed analytic in a neighborhood of A. We define MA = eig-1 (A) to be the set of allmatrices over either R or C with eigenvalue A. MA is a variety which is differentiable almosteverywhere. MA may be thought of as a level set of eig . All A E MA can be characterized bythe constraint equation e(A, A) = 0 for A fixed, where e(A, A) = eig(A) -A in a neighborhoodof A.It is well known (see [39], pp. 320-324) that if A = eig(A) is a simple eigenvalue withright eigenvector x (Ax = Ax) and left eigenvector y (yHA = AyH), then variations in A relateto the variations in A by the formulay HAzy= . (2.2)This formula capt res t e differ ntial b havior of A(t) eig(A(t)) at diagonalizable matrices.This formula captures the differential behavior of A(t) = eig(A(t)) at diagonalizable matrices.More generally, if A(t) is independent of A(t) the derivative of the constraint function, e, isd yHAxe(A(t), A(t))= -A. (2.3)dt yHg:(Note, the condition A(t) = eig(A(t)) implies e = e 0 reducing to Equation (2.2).)The differential behavior of any scalar (analytic) function z = g(w) can be captured bya gradient and a Riemannian (Hermitian) inner product.S= (Vg, b) .Therefore, using the Frobenius inner product on matrices (A, B) = tr(AHB), we may rewritethe perturbation formula (2.2) asyZH(YXA)SHy'from which we see that the gradient of eig is Veig = ~ , and the differential behavior of eise(A(t), A(t)) = (Veig , A)- A.Geometrically, the gradient of any function is perpendicular to that function's level sets.Thus, 1 is normal to the surface MA. If we select y and x such that lyll = IIjxj = 1 thenxHyN = yxH is the unit normal (unique up to sign [or phase]). We may identify s = xHy as theinverse of the length of the gradient (also up to sign [or phase]).There is an unfortunate characteristic of the constraint function e. At points of M\ wherethe eigenvalue, A, is ill-conditioned, the right and left eigenvectors of A are close to beingperpendicular, thus s = xHy is very small and the magnitude of Veig approaches infinity. IfA were a matrix with a nontrivial Jordan block for A then xHy would vanish and the gradientof eig (as well as the differential of e) would be infinite.To obtain a better differential behavior at A with ill-conditioned A, we must dispensewith e and use a constraint equation with a better differential behavior to characterize theMA surfaces. If we let f(A, A) = Umin(A - AI), we may define MA by the constraint equationV eig = N/sFigure 2-2: (Left) M defined by the constraint eig(A) - A = 0 for fixed A. The gradient ofthe constraint function points in the direction of increasing A and has magnitude inverselyproportional to the distance to M+6A. (Right) the MA defined by the constraint amin(A -AI) = 0 for fixed A, giving a gradient of unit length in the direction of increasing A, thedistance to next M\ being proportional to s = vHu.f(A, A) = 0, A fixed. The function f has the differential behaviorf (A(t), A(t)) = uHAv - AuHv, (2.4)where u and v are the right and left singular vectors of A - AI (when f(A, A) = 0 the phasesof u and v are such that the right hand side of (2.4) is real and nonnegative).Observe that u and v are unit length left and right eigenvectors (respectively) of A, andthat setting df = 0 givesuHV = (UVH, A), (2.5)which recovers (2.2). The gradient f with A fixed is the unit normal to MA, while the cosineof the eigenvectors weights the A term. Figure 2-2 summarizes the distinction between thedifferential behaviors of e and f.Another way to view f(A, A) is as the distance from A to M\. Since f is the distancefunction to M, the magnitude of its gradient (with A fixed) must be of unit length andperpendicular to MA. (One can derive 2.5 directly from this interpretation without differen-tiating f.)2.2.2 The double eigenvalue varietyThe double eigenvalue variety, D, is the envelope of the family MA. Recall that the envelopeof a family of surfaces St parameterized by t can be defined as a tangent surface, E, havingthe property that for x E E, there exists a to such that x E Sto and the tangent space of Eat x coincides with the tangent space of Sto at x.The characterization of the MA in terms of f(A, A) = ,min(A - AI) allows us to analyze2D. It is well known ([66], pp. 9-10) that A c D iff there exist A, u, and v such that Av = Av,UHA = AUH, and UHv = 0. Thus, for a differentiable curve A(t) C D with repeated eigenvalueA(t), we have by (2.4)f(A(t), A(t)) = d f(A(t), A(t)) = 0 = uHAv = (uvH A),giving the unit normal N = uvH of D9 at A(t), confirming that 2D is the envelope of the MX.We can also recover an important property of envelopes since f(A, A) is the distancefunction from A to MX. It is the case for all differentiably parameterized families of surfaces,that the envelope is given by the the critical points of the distance function, i.e. A E D~ whenf(A(t), A) = 0 for some A where A is the nearest element of Mx to A (see Figure 2-3).In some sense, we have arrived at a geometric reasoning for the choice of Wilkinson's "lo-cally fastest" perturbation to a matrix A near D which would move A to D nearly optimally.If one suspects that two eigenvalues, A1, A2, of A can be coalesced with a small perturbationE, Wilkinson suggested using a perturbation of the form E = a(y1 + y2) (xl + x 2)H where theyi and the xi are the left and right eigenvectors of Ai and a is adjusted so that A + E E D(Section 9 of [68]).From the preceding results, we see that if A E D then the normal line A + aruvH (pa-rameterized by a) passes through D along a perpendicular, thus the perturbation E = auvHcan be thought of as the perturbation which maximizes the change in the distance to D forsmall a and is then the "locally fastest" perturbation to a nearby diagonalizable matrix A.Figure 2-3: The envelope can be identified by the critical points of the distance from A toMA.If one is close to D then one expects (x 1 + x2)/2 will be close to v and (yi + y2)/ 2 willbe close to u. The "locally fastest" perturbation is then expected to be close to the actualnormal of the variety at the actual minimizing point. We will elaborate on this heuristic inSection 2.3.1.2.2.3 An illustration by characteristic polynomialsFor purposes of "picturing" D and MA geometry, we consider monic nth degree polynomials,which may be thought of as projections of the n x n matrices via the characteristic polynomial.Consider the cubic polynomials of the form p(x) = x 3 +ax +b (corresponding to the traceless3 x 3 matrices). Polynomials with a fixed root A can be represented by a linear relationbetween a and b of the formA3 + aA + b = 0.A collection of these fixed root lines is shown in Figure 2-4.The envelope of those lines is a cusp. By varying the constraint equation3A2A, + aA + tA + b = 0,-0 25 -0.2 -0 15 -0 1 -0.05 0 0.05 0.1 015 0.2 025bFigure 2-4: A graph in the a-b plane of the lines of polynomials of the form p(x) = x3 +ax+bsatisfying p(A) = 0.we obtain the gradient of A as a function of a and b,-(,\, 1) - (a, b)3A2 - aWe may define a unit normal N = ,1) and set the inverse magnitude of the gradient3A2 +aS 2 -+1The envelope of the lines of this single root family of polynomials is given by s = 0, ora = -3A 2 (and thus b = 2A3 by the constraint). This gives a curve of polynomials withdouble root A, seen as the cusp-shaped envelope in Figure 2-4. Observe that the cusp pointis the polynomial p(x) = x3 , the triple root polynomial.More complicated shapes arise for quartic polynomials of the form x 4 + ax 2 + bx + c. Inthis case, the double root variety is a collection of multiple cusps, the swallowtail.-1- " 2. .-2 1-3, " -13 2 -210-2 -3-3Figure 2-5: A graph in a-b-c space of the double root variety, or swallowtail.In matrix space, the cusp may be formed from the 1-dimensional family of matricesv1 - 2A 2 A-2A 1 - 2A20 A0 0 0)One then has det(xl - M(A)) - x3 - 3A 2x + 2A3, and it is clear that uvH= 0 0 01 0 0)One can verify that det(xl - M(A) - auvH) - det(xl - M(A)) = c (Ax + 1), parallel to thenormals of the double root polynomials. For the quartic case, we may use the two parametersurfaceA0001 -A +A-A - fl 1 + A4 -_ A2(1 _)0 -A + ,/A2-A 3 +,1A)for which det(zI - M(A, q)) = (Z 2 - A2) 2 -t(x - A)2./AM(A)= 00These matrix families embed the cusp and the swallowtail in matrix space as subsurfacesin (respectively) 9 and 16 dimensional spaces.2.2.4 A pseudospectra "movie"Since we wish to find the nearest matrix A E D9 to a matrix A, it is useful to study the normalline in matrix space, A + auvH (parameterized by a). We proceed by taking a journey alonga normal to D and plotting the pseudospectra as we go. The spectral portrait will be our 2dimensional window onto n2 dimensional matrix space.We have seen in the previous sections that along this line, the double eigenvalue separatesinto two distinct eigenvalues which move apart initially with great speed. Additionally, theywill leave behind a stationary point whose location and character will remain invariant aswe go up the normal line for some finite length.We start out on the cubic polynomial example matrix(A /1 2A2 AM(A) = 0 -2A V1 - 2A20 0 Awith A = 0.15. The series of spectral portraits presented (Figures 2-6 through 2-10) are thoseof M - auvH for various interesting a.The reader studying these plots may wish to recall a certain "conservation law" whichcomes from planar bifurcation theory (see [59, 16]). The "law" is that if all the criticalpoints of a surface are nondegenerate, they may be assigned charges, all extrema get a+1 charge and all saddles get a -1 charge, which will be locally conserved through allcontinuous deformations of the surface. Additionally, spectral portraits are continuous, andtheir contours approach concentric circles about the origin for large IA . Thus any spectralportrait can be continuously deformed into a bowl and therefore has a total charge of +1.The contour lines in Figures 2-6 through 2-10 are not uniform, but have been chosen tobe dense about the current value of a and to highlight other important qualitative features.We can summarize the sights we have seen on this journey in Table 2.1. The reader-04 -02 0 02real04 06 08Figure 2-6: At oa = 0 we have a single eigenvalue at -0.3 and a parabolic valley for thedouble eigenvalue at 0.15. The conservation law indicates that there is also a saddle in theportrait (it sits between the valleys).-04 -02 0 02 04 06 08real-04 -02 0 02 04 06 08realFigure 2-7: As we begin to move off D, the double eigenvalue splits into two complexconjugates which move off from their origin very quickly. The valley of the double eigenvalue(charge +1) has split into two valleys (+1 each) with a saddle (-1) between them.~0I 1 1real-0.4 -0 2 0 0.2 04 0.6 08realFigure 2-8: We then see a critical change to the left of A = 0.15 as the saddle which initiallysat between the two valleys at a = 0 decays into two saddles and a maximum.0 6 060.4 0402- 0.20 0E E-02 -02-06 =0.30 -06 -0.-04 -0.2 0 02 04 06 08 -04 -02 0 02 04 06 08real realFigure 2-9: The maximum which was spun off from the decay collides with the saddle atA = 0.15. At the focal distance a = 0.45 the maximum and saddle have annihilated eachother, creating a degenerate shoulder.-02-04- _-06- (7=0.7-04 -02 0 02 04 06 08realFigure 2-10: After the collision, the critical points have exchanged charges. The A = 0.15critical point is now a maximum while a saddle heads from it to the right.a observationa = 0.0 double eigenvalue0.0 < a < 0.45 double eigenvalue splits leaving behind a saddlea = 0.45 the saddle becomes degeneratea > 0.45 the saddle has become a local maximumTable 2.1: phenomena observed along the normalmay wish to review this section as he/she reads subsequent sections to understand themathematical significance of the observations listed here and to illustrate the mathematicalphenomena described in subsequent sections.2.3 Finding the nearest double eigenvalueOne way to find the nearest A E D to a given A is by directly finding the repeated eigenvalueA of A. If the repeated eigenvalue A were known then one can take A to be the nearest matrixto A that has eigenvalue A (using the singular value decomposition to find A).This section explores methods of calculating that A. We will begin by exploring thecommon heuristic of averaging two nearby eigenvalues of A. A fixed point algorithm based onan improvement of this idea will be presented. We then show that A can also be determinedby an examination of pseudospectra of A, namely via the critical points of the spectralportrait.2.3.1 A perturbative look at double eigenvalue averagingA small perturbation in the normal direction of a generic point A E D~ results in a diagonal-izable matrix A = A + auvH. The J2 block that is part of the structure of A bifurcates intotwo J1 blocks. In this section, we will examine the heuristic of approximating the doubleeigenvalue of A by the average of the two nearest eigenvalues of A by examining the Puiseuxseries (see [43]) in a.Puiseux series are nothing new in the description of the bifurcation of multiple eigenval-ues, however, the additional contribution that we add is the derivation of the series alongthe geometric normal of D. We will also derive the second order behavior of the series alongthe normal which can be used to measure the error of the averaging heuristic.We quantify the non-surprising fact that the averaging heuristic tends to break sooner(i.e. for smaller a) when A is close to a nongeneric point of D, that is, is close to being moredefective than J2 or derogatory. Later sections will show that one cannot expect to find awell-conditioned double eigenvalue if the nearby A is close to being nongeneric, and thus theaveraging heuristic must suffer from difficulties which will also plague the more sophisticatedalgorithms suggested in this chapter.We adopt the convention, in this section and throughout this chapter (with the exceptionof Section 2.5.2), of using v and u for the right and left eigenvectors respectively, for A Dand of writing the singular value decomposition of A asA = AI + uOv 'Hwhere U = (ul u 2 ... Un-1), V = (v 1 v2 ... vn-1 ), and E = diag(al,..., a_).Additionally, we set r = UHv and 1 = VHu, noting that v = Ur and u = VI since uHv = 0and (U u) and (V v) are unitary matrices.We now consider a Puiseux series expansion of the eigenvalues and right eigenvectors ofA(a) = A + ruvH given by1 3A + A(a) = A + wo A1 + w2 A 2 + wcaA 3 + ---and1 3v v(oa) = v + WUV 1 + w2 U 2 + W 32v3 + ...satisfying(A (A + A (a)) I) (v + v(r)) = 0, (2.6)where w2 = 1 and vHvk = 0. Since vHvk = 0 we will write vk = Vrk for some n - 1 x 1 vectorrk and let r(a) = VHv(). The series for the average of the two eigenvalues is given byA + -((a)(,= 1 + A(a) ,=_1) = A + OA2 + O2 4 + .'.2Thus, averaging will be a good approximation to A as long as the A2 term is relatively small.Substituting A(a) - A = auvH+ UfVH into (2.6) we have(UEVH"+ UUvH A (o))(v + v()) = 0,which can be simplified to(U VH - A(u))v(c) - A(u)v + ou = 0. (2.7)We split (2.7) into two sets of equations, first by multiplying on the left by uHA(a)a = uH(), (2.8)and then by multiplying on the right by UH,tr(a) - A(a)&UHr(a) = A()r. (2.9)We may invert Equation (2.9) to obtainr(u) = A(u)( - ()UH)- 1 .To find A(a) we substitute into (2.8) to geta= (a))uHv(o) _) = A(c)lHr(r) = 2 (a) H(I - H()UH )-1 ror= Z uH(f-l UH)k0 k+ l(u) = Ck k+l(a), (2.10)k>1 k>1which must be inverted to obtain the series for A(a) (the reader may note that VC-iUH =(A - AI)t, the Moore-Penrose inverse of A - AI, and that cl = u(A - AI) t v = uHw where wis the right generalized eigenvector of A).We seek only the first two coefficients of A(a). They are given by1 = 1 (2.11)UH(V -1 )2V2A2 ()2 (2.12)__y__. i- ~n-I ar~arr~--~-----~- L-.~L^m-r~-^^~--ril-- y.-I Frraa~-------llYII ---~l~--*ll---r;iLLr . _ - .1~ ---- ~--111 ~-~ _11_11--~ X- I~S~~-. .~II~C~I*CC--l C- 11~-^--~----7~~-~~~1~-~-~~-From (2.11) and (2.12) we see that the second order term in the expansion becomes significantwhen(uHw)3(UH(VV-1U) 2 V) 2Thus we see that the averaging heuristic can be expected to work poorly if uHw (= lH-lr)is small or if I VIE-1UH 2 = I (A-A)t i2 = 1 is large. There is an algebraic interpretationfor both of these cases. If uHw - 0 then it can be shown that A must have at least J3(A)defectiveness. If I(A - AI)t 12 = oc, then A - AlI must have two vanishing singular values,and thus has a derogatory Jordan structure in A.In summary we see that the averaging heuristic will produce close approximations toA for nearby A so long as the nearby A are not close to any less generic structure than aJ2(A) block. For all of the algorithms presented in this chapter to calculate nearby doubleeigenvalues, analogous conditions will be found which govern their stability and/or accuracy.2.3.2 An iterative fixed-point algorithmEdelman and Elmroth [28] introduced an algorithm to find the nearest A E D to a givenmatrix A. The algorithm is the simultaneous solution of "dual" problems which can bestated as follows:Problem 1:Given a matrix A and an approximation to the perturbation bringing A to D), find therepeated eigenvalue of A.Problem 2:Given a matrix A and an approximation of the repeated eigenvalue of A, find the smallestperturbation bringing A to D.For Problem 1, they use a clustering algorithm on the eigenvalues of A - auvH, whereuuvHis the approximation to the normal bringing A to D (a a scalar and uvH an approximatenormal). For Problem 2, they approximate u and v with the left and right singular vectorsof -min(A - AI) for the given approximate repeated eigenvalue A.This resulted in an extremely compact Matlab routine, which reads roughly aslam = guess;while (1)[u,s,v] = svd(A-lam*I);Ahat = A - u(:,n)*s(n,n)*v(:,n)';e = eig(Ahat);lam = cluster_routine(e);endwhere the cluster routine takes the mean of two nearest elements of e.The resulting algorithm converges linearly much of the time, but it will sometimes go intocycles which are empirically not convergent. In Section 2.5.3, we present a detailed analysisof the convergence of the fixed point algorithm.The Figures 2-11 and 2-12 show the fixed point algorithm both succeeding and failingfor two different random matrices. We have plotted the new value of the variable e againstits value in the previous iteration.2.3.3 The singularities of the spectral portraitWe have identified that the normal space of a matrix A e V is given by the line A + auvH,parameterized by a, where u, v are left and right eigenvectors of A. Thus if A lies in thenormal space of its nearest matrix A E , then A = A+ auvH. There is, however, no obviousway to compute u and v without knowing A in advance. In this section, we present a methodfor this based on the spectral portrait of A.Let a(x, y) = amin(A - (x + iy)I). In the case where a(x, y) is smooth, we have thefollowing theorem about the derivatives of a(x, y) due to Sun (see [60] ).Theorem 2.3.1 Let U(x, y)E(x, y)VH(x, y) = A - (x + iy)I, where an(Xo, yo) is simple, andlet U = (Ul u2 ... Un-1), V = (vi v2 ... Vn-1), and = diag(al,. .. ,n-). Wehave then8x (xo, Yo) = -Real{uHvn},28 y = f(x)26 y=x24,22218-1.6141212 14 16 18 2 22 24 26 28x = e.21205ijt 95S19185181 7517 175 18 185 19 195 2 205 21x = e.Figure 2-11: The fixed point algorithm converging.() 2 2II>,.x = e.Ix = e.IFigure 2-12: The fixed point algorithm approaching a limit-cycle.O (XoYo) = Imag{uHvn},92jX 2 (o, YO) = Real{rHJ1r + 1 H 1 + 21HRir} + Imag{u HVn}2 /Un,02oOxOy (xo, Yo) = -Imag{21HIF r} + Imag UHVn Real uHvn /n,,02 C202 (Xo, Yo) = Real{rHbr + 11 - 21Hr} + Realfvn 2 /u n,where = (cr I - 2)-, - = u2 _ 2)-1, = CjHv, and 1 = V(un.Sun's theorem implies that ax = ay = 0 iff uHv = 0. In this case, the matrix =A - a(o, Yo)uvH E Mxo+,yo is a point on the envelope of the family of MX, i.e. Ac D. Thuswe haveCorollary 2.3.2 If a(x, y) is stationary for some A0 = xo + iyo then A lies on the normalline of the matrix A c D with repeated eigenvalue A0 and left and right eigenvectors u and vwhere A = A - or(xo, yo)uvH and u and v are the left and right singular vectors of A - Aolwith singular value a(x 0 , Yo).Since the stationary points of a(x, y) correspond to normal lines off of D. An examinationof the spectral portrait of A can be used to determine the distance from A to A E D.This would lead to an algorithm which is based on finding the repeated eigenvalue of Afirst, i.e.(xO,yO) = find lowest critical points of the spectral portrait;lam = xO + iyO;[u,s,v] = svd(A-lam*I);Ahat = A - u(:,n)*s(n,n)*v(:,n)';where one finds the critical point of lowest height by Newton's method or some other opti-mization routine over the complex plane.Another result of Sun's formulae is that at stationary pointsdet jrH1 r + 1Hr+ lI 2 _ 121HIgr12.This formula will be useful in the later section on conditioning and stability (Section 2.5).2.4 Nearest A's are basically J2 'sWe have already seen that generically speaking, the matrix A - A is a rank one matrix andA has a J2 structure. Of course, less generic Jordan structures may arise in special cases,but in this section we will show that even when A has a less generic Jordan structure, A - Astill has the form of a generic J2 normal.In Section 2.3.3, we showed that, with perhaps few generic conditions on A, the nearestmatrix, A c D, to A is always given by A - auvH where a = amin(A - AI) is stationary in A,u and v are the associated singular vectors, and A is the repeated eigenvalue of A for whichu and v are right and left eigenvectors.If one were to be cavalier, one might assume generically that the nearest element of Dwas always some defective matrix with a J2 structure. If this were the case, then D can beassumed to be locally differentiable about A. Consequently, there would be a unique normalperturbation that would be of the form uvH where uHv = 0. Thus, the A would correspondto a critical point of the spectral portrait of A - AI since the spectral portrait is a graph ofcrmin(A - AI).However, one might consider the ways in which this generic assumptions could be violated.If A were derogatory or had degeneracy less generic than J2, then there could be multiplenormal directions off of A, and A - A would not be of the assumed rank 1 form. Oneparticular example of where one might really be concerned is in the case where A actuallyhas a J3 structure and A might sit above all of the rank 1 normals of A as shown in Figure2-13 (since a matrix with J3 structure is a cusp point of D, there is no reason to believe thatit has a continuous local normal field).Clearly, one needs to have some lemmae to understand the conditions under which such0 15r01- A005A-0 05-015-0 2-0 25-03-0 25 -0 2 -0 15 -0 1 -0 05 0 005 01 015 02 0.25bFigure 2-13: A possible problem finding a normal from A to A illustrated in characteristicpolynomial space.an event will or will not happen in order to make generalizations which are everywhere valid.The purpose of this section is to clarify these possibilities. We present in this section a seriesof lemmae which will progressively refine the possible form of N = A - A until we find thatit has the form of a normal to a generic J2. The reader who desires to skip the technicalitiesof this chapter is advised to skip this section.It will be useful to take a partial canonical decomposition of A in the repeated eigenvalue.That is, letA= (x X) 0 M ( Y Y)Hwhere (Y Y)H(X X) = I and J(A) is the canonical matrix of the Jordan structure ofA at A.Lemma 2.4.1 Let A be the nearest matrix in E) to A. ThenN = A- A= YWXH,for some W.Proof: We also decompose N = A - A according to the following (slightly twisted) decom-position,NilN = (Y 2) 1one may check that (XN12(X Y)HY)H(y X) = I, and that the displacements YN 11X H, XN 12X HYN 2 1Y H, and XN 22YH are mutually orthogonal.Since A + XN 22 H has the same canonical structure in A as A, it is clear that N 22 =0 for A to be the minimizer. Similarly, by observing that (A + YN 21YH)X = XJ andYH(A + XN 12XH) = yH, we see that N 12 and N 21 must vanish.For the next lemma, we consider the case where J, the Jordan structure of A at A, isderogatory and thus of the form J = j()partitioned into Y = ( Y Y2 ) and X = ( X(2)) and the eigenvectorX 2 ). Decompose N as's are correspondinglyN=(Y Y2) N N2 2 ( X X 2 ) H.Lemma 2.4.2 Suppose that the Jordan structure of A is derogatory. Then N is of the form,N = YNiX Y2N2 2XHProof: Observe that A + X 2WY has the Jordan structure of ((being possibly different from J, still has a repeated eigenvalue A.j(2))which, whileA+X 2WY 1 is closer to A than A if (X 2WY 1 , N) = (W, N 12) is greater than zero for somearbitrarily small W. Thus, N12 = 0.By a similar argument, N21 = 0. 0Finally, we can, in a sense, eliminate the derogatory A matrices by the following lemma.Lemma 2.4.3 If A has a derogatory canonical structure, J, then one of the irreducibleblocks of J, say J(1) is defective, andN = Y1 NI X.Proof: Suppose that all of the irreducible blocks of J were simple, i.e. J = AI. Then, fromthe previous lemma, N = ylnx (where xl and Yl are vectors and n is a scalar).A + pI has a repeated eigenvalue of A + p.A+ p will be closer to A if (pl, N) = Real{pn} is greater than zero for some arbitrarilysmall p. Thus n = 0 and A = A contradicting hypothesis.Combining this with the result of the previous lemma gives the the conclusion. OIWe summarize with the following corollary.Lemma 2.4.4 The nearest matrix AE D to A is reached by a perturbation of the formN = Yp(J(O)T)XHwhere p(x) is a k - 1-degree polynomial with vanishing constant term, and X, Y are the rightand (corresponding dual) left generalized eigenvectors of a canonical k x k Jordan block, J,in the Jordan structure of A, and J(O) is the k x k Jordan block with eigenvalue 0.Proof: Prior lemmae have shown that N is of the form N = YWXH. The similaritytransformation e~Me - ' m E D is closer to A for some arbitrarily small o whenever thecommutator, MA - AM has positive inner product with N.Thus ([M, A], N) = (M, [N, A"]) must vanish for all M. Giving NAH - AHN = 0 whichimplies that WJ(O)T = J(o)TW, and thus W = p(J(O)T). OWe now need only really consider the case where the Jordan structure of A is a singlecanonical Jordan block.Lemma 2.4.5 Let the Jordan block J of A be of size greater than 2. The polynomial p(x)inN = Yp(JT)xH,must be of the form p(x) = axk- 1 + bxk - 2Proof: Suppose that the coefficient, c, of xj (j < k - 2) is nonvanishing. Then we considerthe matrix A + XMYH, where Mkj = s (Mpq = 0 elsewhere).A + XMYH, while not necessarily having the same canonical structure as J still has Arepeated at least twice, since(A + XMYH) ( x l 2 ) =( 1( 2)However, A + XMYH is closer to A so long as (M,p(J(O)T))arbitrarily small s. Thus we must have c = 0.Thus, for all A not in D, N = A - A is of the form that oneeither a J2 or a J3 canonical structure. In what follows, we shalltake values in C then the A - A will be rank 1.Theorem 2.4.6 Let A E D be the nearest matrix withmatrix A with distinct eigenvalues. Then= sc does not vanish forgenerically expects fromsee that if we allow A toa repeated eigenvalue to the givenN = A - A = auvHwhere u and v are the left and right eigenvectors of A and a is a scalar.Proof:Assuming the canonical structure of A is J = Jk(A), we have from the previous lemmae,we have thatN = Y(a(J(O)T)k- 1 + b(J(O)T)k- 2 )XHb : 0.Consider the matrix A + X(a'(J(O)T)k- 1 + b'(J(O)T)k-2)yH, where b' = ky k - 1 and a' =-(k - 1)p k . It is easy to show that J + (a'(J(o)T)k- 1 + x(J(o)T)k- 2 ) has a repeated (J 2 )eigenvalue A + p.if for some arbitrarily small p (X(a'(J(O)T)k- 1 +b'(J(0)T)k-2)YH, N) is greater than zero,A + X(a'(J(0)T)k - 1 + b'(J(O)T)k-2)YH will be closer to A than A.(X(a'(J(O)T)k-1 + b'(J(O)T)k-2)yH, N) = Real{b'b + a'a}. Substituting for a' and b' wehave Real{b'b+a'a} = Real{pk-l(kb - (k- 1)fPa}. If we take [p| < and let select the phaseof p such that pk-lb is positive real, then (X(a'(J(O)T)k- 1 + b'(J(O)T)k-2)YH, N) > 0.Thus, we have A + X(a'(J(O)T)k- 1 + b'(J(O)T)k-2)YH closer to A than A.We must conclude that b = 0 and therefore that N is rank 1 of the form asserted. DThe reader will note that it is not always possible to find such a' and b' if one is constrainedto real matrices, but one can always do so over complex numbers. In fact, the illustration atthe beginning of this section demonstrates this shortcoming of the reals. If this illustrationhad been drawn with complex a and b axes, then it would have been apparent that A wascloser to a complex A.We have now established rigorously that any matrix A with distinct eigenvalues must liealong a rank 1 normal from the nearest A E D. Thus A can be determined by examiningthe singular values of A - AI.What remains are some lemmae concerning the singular values themselves so that itcan be determined what sort of behavior one must look for in the singular values in orderto locate A. The previous section shows that the critical points of Ormin are the points ofinterest in determining A. However, one may also wonder what happens when 'min is notsimple (and hence not rigorously differentiable), and whether or not one must examine thecritical points of the other singular values of A - AI. We assure the reader that this will notbe trouble for us. In Section 2.6, we will be able to show that the critical points of othersingular values will not produce minima, and that, the minimum will occur on a criticalpoint where -min is almost always simple.2.5 Conditioning and stabilityIn this section, we wish to examine the conditioning of the problem of finding the A E Dnearest to A. The conditioning of A is closely related to the behavior of the singular valuedecomposition of A - AI under first order variations. We will find that the matrices withill-conditioned A on D are sitting close to "focal points" of D analogous to focal points ingeometric optics, explained in Section 2.5.2.-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1bFigure 2-14: The normals of the cubic double root variety forming their parabola-like enve-lope. The concave down arc in the plot is another part of the envelope which occurs overcubics with complex coefficients.2.5.1 Sensitivity analysis of A = A + auvHConsider the normal lines passing through a matrix A E 2D. We can parameterize these linesby the real or complex line A = A + auvH, where a is real and the phases on u and v arearbitrary (alternatively, we could fix the phase of uvH and let a be complex). For most A weexpect this decomposition (A = A + auvH) to be locally invertible. We define a focal pointas some matrix A for which the decomposition is not locally unique, i.e., where the Jacobianof this decomposition is singular.Geometrically, we may think of the set of focal points as the evolute of D2, i.e., the envelopeof normals of D. A point which is in the evolute can be thought of as the intersection of twoinfinitesimally separated normals.It will be useful to do some dimension counting in both complex and real versions of D2so that we might better understand where the degrees of freedom are. (We shall always referto one real parameter as one degree of freedom.)-1F 2....i-n s . A-- f t d e p g-3Figure 2-15: The normals of the swallowtail. A faint hood-like shape begins to emerge.Figure 2-15: The normals of the swallowtail. A faint hood-like shape begins to emerge.Figure 2-16: The swallowtail and its hood-like evolute of ill-conditioned problems.The solution to the problem of finding a stationary point A to A can be expressed interms of decompositions of A as follows:A = auvH + A, (2.13)A = auvH I+ UEV H (2.14)A = AI+ UEVH (2.15)where uHv = 0, = diag(al,..., a- 1), E = diag(al,..., o-_, o-) = ( ), U = (U u),V= (V/ v), and U, V are unitary matrices. For the real case, Equation (2.13) decomposesthe n 2 degrees of freedom of A into n(n - 1)/2 for U, n(n - 1)/2 for V, n for E, and 1for A, which leaves exactly 1 degree of freedom extra to satisfy the constraint uHv = 0.For the complex case, A has 2n 2 degrees of freedom, U has n 2, V has n2, has n, and Ahas 2, leaving n + 2 surplus degrees of freedom. Since the complex decomposition is onlyunique up to transformation of the form U, V - UA, VA (A a diagonal unitary matrix),n of those surplus degrees of freedom are ineffective, leaving two to satisfy the constraintsReal{uHv} = 0, Imag{uHv} = 0.Thus, by counting dimensions we see that the decomposition which reveals the nearestAG eD to A is really a singular value decomposition in which one dimension has been removedby the constraint uHv = 0 and one dimension has been added by the presence of the AI term.Our strategy to vary A will be to vary the SVD (as in [60]), and then to use the variationin A to satisfy the varied constraint (uHv = 0) equation.Note, if we were to vary this decomposition for an element of A C D, we set a = 0 inEquation (2.13). In the real case, this eliminates one degree of freedom corresponding to D'sreal co-dimension of 1. In the complex case, this eliminates not only the degree of freedom ina, but also the extra unitary freedom of u and v corresponding to A's complex co-dimensionof 1.Theorem 2.5.1 Assuming that the singular values of A - AI are distinct, the first orderbehavior of the decompositionA = A - auvHis given byA = A + (uHAv)uvH + U(UHit)vH+ au(VHi)HVHwithUH ( HA - I) vVH (a) 4 (a) V H (AH- Ai)where 1() = a(a 21 - Z 2 )-, () = (_[ -_ 2)-1, and A is given by( Oxx O Uxy Areal (Real{p + lHIp lq + qHir + pH4r}Oy Oyy Azmag Imag{lHIp + 1'Hq + qHHir + pHr} Jwhere axx, axy, and ayy are given by Sun's theorem.Proof: Varying the decomposition (2.13) along a differentiable curve A(t), we haveA = A + 6uvH+ aUivH+ aUU1Hwhere iHv + uHi = 0. We will proceed by solving for an expression for iHv + uHi, findingvalues for a and the phase of uvH for which the variation of the constraint is satisfied (it isconvenient for us to take a to be real and to adjust the phases of uvH implicitly).Differentiating (A - AI)v = an and uH(A - AI) = avH, we have that(A - AI)v - &u = ai - (A - AI)iand(AH - I)U - v = aUi - (AH -AI)i.We may substitute A - A = UEV H, obtainingUH(A - AI)v - &en = aUHit - EVHUandVH(AH- AI)u - Tre, = aVHi - EUHiL,... 0 1)H, which can be collected into 2n - 2 equations,(iaS UHl UH(A - (2.16)and 2 identical equations,UH(A - AI)v - & = a(uHit - vH'). (2.17)If a = an is distinct from the other ai then we can invert the matrix in (2.16) obtaining( UH V () T(U) UH(A - AI)vT(a) 4(or) \ VH(AH- AI)u'(2.18)where ((a) = a(a 2 1 - 2 )-1 and T(a) = t(a 21 - t2)-1. Since UHj andskew-symmetric, (2.17) separates into real and imaginary parts& = Real{uHAv}- i(UHiL - vH) = Imag{uHAv},VHV must be(2.19)(2.20)where we have used uHv = 0.We now apply the differentiated constraint,*HV + Uh = (UHit)H + H(VZ) = 0,where we have used v = Ur and u = VI (r = UHv and I = VHu). Substituting (2.18) intothis, we haveA(21H'Fr) + A(lHl + rH(r) = 1H p + lHq + qH4r + pHr,_1_ ~__~~~__~I1I _I11llr_ I / II I YIII 1111 111 1111 111111111where e n = ( 0J)u(2.21)where p = UHAv and q = HAHuH, which one can solve for A solving(u axy Areal Real{lHp + lHbq + qH4Jr + pHIr}Oxy O iyy " mag Imag{lHIJp + 1Hq + qH4r + pHbr}JOne then hasA = A - Real{uHAv}uvH - aU(UHUi)v H - (u(VHi)HVH (2.22)A = A - Real{uHAv}uvH - aU(UHi1)vH - u(VH )HVH oU(uHji - ilHv)UvH (2.23)A = A - (uHAv)uvH- OU(UHi)vH- Ou(VHj)HvH, (2.24)where we have used Equation (2.20) to get (2.24) from (2.23). EOOne can interpret the formula for A as displacement resulting purely from the change inthe phase and length of the normal (second term) plus a tangential variation of A (third andfourth terms) which results from a change in the direction of the normal from A to A. Themagnitude of the tangential part of the variation is ~i 2 + 2.The set of A with A having infinite condition number is given by finding those A forwhich A can be nonzero when A vanishes. In this case, we set the p and q terms in (2.21)to 0 and solve(cxx cxy (ArealUxy y'y Azmagas an equation of o equivalent to21H) lH(o)l + rHI(u)r,which gives the following corollary.Corollary 2.5.2 The nearby A E D to a generic matrix A is ill-conditioned whenever2 -1H + r is small (i.. U y is nearly singular).O-yx /yy2.5.2 Ill-posed problems and principal curvaturesIn this section we explore the differential geometry of D itself, as a co-dimension 1 manifoldembedded in Rn 2 or C2 . We will give an additional interpretation of the focal points of Din terms of principal curvatures.The classical study of principal curvatures is concerned with surfaces in R3 (as in [25]).In this setting, one uses the unit normal field i of a surface M to construct a map from Mto S2 . Studying the local behavior of this map ultimately results in a quadratic form on thetangent space of M at a point. The reciprocal lengths of the principal axes of this quadraticform are called the principal curvatures, their sum the mean curvature and their productthe Gaussian curvature.One way to generalize this approach which is suitable for any manifold of co-dimension1 is to use the normal to construct a natural quadratic form on the tangent space and thenexamine the eigenvalues of this form.We first look at an n-sphere of radius 1 embedded in Rn+1. Geodesics on S n satisfy theequationsdp- Vdv v -vdt Lwhere p is a point on the sphere and v is the velocity, satisfying pHv = 0, and h = p/1. Onemight recognize the derivative of v as the centripetal acceleration in classical mechanics.For a more general manifold of co-dimension 1, one has (at least locally) a unit normalh and geodesic equations of the formdp = v (2.25)dtdvd = -(v, v) (2.26)dtwhere F(v, w) is the Levi-Civita connection. However, for manifolds of co-dimension 1, theLevi-Civita connection must be of the form F(v, w) = a(v, w)i where the aceleration, a(v, w),a-.I.XI?---~PII-----' ~~~-~ll-.--.~.I ~ I)^ -- .-.-.I- -~ ---^-~I~~~~-~ --- ~- ------P~-L ~ - -~ ----~IIIII aa__ --C-.0.020--0.02--0.04--0.06 --0.08 a(v,v)-0 1-012-0.14-0.16 --0.1 -008 -0 06 -0.04 -0.02 0 0.02 0 04 0.06 0.08 0 1bFigure 2-17: In the cubic polynomial space, a particle moves along the double root variety,with velocity v and centripetal acceleration a(v, v).is a symmetric bilinear (not complex sesquilinear) function defined on the tangent space atp and determined by d v + v = 0.If we consider a(v, v) to be an instantaneous centripetal acceleration along a geodesicwith velocity v through p (see Figure 2-17), then we may define the principal curvatures at pas the extrema of the function la(v, v) I over the unit ball in the tangent space at p. As withthe sphere, the center of curvature is located at p - a(,) . We may refer to the reciprocalsof the principal curvatures as the focal lengths along f from p (see Figure 2-18)With this in mind, we now consider the differentiable subset of D. For the remainder ofthis section any matrix A is a point in a differentiable neighborhood of D with decompositionA = J tV H Aand normal N = uvH (uHv = 0). Any element H of the tangent space of A must satisfyuHHv = 0.'ai-0.2--0.4--0.6--1.21--1 -0.8 -0.6 -0.4 -0.2Figure 2-18: The centripetal acceleration ofcan be associated with a circle centered at0 0.2 0.4 0.6 0.8 1the particle moving along the double root varietythe focal point of the double root variety.LLLX--I_~I_ _ __ __ IIC~l-^--_ -----~. .. .0064" 1 NF A---"-- C- - - , _ _~.- -- , ..If A(t) is a differentiable curve in D through A then we differentiate uHHv = 0 to getdd (uHHv) = iHHv + uHHi - a(H, H)dtHit) )Hp + q HfZv= 0= a(H, H),where q = UHHv and p = V7HHv (note that since uHHv = 0 we have Hv = Up anduHH = qHfH).Since A E D we may set a = 0 and employ Equations (2.18, 2.21, 2.24) to see that(tUHit, H'k,)0 p - Arq -with A given byA(21H-lr) = 1H-1lpqHt-lr.Substituting and simplifying, we have(1Ut-lp + qH-1-lr) 2a(H, H) = -2pH-lq +21HE-1 rSince UpvH + uqHVH is a tangent vector for any p and q, we may consider a(H, H) as afunction of 2n - 2 independent p and q vectors.We see that the quadratic form a(H, H) can be written as zTWz, wherez = (2.27)qand0-t-11+ 21T r C-l ri-1T(r T t-1 (2.28)By a simple orthogonal similarity transformation we haveQWQT = Z + aaT,-1 -10(2.29)where Z = - and a =( 1 , a real diagonal matrix plus a0 E-1 2 r (-1(r -rank-1 matrix.We now seek the extrema of zTWz where HI 2 = Ix12 + 1y12 = 1. Note that WT = W.Varying IzTWzI we see that the condition for z to be an extremum of IzTWzl isWz = ( ,a somewhat unusual eigenvalue equation over C.Over R this reduces to the problem of finding the eigenvalues of a rank-1 update, Z+aaT.An eigenvalue ( of this matrix then satisfiesaT((I - Z)-la = 1,which, in our terms, is((r + l) T t- 1 (r - l)TE-1) + 1)1) ( r + 1 = 41T t- 1 r.O ((1 - t-)-1 -(r 1)Multiplying this out and simplifying gives-1 1 -1 1 1rT (I - 2)-1r +T -_ 2)-11 + 2 1 T( -_2 ) -r = 0.( C2 ( (2 (2-- 1 -- 1 -- 1r T(l )r + 1T4( )l+ 21T ( )r = 0,where )(a) = a(a2I -_ 2)- 1 and I(a) = t(a2I- t2)-1. The center of curvature is thenA + auvH where a = 1 We may now interpret the vanishing of the second real-axisderivative at a stationary point of the spectral portrait as equivalent to being on a focalpoint of the real double eigenvalue variety.Over C, we solve QWQH(Qz) = (Qf by first observing that given an eigen-pair z,( , theeigen-pair zei2Z , (e-iz is also a solution. Thus, while acknowledging that A can have anyphase. we may restrict ( to the reals for purposes of computation, in which case we solve,QWQH-(I QWQH) ()Following a derivation of the usual rank-i update formula, we rewrite this asz-(I (zf)0aa)Letting k = aTz we have(z Z(Z2 _ 21)-la( ((Z 2 (2 )-lak a+((Z2 - 2)-1 aZK-(((z(aT((aTZ(-(I z akZ akZ2 - 2Ia)-(lZ2 - (2I)-la KNote that the offjugates, and that (( diagonals are real, while the diagonals are each others) must be a fixed point of this matrix. A necessaryk/complex con-and sufficientcondition for such a k to exist isdet (aTZ(Z2 - (2 I)-a +- 1aT( Z 2 _ (2) -l aaT((Z 2 _ c2)-1rCTZ(Z2 _ (2i)-1 += 0,or equivalently,laTZ(Z 2 _- (2+)-la + 11 = laT((Z 2 - (21)-lal.Substituting and simplifying we have-1 -1 -1 1r ~( )r +t 1F( )1=|Fi/ rThus we have seen that an additional way to interpret these focal points is to imaginean osculating circle of radius !, centered at the focal point and passing through D along a=0S0- 0.aaT+ ( 0- (I z)Z ) (fprincipal direction.2.5.3 Stability of the fixed point algorithmThe fixed point algorithm can be viewed as a series of guesses at the repeated eigenvalueA = limi,, ei, where ei+l = f(ei). In iterated systems of this sort, convergence is to somefixed point x0 = f(xo), and the linear multiplier governing convergence is f'(xo).If f'(xo) > 1 then the system diverges, either to another fixed point or to a limit cycle. Iff'(xo) < 1, but close to 1, then convergence can be very slow. We therefore need to examinef'(xo) to understand the stability of this algorithm. We linearize f(x) about a fixed pointA.f(A + 6A) = cluster(A - (a + 6r)(u + 6u)(v + 6v)H)where u + 6u, v + 6v, and a + 6o are the singular vectors and value of A - (A + 6A)I and u,v, and a are the singular vectors and value of A - AI. This may be rewritten asf (A + 6A) = cluster(A - 6auvH - aSuvH - auvH).Assuming that the cluster function simply averages the two nearest eigenvalues, we havea first variation givingf(A + 6A) = A + -tr((yl Y2 )H - 6 uvH - uvH - aU6vH) (X1 2 )),where Y2 and x, are the left and right eigenvectors of A, yHA = yH, and Ax 2 = X1.Since v is a right eigenvector and u is a left eigenvector, we may take xl = 1 _lV,= U X2 -lr, and YI = U -1 1. Thus we see that6 f (aOiH-- UHu + a6vVt-lr).21Ht 1rdll l IIIY I1141 N I, id'. 1. 1 Ji NI Y NN Eno NM I@IbI , JSubstituting from Theorem 2.5.1 (with A fixed), we haveSf = 1 (6A21H(E - l + T)r + 6A(lHl + rHDr)),21HE-1rsimplifying to16f = 6A + (6A21Hr + 6A(IHl + rHr)),which is the linearization of the iterated map f on a small perturbation 6A.The behavior of this map will determine the nature of the fixed point. We recognize ele-ments of the second derivative of the spectral portrait in this formula. Using Sun's formulae,we have1 UxcyyY9QC + Y+yy6f 6A- + - r (6A( - iaxy) + 6A( 2 ,2lH-r 2 2from which we may assert the following theorem.Theorem 2.5.3 The fixed point algorithm iteration linearizes torea a rmag ( xy \ (6Areal (2.30)6Aimag i+1 1H -rimag -- IH-lrreal ay y 2.30)imagWe can make a number of qualitative statements about Equation (2.30). When theHessian of Omin is singular (i.e., near a focal point) this iteration can stagnate. Secondly,if IH-lr is small (corresponding to the Puiseux expansion in Section 2.3.1 having poorbehavior), then the iteration is more likely to diverge. Lastly, since(U Xz Y (yH-lrreal H -1rimagUaxyy YlH-lrimag -Ht-lrreal)as a -+ 0, we see that fast convergence is assured for small enough a.2.5.4 Sensitivity of the critical points of the spectral portraitWith results which will be derived in Section 2.6, the problem of finding the nearest A E Dto A can be reduced to locating the critical points of the spectral portrait of A, a(x, y) =amin(A - (x + iy)I)) as suggested by Section 2.3.3. We now turn our attention to thesensitivity of this problem, investigating the sensitivity of the critical points of a(x, y) tosmall perturbations in A. To first order, a perturbation 6A in A will cause a perturbationin the gradient of a(x, y) = (ax, ry) at the point (x0 , yo), which we may call (6ax, Say). Byinverting the Hessian (given by Sun's formulae) on this perturbation, one then has the firstorder correction to (x0 , yo), which we write as (6x 0, 6So).Theorem 2.5.4 Let (x0, Yo) be a stationary point of a(x, y) and 4, I, r, and I be definedas above. Then let p = UH6Av and q = VHSAHu. Then6xu(xo, yo) = Real{rHop + rHTq + pHHq1 + qH4l}6ay(xo, Yo) = Imag{rHDp + r4q + p H T + qH4l}.Proof: From Sun's theorem, we haveaX + iay = -vHuand thus6Ux + i6S, = -6vHU - vHu.Since uHv = 0, we may substitute v = UUHv = Ur and u = VVHu - V1,6a + i6Su = -- rHH6u - (VH6v)Hl.From Theorem 2.5.1 we have a formula for UH6v and VHSu, which we substitute to giveax + iay = rHp + rHq + pH 1 + qHl.OFrom this lemma and Sun's theorem, we see that the condition number on the criticalpoints of the pseudospectra is given by largest singular value of the 2 x (4n - 4) real matrixO-yx O-yy 6a ywhere 6a, and 6ay are the real linear functions of the real and imaginary parts of p and qdefined by the lemma.2.6 Nearest A's come from simple aminIn Section 2.5.2 we derived expressions for the centripetal acceleration of a geodesic along D.We will now use this machinery to give us the prove the claims made at the end of Section2.4, namely, that the only critical points of singular values one need examine are those ofamin and not any other ak. Once again, the reader who wishes to skip over technicalities isadvised to skip this section.First we will make some fairly general derivations over co-dimension 1 manifolds. Thesederivations will relate the extrema of centripetal acceleration through a point to the characterof that point as a minimizer of euclidean distance. Then we will apply these results to thespecific problem of minimizing I A - All.We first consider a general manifold M of co-dimension 1 in 1nV or Cn. Suppose thatp C M locally minimizes the function1 1D(p, q) = -(p- q,p-q) -d(p,q) 2 (2.31)2 2for some point q in ER or C". Let p(t) E M be a geodesic passing through p at t = 0 satifying(2.25, 2.26) with unit velocity v ((v, v) = 1). Differentiating (2.31) with respect to t, we havedtD(p(t),q) = Real{(p-q,v)}d2dt D(p(t),q) = Real{(v, v) + (p - q, -F(v, v))}d 2dt2D(p(t),q) = 1 - Real{(p - q, f(v, v))}.Setting the first of these equations to zero at t = 0 for all v gives p - q = au where a is ascalar (Ja = d(p, q)) and h is the unit normal of M at p. This givesd2dD(p(t),q) = 1 - Real{(au, F(v,v))} (2.32)d2dt2 D(p(t), q) = 1 - Real{aa(v, v)}, (2.33)where we have substituted the Levi-Civita connection, F(v, v) = a(v, v)i.In order for p to be a stationary point of D(p, q) it is necessary that q = p+a U. However,in order for this stationary point to be a local minimum, the second derivative must benonnegative for all v. Since a(v, v) is bilinear in v we may select the phase of v so that thecondition (2.33) becomes1 - alIa(v, v) I> 0,or1- > a(v, v), (2.34)where we understand (2.34) as a comparison between magnitudes. Since (2.34) must betrue for all unit v, and since the principal curvatures of M at p are given by the extrema ofla(v, v)l over the set of all unit v, we haveLemma 2.6.1 If d(p, q) is minimized by p then it is necessary that the focal lengths of M atp be greater than or equal to d(p, q); i.e., if there is a focal point on the line segment p + sau(0 < s < 1) connecting p to q, p cannot minimize d(p, q).This relation of focal points to the characterization of the stationary points of d(p, q) canbe illustrated clearly in two dimensions. In Figures 2-19 and 2-20, we use a hemisphericalsurface to illustrate these lemmae.We may now move to the specific case of minimizing d(A, A) for A E D.Lemma 2.6.2 LetA =Uv H +ASfocus#q4local minimumFigure 2-19: Below the focal point of a hemispherical surface, the distance minimizing p onthe hemisphere to a given q is at the bottom of the surface.focuslocal maximumFigure 2-20: Above the focal point of the hemispherical surface, the bottom point of thehemisphere is now a maximum of distance (the nearest p on the hemisphere to q has becomeeither of the end points of the hemisphere).be an element of D. Let a(H, H) be the centripetal acceleration defined on the tangent vectorsof A by Equation (2.27)., Then, taking maximums over the unit tangent ball ((H, H) = 1),1 0 thenI 107n-1 01 + aaT)z| 0. Thus,1 < zT( (1Un1 ( 10 + aaT)z Oan(A - Aol)) then(using singular vectors u and v) A = A - kuvuH cannot locally minimize I A - All over allA E D.Proof: Since the smallest focal length (reciprocal principal curvature) of A is less than on,it is impossible for d(A, A) = Uk > n to be locally minimized. OWe can also follow this sort of reasoning to establish that the points in the spectralportrait which are crossing points of singular values will usually not produce local minimizersof distance.Corollary 2.6.4 Generically, oa,(A - AI) are simple in the neighborhood of A0 .Proof: Generically, if Unl- = a, at A = A0, then A cannot be a local minimum of d(A, A).2.7 Triple eigenvaluesThe variety of matrices with triple eigenvalues, T, is a subvariety of D). It is the varietyof cusp points of D, the points where D) is not smooth. This might lead one to think thathigher order derivatives of Omin might shed light on the distance to T. This is not quitethe case, however, though some bounds are possible. This section will investigate relationswhich can be made between T and the focal points of D.It is easy to see that the intersection of D with its evolute is T. For A E D with repeatedeigenvalue A, we may take xz = v, x 2 = UZ-1r as the right eigenvector and generalizedeigenvector and y2 = U, 1 = V-11 as the left eigenvector and generalized eigenvector, with(Yl y 2 )H(xI x 2 )=lH -r( 0Note that if lHi-lr = 0 then the left and right eigenvectors cannot be duals and thus theremust be additional right and left generalized eigenvectors of eigenvalue A. Let us now supposethat A has a focal length of 0 (thus A is in the evolute of D). Then, by the definitions ofprevious sections,IrHI(O)r + 1H(0)11 = 121H (0) r.However, the LHS is zero. Thus we have 1HJ(O)r = lHY-ir = 0.This is equivalent to the condition det )x axy = 0. Thus T can be thought of asthe locus of points satisfyingamin(A- A) = 0dA- min(A - A) = 0for some A, with the additional constraintd 2dA2 min(A - A) = 0,in contrast with the focal points which satisfyamin(A - A) = 0damin(A - A) = 0d2dA2 m i n ( A - A) = 0.Unfortunately, while the first derivative of omin is independent of the value of Omin, thesecond derivative is not, preventing us from being able to use this set of equations as amethod for finding the nearest element of T.We can bound the distance to the nearest AE T from A if A e D is ill-conditioned. Itis easy to show that one can always find a A within |A - AiI(21HE-r) of A, where Ai is aneigenvalue of A distinct from and closest to the double eigenvalue A.68Chapter 3The Computation of Nearest k-foldEigenvalues3.1 IntroductionA curious formula presented by Malyshev [49] has the property that it will, in some cases,determine the 2-norm distance of a matrix to the set of matrices with a given double eigen-value. These results can be extended to a formula which gives the distance to the set ofmatrices with a given eigenvalue of any multiplicity.There is a sense in which our extension (when successful) decomposes the given matrixA into A = A + aN where A is a matrix having Jk(A) structure for some given k and A, N isa member of a set of transverse directions dependent only upon A, and a is a real number.We can view the set of N for a given A as a set of normals to the set of all Jk (A) matrices,although such an identification is not rigorously possible without an inner product betweenmatrices (which the 2-normed matrix space lacks).The formula we present does not always produce the aforementioned decomposition.Examples can easily be generated in which our formula fails. While we have reason tobelieve that the nearness of A to the set of Jk(A) is an important factor, it is currentlyunclear under what conditions our formula is guaranteed to succeed.- hIter . . ................ . 1 I~ii itI~ IYii II , i hi1''111111 oilII3.2 Malyshev's formula for dist(A, J2)We let dist(A, Jk(A)) denote the 2-norm distance from A to the set of matrices with a k-thorder Jordan block for the eigenvalue A. A tight estimate of this distance can be of greatvalue in regularizing the numerical computation of the Jordan canonical form of a givenmatrix [40], as well as understanding the conditioning of the eigenvalues of A [67, 68].Since dist(A, Jk(A)) = dist(A - A, Jk(O)), we need only consider nilpotent Jordan struc-tures in what follows, setting Jk = Jk(O). Malyshev, in [49], proposed the following formulafor the 2-norm distance of a given matrix A to the set of matrices with a double 0 eigenvalue:dist(A, J2) = maxo 2n_( , (3.1)where 92n-1 is the second smallest singular value. In what follows, we see that this formulahas a natural generalization to more complicated Jordan structures. We do not give a proof ofthe correctness of the generalized statement, and in fact it is does not take long to randomlygenerate a failure. However, it may be that there is some appropriate generalization of (3.1)which will reliably give dist(A, Jk) for any A.3.3 Formula for dist(A, Jk)Let A be an n x n matrix and B be k x k. We consider matrices of the form I 0 A - BT® I,where 0 is the Kronecker product ((XT® Y) . vec(Z) = vec(YZX)). In (3.1) A is the givenmatrix and B is the matrix with y indeterminant.0 0In general we will consider B to be a k x k strictly upper triangular matrix of indeter-minants. Then we examine the generalizationdist(A, Jk) = max a(IO A - BT I), (3.2)Bwhere a = Onk-k+1 (singular values in decreasing order). This equation, unfortunately, isnot always true. However, sometimes it is.Lemma 3.3.1 Let A have the eigenvalue 0 with multiplicity of at least k. Thenank-k+l (I0 A - BT0 I) = 0for all k x k strictly upper triangular B.Proof: Let X be a (n x k, full rank) matrix of right eigenvectors and Y a (n x k, full rank)matrix of left eigenvectors such thatAX-XJk = 0YHA - jYH = 0.We proceed by showing that the equation AV - VB = 0 has at least k independentsolutions. Assuming that B is generic (i.e. all elements on the first superdiagonal arenonzero), we have B = R- 1JkR for some upper triangular R andAV - VB,for V = XR. Since X is full rank, V is full rank.Let V (j ) = VB j -1 (where 1 < j < k) thenAV (c ) - V ( )B = 0for all j. Since V is full rank and the superdiagonal elements of B are all nonzero, the vectorsvec(V( j )) are linearly independent solutions of the equation (I 0 A - BTO I)v = 0.Since a is a continuous function in B, it must vanish everywhere if it vanishes on genericB.Lemma 3.3.2 9nJk-k+1(I 0 A - BT® I) is maximized by some finite B.Proof: It suffices to show that ank-k+l(10 A - sBT® I) -+ 0 as s - c00 for all B such thatl1B112 = 1-, r - ~-x-- -ll i -iJ i l----" " .. ... . . . . lN lill I llmll0lhkIf A is nonsingular then we writeUnk-k+l(I 0 A - sBT® I) = 1/9k((I 0 A - sBTO I)-1)Let m be the largest power of BT such that I I(BI)m I 0. We may substitute(I 0 A - sBT® I)-1 = (I A-1)(I 0 I + sBT A -1 + --- + (sBT A-1)m),to give (by perturbation theory)ak((IA - sB TI) - ) > S'~k((BT)m A-m-1) _ m- I I(B)m-1l A-m 2-- -- I9A- 1 .Thus we see that as s - oc, ak((I 0 A - sBT& I)- 1 ) _ 00.If A is singular then we may modify this argument slightly. For any C > 0, let A, bea nonsingular matrix such that A - A,1|2 < e. Then by the above argument there existsso > 0 such thatUnk-k+1(I o AC - sBT I) < Efor all s > so. SinceOnk-k+1(I 0 A - sBTO I) < Unk-k+l(I 0 Ac - sBTo I) + C,givingank-k+l(I 0 A - sBT® I) < 2e.for s > so.Theorem 3.3.3 Let dist(A, Jk) be the 2-norm distance from a matrix A to the set of matriceswith 0 as an eigenvalue of multiplicity at least k. Then dist(A, Jk) < dist(A, Jk).Proof: Let A - A = A where A has at least k zero eigenvalues. By lemma 3.3.1,O'nk-k+l (I o - BTo I) = 0for all B. Substituting A = A - A we haveUnk-k+l(I A - I A - BT I) = 0,and thusank-k+l(I 0A - B T I) < IAll 2for all B. Since IIA 12 = dist(A, Jk) we have dist(A, Jk) = dist(A, Jk).DWe see that dist(A, Jk) is always an underestimate of dist(A, Jk). If one could exhibit a Asuch that A-A = A for some A with Jk(O) structure having the property IlAIl 2 = dist(A, Jk)then optimality would follow.As we shall see, finding such a A is often possible.3.4 Finding AWe assume that Onk-k+1 (10 A - BTo I) is simple about the maximizing B. We also assumethat the elements of the first superdiagonal of B are nonzero (i.e. B is nonderogatory). Thereader should be warned that these are very questionable assumptions, and that in numericalexperiments one can readily observe violations of both. When these assumptions are made,the construction of A is straightforward. When these assumptions are violated, it is unclearhow to construct A.If one forms n x k rectangular matrices from the vectors u and v by takingU = ( 1 2 ... Uk)andV = ( v v 2 ... Vk),whereandthen the definition of singular vectors implies the following identities:AV - VBUHA - BUH= aU= aV H.One consequence of the above identities is thatu(VHV - UHU) + (BUHV - UHVB) = 0. (3.3)Since B maximizes a, a is stationary with respect to variations in B. The variation of a(when simple) with respect to B, by Sun's formula [60], is given bytr(UHV6B) = -6a.Assuming that a is simple and is stationary with respect to variations in B, we will havetr(UHV6B) = 0.Thus, UHV must be upper triangular. Since the first term in (3.3) is symmetric and thesecond term in (3.3) is strictly upper triangular, both terms must vanish, givingUHUBUHV= VHV= UHVB.(3.4)(3.5)From (3.4) we have we have the following singular value decompositions for U and V:U = PSZHV = QSZH.(Where we have taken the "economy sized" SVD, i.e., all entries on the diagonal of S arenonzero, S is p x p, and Z is k x p for p < k.) We may then construct a matrix A = PQHwhich has the propertiesAV = UUHA = V.From these properties one can then construct a matrix A = A-aA, having the propertiesAV-vB = 0UHA - BUH = 0IIA-All = o.What remains is to show that A does have a 0 eigenvalue of multiplicity k.Since B is nonderogatory, B is similar to a canonical Jordan block of size k by someupper triangular similarity transformation r.010 ... 0r-4Br = J= 0 0 1 ... 00 0 0 ... 0/The r matrix is not uniquely determined. Any matrix of the form r 2 = rlp(Jk), wherep(x) is any polynomial with nonvanishing constant term and rl = r will also satisfy Br2 =r2 Jk(0). We will therefore leave r 2 undetermined until later, when it will become clear how toexploit this freedom. For now, the reader should just keep in mind that the ri are similaritytransforms of B and that they are always upper triangular..^ _I--.~LI-_ilil -II^ ----~-.1~ -IYII~PCIC~---^IC~---We can employ rl and r2 to getA(Vri) - (Vrl)Jk(rT21UH)A - Jk(T21 UH) - 0,or the canonical chain conditionsAX- XJ =yHA - YH(3.6)(3.7)= 0,where X = VrI and Y = Ur2HWe turn our attention to (3.5) which, in terms of X and Y, can be written asr 2YHXrIlB = BT 2YHXrl.Since B = rlJkrl = r2Jkr 1 , we haveY HXJ k JYk HX.Thus YHX must be an upper triangular matrix with constant diagonals i.e. YHX = q0o +q1Jk + q2J + ... + qk- 1 J k- . Since r2 can be written as r 2 = rlp(Jk), we may fix r 2 suchthat YHX = Jk-m where m is the rank of YHX. For the remainder of this section, we taker2 to be fixed in this manner.3.4.1 Generic success: tr(UHV) -/ 0Substituting from the definitions in Equations (3.6) and (3.7), we havetr(UHV) = tr(p(Jk)Jk-m).Thus if tr(UHV) :- 0, we may infer that m = k and YHX = I, and may interpret X and Y asthe dual right and left eigenvectors of A for a Jordan block Jk in the canonical decompositionof A.3.4.2 Nongeneric success: tr(UHV) = 0Lettr(p(Jk)J - m ) = 0.Since p(Jk) is invertible, it must be that m < k. We will delay the proof that X and Y arelinearly independent to first discuss the interpretation of this case.Suppose 0 < m < k. Then we must conclude that A actually has Jordan blocks of sizeJ2k- m (0) with right and left eigenvectorsP = (pl P2 P2k-mr)Q=(q, q2 ' q2k-m),where pi = xi and qk-m+i = yi for 1 < i < k. This follows from a fairly obvious considerationof the canonical decomposition.If m = 0 then we can say that A has a Jordan block with eigenvalue 0 of order greaterthan or equal to 2k. By taking repeated powers of a regularized A-1 on the X and Y matricesone could determine the size of the blocks and fill in the rest of the eigenvectors.Thus, when a is stationary in the B variables and B - Jk, then A must be a matrix thathas Jordan structure Jm(O) for some m 2 k.These conclusions rest on X and Y being full rank, however. For the remainder of thissection we will establish that they are (still assuming that B is nonderogatory).Define rank(X) = rank(V) = ji and rank(Y) = rank(U) = j2. By (3.4) ji = j2 = j.Since the U and V matrices are actually the unit vectors, they must have at least rank 1.We now assume 0 < j < k and show that this leads to a contradiction.By the chain conditions, one hasX1 = X -- =k- j = 0andY+1 = "'" = yk = 0.From the definitions of X and Y, we haveV 1 = - = vk 3 = 0andUj+1 - - - - = Uk = O.However, from (3.4), this also givesU1 = - - - = Uk-j = 0andVj+1 = - - -= .k = O.Thus the rank of U and V must be max(0,j - (k - j)) which cannot be satisfied for any jwhere 0 < j < k, a contradiction.3.5 The 2-norm "normals"In the Frobenius norm, which is derived from the inner product (A, B) = Real{tr(AHB)}, itis possible to define the normals to the various repeated eigenvalue varieties at their genericpoints from the equations (N, AX - XA) = 0 and (N, xiy) = 0 for each distinct eigenvalueAi of A with right and left eigenvectors xi and yi. For a generic point A on the variety ofmatrices with a Jk(A) block (A fixed), this givesN = Yp(JkT(o))XHwhere p(x) is a k - 1 degree polynomial in x, and X and Y are canonical, dual right and lefteigenvector matrices satisfyingAX-XJk = 0YHA - JkYH = oYHX = I.From this it is clear that the set of normals to A is a k dimensional vector space, comple-menting the n2 - k dimensional tangent space at A, and we see that if the nearest Jk matrixto A is A thenA = A + Yp(JT(0))XH. (3.8)Additionally, we may recognize that if p(x) has nonvanishing constant term, then it canactually be absorbed into the Y matrix, giving,AX-XJkk = 0YHA_ jk YH = 0YHX = p(JA).Since the third equation is implied by all X and Y satisfying the chain conditions, we seethat we may write any rank k normal as N = YXH where Y and X are any left and righteigenvector matrices satisfying the chain conditions. This is a subtle shift in point of view.The vectorial nature of the normal space has been obscured. Also, since p(x) could have avanishing constant term, we can only absorb p(J[T(O)) into Y by taking X and Y to be n x mmatrices (m = rank of p(Jk)) satisfying reduced chain conditions,AX - XJm(0) = 0YHA - jm()YH = 0Y HX = P' (Jm(0)),where P'(x) = X2(k-m)p(x), allowing us to write any rank m normal as YXH for matrices Yand X satisfying both reduced chain conditions and the additional requirements of the thirdequation.Unfortunately, the matrix 2-norm does not arise from an inner product. It is therefore notnecessarily possible to construct normals and investigate decompositions similar to (3.8). Yet,when the Malyshev method does give a simple maximum singular value with a nonderogatoryB, there is a decomposition of the formA = A + aN. (3.9)We then seek to understand how to construct, for a given generic A with a Jk block, the setof all possible N such that A is the nearest matrix constructed by the Malyshev method toA = A + aN for some or.For a fixed A assume that there exists an A such that A = A + aN, constructed as perthe Malyshev method. We have n x k matrices, U and V satisfyingAv-VB = 0UHA - BUH = 0UHU = VHV = Rfor some upper triangular, nonderogatory B, and the normal N given by N = UR-1VH. Wemay decompose B = rJkr- 1 obtaining for X = Vr and Y = Ur - H the chain conditions on Xand Y with YHX = p(Jk) for some polynomial p(x) assumed to have nonvanishing constantterm.Working this construction backwards, we see that given X and Y satisfying the chainconditions for A, our task is to construct Xr - 1 = V and YrTH = U such that UHU = VHV.This r determines B = rJkr- 1 and N = Y(YHY)- 1 (rHr)-lXH N = YrHr(XHX)-IXH. Tofind r we first solve for rHr and then perform Cholesky factorization. rHr must satisfy(XHX)(rr)- 1 = (rHr)(yHY). (3.10)If we let X = Qu be the QR decomposition of X with upper triangular u and Y = P1 be theQR decomposition of Y with lower triangular 1, we have rHr = -1(luHulH) - H satisfying(3.10), from which we constructN = P(luHulH)-luHQH(note that luHulH = (uH)(luH)H is symmetric positive definite, allowing us to take theunique symmetric positive definite inverse square root). This cumbersome expression can bereplaced withN = PWHQHwhere W is defined as the polar matrix in the polar decomposition ulH = WS (S =(luHulH) ).Since X and Y have not been assumed to be dual, we may replace Y with Yp(JT(o)) tomake this nonuniqueness manifest. This does not affect the P and Q matrices which resultfrom the QR factorizations of Y and X, but modifies W by giving the equationWS = up(J)l1H. (3.11)We see that W depends on the k parameters of p(x). However, since W does not depend onthe absolute scaling of the coefficients of p(x), we see that the set of W is of dimension atmost k - 1. The set of all N is the cone formed by N = aW for all a, and thus is expected,in the generic case, to have k dimensions.For comparison, we may rewrite the Frobenius normals asN = PWHQHwhere W = up(Jk)H. It is apparent that in some sense the 2-norm normals selected byMalyshev's method correspond to the Frobenius normals, but with a polar decompositionas an additional step, and it is that polar decomposition which destroys the vectorial natureof the normal space.Lastly, one can show that the lower rank normals can be written asN = PWH Q",where WS = up(Jm(O))lH and P1 = Y, Qu = X are the left and right partial eigenvectormatrices. In contrast to the Frobenius normals, the polar decomposition step in the 2-normnormals makes the normals of different rank in the 2-norm case separate from the normalsof lower rank, i.e. they are not the limit points of the set of normals of higher rank.3.6 A possible rationale for the 2-norm normalsWhen one has an inner product, one can derive the condition for normality by I N + T IFhaving a local minimum at T = 0, where T is any first order variation in a direction tangentto A. Expanding the norm in terms of the inner product and differentiating gives the familiarformula(N, T) = 0.Although we cannot define normals with the inner product against tangent vectors, wecan, however, adopt the condition of IN+TI2 having a local minimum at T = 0 and see thatthe set of N defined in the previous section satisfies this requirement. For rank k normals,we may take W = ulH. Substituting, we see that we must examine IPWHQH + TI for alocal minimum at T = 0 (where we have suppressed the subscript on the norm). We canwrite T = PT 1 QH+ PT 12 (Q)H+ PIT 21QH+ PlT 22(Q) H. By the unitary invariance ofthe 2-norm, we haveIIPWHQH+ TI = IIWH + T 11 I,where the other parts of T, being infinitesimal do not contribute to the norm.The tangent space of A is independent of the norm being employed. Thus, since everytangent of A must satisfy (T, Yp(JT(0))XH) = 0 for all p(x). Thus, we have(PHTQ, lp(JT(0))uH) = Real{TlUp(Jk)lH} = 0. (3.12)Using unitary invariance again, we haveIIWH+ T11I = II + T11WII.From the perturbation theorem of singular values, the singular values of I will perturb to firstorder by the eigenvalues of the symmetric part of T11W. Thus, the 2-norm is guaranteed toincrease if the symmetric part of T 11W does not have all positive or all negative eigenvalues.Consider the equationIxH(T11W + WHTH)x = RealxHT11Wx} = 0. (3.13)2(3.13) is satisfied by some x - 0 if and only if the symmetric part of Tl1W does not have allpositive or all negative eigenvalues. Expanding W = ulHS- 1 and letting y = S-2x, (3.13) isequivalent to finding a y # 0 such thatReal{yHS TliulHS- y} = 0. (3.14)But, we haveReal{trS TiulHS-} = Real{trTiiulH} = 0,by (3.12). Thus, the sum of the eigenvalues of the symmetric part of SXT 1 ulHS- 1 vanishes,so we may find y satisfying (3.14), and then an x satisfying (3.13). We see, that III+Tx WII >1 for all tangent vectors to first order.In conclusion, we see that the perturbations selected by the Malyshev method are all2-norm minimizers. Thus, they can be considered "normal" even though the 2-norm doesnot have an associated inner product. We also see that these normals are determined only byproperties of A. In this sense we say that the normal bundle to A is given by Yl-1WHu-HXHwhere u and I are the upper triangular and lower triangular matrices of QR(QL) decompo-sitions of X and Y, and W is the polar part of up(Jk)lH.3.7 Concluding remarksInvestigation of the properties of formula (3.2) has not yet revealed the conditions for beingable to construct A after maximizing Unk-k+1. One can show for any A with Jk(A) structureand normal N, that (3.2) recovers A from A = A + aN when1 p)matrices Y with orthonormal columns such as the orthogonal Procrustes problem.Furthermore, some functions F may have the symmetry property F(Y) = F(YQ) for allorthogonal (unitary) Q with a specified block diagonal structure, which causes some searchdirections to have no effect on the value of F. For example, if A is a symmetric matrix,n > p, and F(Y) = tr(Y*AY), then F(Y) = F(YQ) for all p x p orthogonal Q. Moregenerally, some functions F have the property that F(Y) = F(YQ) where Q is orthogonal(unitary) with the block diagonal formQi 0 ... 0Q= 0 Q2 --- 00 0 ..- QPwhere the Qi are orthogonal (unitary) matrices. More complicated problems with blockdiagonal orthogonal (unitary) symmetry can arise in eigenvalue regularization problems (anexample of such a problem is the sample nearest Jordan block problem solved in the exam-ples).This chapter is organized as follows. Section 4.2 will discuss the basics of calling thetemplate code. Section 4.3 discusses the objective functions and their derivatives for theexample problems explored in this chapter. Section 4.4 contains sample runs of instancesof the examples problems. Section 4.5 explains the code structure and the places where auser may wish to make modifications. Section 4.6 will cover some of the basic mathematicsconcerning the geometry of the Stiefel manifold.4.2 Calling sgminThe matlab templates are ready for immediate use or as a departure point for generalizatione.g. problems over multiple variables with orthogonality constraints, or code optimizations.In the simplest mode, the user need only supply a function to minimize F(Y) (and firstand second derivatives, optionally) in F.m (in dF.m and ddF.m) and an initial guess YO. Thecalling sequence is then a single call to sg_min (named in honor of Stiefel and Grassmann).Efopt, yopt] = sg_min(YO)For example, if the function F.m has the formfunction f=F(Y)f=trace( Y' * diag(1:10) * Y * diag(1:3) );we can call sgmin(rand(10,3)) which specifies a random starting point.We strongly recommend also providing first derivative information:function df=dF(Y)df = 2 * diag(1:10) * Y * diag(1:3);The code can do finite differences, but it is very slow and problematic. Second derivativeinformation can also be provided by the user (this is not nearly as crucial for speed asproviding first derivative information, but may improve accuracy):function ddf=ddF(Y,H)ddf = 2 * diag(1:10) * H * diag(1:3);A sample test call where F(Y) is known to have optimal value 10 is>> rand('state',0); . initialize random number generator>> fmin = sg_min(rand(10,3))iter grad F(Y) flops step type0 1.877773e+01 3.132748e+01 2192 noneinvdgrad: Hessian not positive definite, CG terminating early1 1.342892e+01 2.011465e+01 146438 newton stepinvdgrad: Hessian not positive definite, CG terminating early2 1.130914e+01 1.368046e+01 290568 newton stepinvdgrad: Hessian not positive definite, CG terminating early3 5.974785e+00 1.063044e+01 434822 newton step4 1.135397e+00 1.006835e+01 681414 newton step5 5.526059e-02 1.000009e+01 929492 newton stepTable 4.1: A short list of the optional arguments for sgmin6 5.087895e-05 1.000000e+01 1203160 newton step7 1.706261e-05 1.000000e+01 1518165 newton stepfmin =10.0000The full calling sequence to sgmin is[fopt, yopt]=sg_min(Y0,rc,mode,metric,motion,verbose,gradtol,ftol,partition),where YO is required and the optional arguments are (see Table 4.1):* rc The rc argument specifies whether the matrices will have complex entries or not.Although most of the code is insensitive to this issue, rc is vital for counting thedimension of the problem correctly. When omitted, sgmin guesses based on whetherYO has nonzero imaginary part.* mode The mode argument selects the optimization method will be used. 'frcg' selectsFletcher-Reeves conjugate gradient (popular in the electronic structures community),and 'newton' selects Newton's method with a conjugate gradient Hessian inverter (fol-lowing steepest descent directions whenever the steepest descent line minimum is lessthan the Newton line minimum). While 'newton' is the default, the user may wish totry the 'frcg', which can be less expensive computationally.argument descriptionrc {'real', 'complex' }mode {'frcg', 'newton' }metric {'flat', 'euclidean', 'canonical' }motion {'approximate', 'exact' }verbose {'verbose', 'quiet'}gradtol first convergence toleranceftol second convergence tolerancepartition partition of p indicating symmetries of f* metric The metric argument selects the kind of geometry with which to endow theconstraint surface. This ultimately affects the definition of the covariant Hessian.'flat' projects the result of applying the unconstrained Hessian onto the tangent space,while 'euclidean' and 'canonical' add on connection terms specific to their geometries.'euclidean' is the default.* motion The motion argument selects whether line movement along the manifold willbe by the analytic solution to the geodesic equations of motions for the metric selected,or by a computationally less expensive approximation to the solution (default). (For a'flat' metric, there is no geodesic equation, so this argument has no effect in that case.)* verbose The verbose argument determines whether the function will display reportson each iteration while the function executes. When this argument is 'verbose' (thedefault) data will be displayed and also recorded in the global SGdata. When thisargument is 'quiet' then no convergence data is displayed or recorded.* gradtol and ftol We declare convergence if either of two conditions are true:grad/gradinit < gradtol (default le-6) or (f-f old)/f < ftol (default le-8),where gradinit is the initial magnitude of the gradient and fold is the value of F(Y)at the last iteration.* partition The partition is a cell array whose elements are vectors of indices thatrepresent a disjoint partition of 1:p. If F has no symmetry, then the partition isnum2cell(l:p). If F(Y) = F(YQ) for all orthogonal Q, then the partition is {l:p}.The partition is {1:2,3:5} if the symmetry in F is F(Y) = F(YQ) for orthogonal Qwith sparsity structurexxxxxxxxxxX X XThe user could equally well specify {3: 5, 1: 2} or { [5 3 4] , [2 1] } for the partition,i.e. a partition is a set of sets. The purpose of the partition is to project away thenull-space of the Hessian resulting from any block diagonal orthogonal symmetries ofF(Y). If the argument is omitted, our code will pick a partition by determining thesymmetries of F (using its partition.m function). However, the user should be awarethat a wrong partition can destroy convergence.4.3 Sample problems and their differentialsThis section serves as a sample guide on the manipulation of objective functions of matriceswith orthonormal columns. We have found a few common tricks worth emphasizing.Once one has a formula for the objective function F(Y), we define the formula for dF(Y)implicitly by tr(V*dF(Y)) = _F(Y(t))t=o where Y(t) Y + tV (or any curve Y(t) forwhich Y(O) = V). The reader may recall that tr(A*B) = Ei, A3,B 3 , so it functions justlike the real inner product for vectors1 and the implicit definition of dF(Y) is actually thedirectional derivative interpretation of the gradient of F(Y) as an unconstrained function ina Euclidean space.For most of the functions we have used in our examples, the easiest way to obtain theformula for dF is to actually use the implicit definition.For example, if F(Y) = tr(AY*BY) one then hastr(V*dF(Y)) = tr(AV*BY + AY*BV).Since the value of the trace is invariant under cyclic permutations of products and transposes,we may rewrite this equation astr(V*dF(Y)) = tr(V*BYA + V*B*YA*),and, since V is unrestricted, this implies that dF(Y) = BYA + B*Y*A*.lover complex numbers, we always take the real part of tr(A*B).The process we recommend is:* try to write F(Y) as a trace* compute -F(Y(t))t=o where we let V = Y(t)* use trace identities to rewrite every term to have a V* in the front* strip off the V* leaving the dF(Y)As a check, we recommend using the finite difference dF.m code supplied in the subdi-rectory finitediff to check derivations before proceeding.The software needs a function called ddF.m which returns -4dF(Y(t)) t=o for Y(O) = H.The sort of second derivative information required by the software is easier to derive thanthe first. If one has an analytic expression for dF(Y), then one need only differentiate.If, for some reason, the computation for ddF.m costs much more than two evaluations ofdF(Y) with dF.m, the reader may just consider employing the finite difference function forddF.m found in finitediff (or simply use ddF.m as a check).It is, however, strongly suggested that one use an analytic expression for computingdF(Y), as the finite difference code for it requires a large number of function evaluations(2 np).4.3.1 The Procrustes problemThe Procrustes problem (see [32]) is the minimization of IAY - YBIIF for constant A andB over the manifold Y*Y = I. This minimization determines the nearest matrix A to A forwhichQ*AQ= )i.e. the columns of B span an invariant subspace of A.The differential of F(Y) = I IAY - YB 2 = ltr(AY - YB)*(AY - YB) is given bydF(Y) = A*(AY - YB) - (AY - YB)B*.This can be derived following the process outlined above. Observe thatd 1dF(Y(t)) t=o = -tr((AV - VB)*(AY - YB) + (AY - YB)*(AV - VB))dt 2= tr((AV - VB)*(AY - YB))= tr(V*(A*(AY - YB)) - V*(AY - YB)B*).The second derivative of F(Y) is given by the equation,ddF(Y(t)) t=o = A*(AH - HB) - (AH - HB)B*,where Y(O) = H, which can be obtained by varying the expression for dF.4.3.2 Nearest Jordan structureNow suppose that the B block in the Procrustes problem is allowed to vary with Y. Moreover,suppose that B(Y) is in the nearest staircase form to Y*AY, that is:B(Y) = 0 A2 *0 0 ...for fixed block sizes, where the *-elements are the corresponding matrix elements of Y*AYand the Ail blocks are either fixed or determined by some heuristic, e.g. taking the averagetrace of the blocks they replace in Y*AY. Then a minimization of I|AY - YB(Y)H F findsthe nearest matrix a a particular Jordan structure, where the structure is determined by theblock sizes, and the eigenvalues are Ai. When the Ai are fixed, we call this the orbit problem,and when the Ai are selected by the heuristic given we call this the bundle problem.Such a problem can be useful in regularizing the computation of Jordan structures ofmatrices with ill-conditioned eigenvalues.The form of the differential of F(Y), surprisingly, is the same as that of F(Y) for theprocrustes problem.dF(Y) = A*(AY - YB(Y)) - (AY - YB(Y))B(Y)*.This is because tr(-(AY - YB(Y))*YB) = tr((B(Y) - Y*AY)*B) = 0 for the B selectedas above for either orbit or bundle case, where B = $B(Y(t)) t=o.In contrast, the form of the second derivatives is a bit more complicated, since the Bnow depend on Y.d dddF(Y(t)) t=o = -dFproc(Y(t)) t=o - A*YB - AYB* + YBB* + YBB*,dt dtwhere Y(0) = H, dFproc is just short for the Procrustes (B constant) part of the secondderivative, and B = -4B(Y(t)) lt=o which is the staircase part of Y*AH + H*AY (with traceaverages or zeros on the diagonal depending on whether bundle or orbit).4.3.3 Trace minimizationIn this case we consider a "Rayleigh quotient" style of iteration to find the eigenspaces ofthe smallest eigenvalues of a symmetric positive definite matrix A. That is, we can minimizeF(Y) = ltr(Y*AY) [2, 36, 54]. The minimizer Y* will be an n x p matrix whose columnsspan the eigenspaces of the lowest p eigenvalues of A.For this problemdF(Y) = AY.It is easily seen thatddF(Y(t))Jt=o = AH,dtwhere Y (0) = H.4.3.4 Trace minimization with a nonlinear termWe consider a related nonlinear eigenvalue problem, to the trace minimication of section 3.3,minimize F(Y) = 1(tr(Y*AY) + g(Y)), where g(Y) is some nonlinear term.2\11~ 1/'Y1/,'"' \ UU"""""" ""This sort of minimization occurs in electronic structures computations such as LocalDensity Approximations (LDA) where the columns of Y represent electron wavefunctions, Ais a fixed Hamiltonian, and g(Y) represents the energy of the electron-electron interactions(see, for example, [4, 57]). In this case, the optimal Y represents the state of lowest energyof the system.The only additions to the differentials of the trace minimization are terms from g(Y)which vary from problem to problem. For our example, we have taken g(Y) = p Zi i forsome coupling constant K, where p, = j Yij12 (the charge density).In this casedF(Y) = dFtracemin + Kdiag(p)Y,andd d dddF(Y(t)) t=o = -dFtracemin(Y(t)) t=o + Kdiag(p)H + Kdiag(-p)Y,where Y(O)= H, and p t=0 = , YijH j.4.3.5 Simultaneous Schur decomposition problemConsider two matrices A and B which in the absence of error have the same Schur vectors,i.e. there is a Y E O(n) such that Y*AY and Y*BY are both block upper triangular. Nowsuppose that A and B are somewhat noisy, from measurement errors, or some other kind oflossy filtering. In that case the Y that upper triangularizes A might not upper triangularizeB as well. How does one find the best Y?This is a problem that was presented to us by W. Schilders [55] who phrased it as a leastsquares minimization of F(Y) = '(|Ilow(Y*AY) I + Ilow(Y*BY)| ), where low(M) is amask returning the block lower triangular part of M, where M is broken up into 2 x 2 blocks.For this problem the differential is a bit tricky and its derivation instructive,tr(V*dF(Y)) = tr(low(Y*AY)*low(Y*AV + V*AY) + low(Y*BY)*low(Y*BV + V*BY))tr(V*dF(Y)) = tr(low(Y*AY)*(Y*AV + V*AY) + low(Y*BY)*(Y*BV + V*BY))tr(V*dF(Y)) = tr(V*(AYlow(Y*AY)* + A*Ylow(Y*AY) + BYlow(Y*BY)* + B*Ylow(Y*BY)))dF(Y) = AYlow(Y*AY)* + A*Ylow(Y*AY) + BYlow(Y*BY)* + B*Ylow(Y*BY)where the second equation results from observing that tr(low(M)*low(N)) = tr(low(M)*N),and the third from properties of the trace.With second derivatives given by,d dF(Y(t))|t=o = AHlow(Y*AY)* + A*Hlow(Y*AY) + AYlow( d(Y*AY)d(Y*AY)+A*Ylow( d(YAY) + BHlow(Y*BY)* + B*Hlow(Y*BY)dt+BYlow( Y*BY)* + B*Ylow( d(YBY),dt dtwhere Y(O) = H, d(Y*AY) It=o = H*AY + Y*AH, and d(Y*BY) It=o = H*BY + Y*BH.4.3.6 INDSCALAnother problem like simultaneous Schur problem, involving sets of matrices with similarstructures, is the problem of finding the best set of eigenvectors for a set of symmetricmatrices given that the matrices are known to be simultaneously diagonalizable, but mayhave significant errors in their entries from noise or measurement error. Instances of thisproblem arise in the Psychometric literature, where it is called an INDSCAL problem [21].Phrased in terms of a minimization, one has a set of symmetric matrices Ai and wishesto find Y E O(n) that minimizes F(Y) = jEi off(Y*AiY)| 2, where off(M) returns theoff-diagonal entries of M.Similar to Schur problem (though with symmetric matrices and a somewhat differentmatrix component) we havedF(Y) = Z2YAioff(Y*AiY),andd dF(Y(t)) t=o = Z 2HAioff(Y*AiY), Z 2YAoff(H*AiY + Y*AH)dti iwhere Y(O)= H.4.4 Running our examplesWe have provided implementations for the sample problems described.We present the examples and their verbose outputs. We have included these outputs sothe reader can check that his/her copy of the package is executing properly. The user shouldbe running matlab (version 5 or compatable) in the following directory:>> !lsFline.mREADMEconnection.mdgrad.mdimension.mdtangent.mfinitediffgrad.minvdgrad.mip.mmove.mnosym.mpartition.msg_frcg.msg_min.msg_newton.mtangent.m(note, there may also be the additional subdirectories example/, docs/, or @cell/).Each example problem has a subdirectory in the examples subdirectory.>> !ls examplesindscal jordan idatoy procrustes simschurEach of these subdirectories contains an implementation for F(Y) (dF(Y) and ddF(Y, H)),a parameters function to set a global, FParameters, of fixed parameters of F(Y) (thisfunction must be called before any other), a guess function to generate an initial guess, andpossibly some auxiliary functions for computing F(Y), related quantities, or instances of aspecific problem.>> !is examples/idatoyF.m dF.mKinetic.m ddF.m>> !is examples/jordanBlock.m dBlock.mF.m dF.mguess.mparameters.mddF. mfrank .mguess.mparameters.mtraceminpost .mWe executed the examples in the two supported optimization modes ('frcg', and 'newton')and either 'flat' or 'euclidean'. For different instances of the same problem, different modesmight perform best, so the reader should not feel that a particular mode will exhibit superiorperformance in all instances.4.4.1 A sample Procrustes problemIn this example, we attempt to find a Y which minimizes I IAY-YBIlF for a pair of randomlyselected A (12 x 12) and B (4 x 4).>> !cp examples/procrustes/*.m>> randn('state',0);>> [A,B] = randprob;>> parameters(A,B);>> YO = guess;>> [fn,Yn] = sg_min(YO,'newton','euclidean');iter grad0 1.967151e+01invdgrad: Hessian not1 1.201377e+01invdgrad: Hessian not2 6.998543e+00invdgrad: Hessian not3 6.754473e+00invdgrad: Hessian not4 7.037159e+00invdgrad: Hessian not5 2.716002e+00invdgrad: Hessian not6 2.611537e+00F(Y)3.145130e+01positive definite,1.334642e+01flops step type4365 noneCG terminating early299347 newton steppositive definite, CG terminating early7.885178e+00 549026 newton steppositive definite, CG terminating early4.082520e+00 828756 newton steppositive definite, CG terminating early2.156004e+00 1143569 newton steppositive definite, CG terminating early1.264874e+00 1481879 steepest steppositive definite, CG terminating early7.828101e-01 1779148 newton step3.540531e-014.749149e-01 2339394 newton stepL -5a'a).II)6-10oflops x 106Figure 4-1: Procrustes Problem2.098621e-027.350907e-051.832068e-053.361052e-013.360631e-013.360631e-01285866434260593991138newton stepnewton stepnewton stepFigure 1 shows the convergence curve for this run.4.4.2 A sample Jordan block problemWe now look for a nearby matrix with a Jordan structure J4 (a) E J2 (3) E (don't care) tothe 12 x 12 Frank matrix, a matrix whose Jordan structure is known to be very difficult.>> !cp examples/jordan/*.m>> A = frank(12);>> evs = sort(eig(A)); eigvs = [mean(evs(1:4)) mean(evs(5:6))];>> parameters(A,eigvs, [4; 2],'bundle')>> YO = guess;>> [fn,Yn] = sgmin(YO,'newton','euclidean');iter grad F(Y) flops step typeflopsFigure 4-2: Jordan Problem0 1.775900e-02invdgrad: Hessian not1 1.797005e-03invdgrad: max iterati2 1.615236e-03invdgrad: max iterati3 2.100793e-044 7.793143e-055 3.826000e-076 1.766973e-089.122159e-06 15535 nonepositive definite, CG terminating early4.382300e-07 4539120 newton stepons reached inverting the hessian by CG1.320791e-07 9451244 newton stepons reached inverting the hessian by CG1.364935e-08 14602068 newton step3.098557e-09 19631237 newton step2.924257e-09 24466406 newton step2.924253e-09 26592864 newton stepFigure 2 shows the convergence curve for this run.4.4.3 A sample trace minimization problemIn this problem we minimize trace(Y*AY) where A is a 12 x 12 second difference operatorin 1-dimension and Y is 12 x 4.- 1 5 L0 050 0.5 2.5x 10Note: we do not recommend that this problem be solved this way, since, for constant A,it can be done faster with conventional eigenvalue solvers.>> !cp examples/tracemin/*.m>> A = Kinetic(12);>> parameters(A);>> YO = guess(4);>> [fn,Yn] = sg_min(YO,'frcg','euclidean');grad2.236291e+001 .413602e+009.159910e-015.663336e-014.366693e-012.676809e-011.943933e-011.200285e-015.106617e-022.633155e-021.009583e-024.054257e-031 .416079e-037.175948e-044.258507e-043.065818e-041 .260472e-043.596692e-05F(Y)3.833127e+001 .903903e+001.307258e+001.057510e+009.440163e-018.724658e-018.462110e-018.319293e-018.279622e-018.272484e-018.270559e-018.270313e-018.270275e-018.270269e-018.270268e-018.270267e-018.270267e-018.270267e-01Figure 3 shows the convergence curve for this run.100iter01234567891011121314151617flops33572222953502644903986221777634188993911054873120556613704301547221173751319403192154843242807027008162929162317079505 1 1.5 2 2.5 3flops x 106Figure 4-3: Trace Minimization Problem4.4.4 A sample LDA toy problem (trace minimization with non-linear term)We attempt to solve a toy Local Density Approximation problem described in section 3, forfour electrons on a 1D grid with twelves points and a coupling constant of 20. In this case,we use the first four eigenvectors of the linear problem as a starting point.>> !cp examples/ldatoy/*.m>> K = Kinetic(12);>> parameters(K,20);>> YO = guess(4);>> [fn,Yn] = sg_min(YO,'newton','euclidean');iter grad F(Y) flops step type0 2.432521e+00 7.750104e+00 4104 noneinvdgrad: Hessian not positive definite, CG terminating early1 4.769809e-01 7.546426e+00 303003 steepest step2 4.665924e-02 7.531450e+00 725498 newton step101flopsFigure 4-4: LDA Toy Problem3.363281e-041.678208e-067.531232e+007.531232e+0012437071868470newton stepnewton stepFigure 4 shows the convergence curve for this run.4.4.5 A sample simultaneous Schur decomposition problemIn this case we use the sample problem given to us by Schilders.>> !cp examples/simschur/*.m>> [A,B] = noisy;>> parameters(A,B);>> YO = guess;>> [fn,Yn] = sg_min(YO,'newton','flat');iter grad F (Y) flops2.205361e-016.688045e-031.386562e-041.687421e-024.928702e-034.913259e-035476348084458560957nonenewton stepnewton step10222 '0 2 10 12 14x 10flopsFigure 4-5: Simultaneous Schur Problem2.377737e-062.314313e-074.913255e-034.913255e-031189102514846552newton stepnewton stepFigure 5 shows the convergence curve for this run.4.4.6 A sample INDSCAL problem (simultaneous diagonaliza-tion)We attempt to simultaneously diagonalize two noisy versions of diag (1: 10) and diag (1: 10) ^2.Our initial point is the QR-factorization of the average of the eigenvectors of the two matri-ces.>> !cp examples/indscal/*.m>> randn('state',0);>> [A,B] = noisy;>> parameters(A,B);>> YO = guess;>> [fn,Yn] = sg_min(YO,'newton','euclidean');10312x 106iter grad0 1.871820e+03invdgrad: Hessian not1 1.141235e+03invdgrad: Hessian not2 6.017399e+02invdgrad: Hessian not3 3.157921e+02invdgrad: Hessian not4 1.908635e+02invdgrad: Hessian not5 1.385145e+02invdgrad: Hessian not6 1.274498e+027 2.162910e+008 3.022615e-039 1.708804e-03F(Y) flops5.837012e+02 29736 nonepositive definite, CG terminating early1.996926e+02positive definite7.364112e+01positive definite2.351357e+01positive definite1.450185e+01positive definite1.136202e+01positive definite9.160689e+005.662713e-032.224608e-062.223022e-062514846 steepest stepa, CG terminating early8191352 newton stepi, CG terminating early11094356 newton step, CG terminating early15838109 steepest st, CG terminating early20790540 steepest st, CG terminating early25701418 newton step34638359 newton step42237728 newton step43731038 newton stepepepFigure 6 shows the convergence curve for this run.4.5 Making modifications (the structure of the code)This template is written to function as is, but can be tailored to specific problems, resultingin substantial improvements in performance. There are also functions which could be in-lined and computations which can be reused of which we have not taken advantage for thesake of having a more readable code. Many components are the usual routines found inthe literature (our line minimization routine is matlab's fmin for example). Certainly itis possible for improvements to be made in the optimization routines themselves and wewelcome suggestions from optimization experts.To make modifications as painless as possible, we present a section sketching the basic104o0-5--10 i I I i I0 0.5 1 1.5 2 25 3 35 4 4.5flops x 107Figure 4-6: (INDSCAL) Simultaneous Diagonalization Problemstructure of the sgmin code. Readers who do not wish to modify any of the functionsbeyond F.m, dF.m, and ddF.m can skip this section.4.5.1 Subroutine dependenciesFigure 7 shows the dependencies of the various subroutines on each other. We have groupedtogether the routines based on four loosely defined roles:* Objective Function These are the routines which implement the objective functionand its derivatives. This group is the only group that a user interesting only in thebasic functionality need ever modify to adapt the template. Functions in this groupare used by functions in the Geometrized Objective Function group and High LevelAlgorithm group.* Geometric Implementation These routines implement the geometric features of theStiefel manifold. This includes the projection of unconstrained vectors onto the con-straint (i.e. tangent) space (tangent), line movement (move), and the connection termused for covariant differentiation (connection). These routines are independent of the105Figure 4-7: Our code tree: dependencies of various modules106HIGH LEVEL ALGORITHMsg_minsgnew on sg_frcgFline.mGEOMETRIZED OBJECTIVE FUNCTIONinvdgrad.m grad.mdgrad.m dimension.mnosym.mpartition.mOBJECTIVE FUNCTIONF.mdF.mddF.mGEOMETRIC IMPLEMENTATIONmove.mconnection.mtangent.m1p.mdtangent.mObjective Function group, but are essential to the Geometrized Objective Functiongroup.* Geometrized Objective Function This group uses the routines of the GeometricImplementation to project the differential and second differential of F(Y) onto theconstraint (i.e. tangent) surface producing the geometrically covariant gradient (grad)and Hessian (dgrad). Additionally, the routine for inverting the covariant Hessianby conjugate gradient (invdgrad) is located here. This group provides the raw toolsout of which one builds implementations in the High Level Algorithm group. Lastly,functions which detect (partition) and remove (nosym) block diagonal orthogonalsymmetries are found in the group.* High Level Algorithm This group implements the constrained versions of variousoptimization algorithms. It is here search directions are selected, line minimizationsare performed (using Fline), and convergence criteria are defined.Lastly, every function reads from a global SGParameters structure whose fields containinformation about the manifold and the computation. In our examples we have used aseparate global FParameters to store information related to the computation of F(Y). Thefields of SGParameters are set in sg_min.4.5.2 Under the hoodsgmin is built from a couple fairly elementary optimization codes which have been putthrough a geometrization that allows them to sit on the Stiefel manifold. The basic elementsof a geometrization are the rules for how to take inner products, how to turn unconstraineddifferentials into constrained gradients, how to differentiate gradient fields covariantly, andhow to move about on the manifold.The sgmin routine parses the arguments and sets the defaults. Finally, it calls eithersgnewton or sgfrcg. The source for these routines (in pseudocode) follows:107function [fn,Yn]= sgnewton(Y)ginitial = grad(Y); g = ginitiai;f = F(Y); fold = 2*f;while (llgll > eps) & (lglI/lgintiaill >gradtol) & (Ifold/f-1 I > ftol)sdir = -g;% steepest descent direction and line minimization on different scalesfsa = minimize F(move(Y,sdir,sa)) for sa E [-1, 1] I i-fsb = minimize F(move(Y,sdir,sb)) for sb E [-1, 1] I (gsdisdr4A -rdndir = -dgrad - 1 - g;% newton direction and line minimizationfna = minimize F(move(Y,ndir,na)) for na E [-1, 1]fnb = minimize F(move(Y,ndir,nb)) for nb E [-1, 1]r)Si f(g,ndir)" I (ndir,dgrad.ndir) I;% compare the best newton function value with the best steepestif (fsa < fsb) st=sa; fst=fsa; else st=sb; fst=fsb; endif (fna < fnb) nt=na; fnt=fna; else nt=nb; fnt=fnb; endif (fst eps) & (lglll/lginitialll >gradtol) & (Ifold/f-1 I > ftol)dir = -g;% get a Hessian conjugate directionif (not first iteration)dir = dir - dir.(draddir, ) *dirod;dirld.(dgrad-dirld)end% line minimizationfcga = minimize F(move(Y,dir,cga)) for cga E [-1, 1] - (g,dir)fcgb = minimize F(move(Y,dir,cgb)) for cgb E [-1, 1]. - (dir,d ir) ;if (fcga < fcgb) t=cga; fnew=fcga; else t=cgb; fnew=fcgb; end% move to the new point[Y,dirold] = move(Y,dir,t); fold= f; f fne,,,;g=grad(Y);endfn =f;Yn = Y;-- -1-11--M-NN IIIII 01 MINHII iIlr1r 10 I i I IY1 1 Ifunction to employ for line minimization. If one wished to use a function other thanfmin then one could modify sgnewton.m and sg_frcg.m appropriately." High Level Algorithms The user may wish to use a minimization routine differentfrom or more sophisticated than sgfrcg or sg_newton, or one may just have a differentstopping criteria [6] that we did not anticipate. It is possible to adapt flat optimizationalgorithms to new sg algorithms without having to know too many details. Mostunconstrained minimizations contain line searches (y -+ y+t*v), Hessian applictions(H*v) and inner products (v'w). To properly geometrize the algorithm, one has to do nomore than replace these components with move (Y, V, t), dgrad (Y, V), and ip (Y, V, W)respectively.* Hessian Inversion We currently provide codes to invert the covariant Hessian in thefunction invdgrad (which uses a conjugate gradient method). Any matrix inversionroutine which can stably invert black-box linear operators can be used in place of theseand no other modifications to the code would be required. Another feature the usermay wish to add is an inversion routine to produces dog-leg steps.* Workspaces To make the code more readable we have not implemented any workspacevariables even though there are many computations which could be reused. For ex-ample, in the objective function group, for many of the sample problems, products ofthe form A*Y or Y*B or complicated matrix functions of Y such as B (Y) in the Jordanexample could be saved in a workspace. Some of the terms in the computation of gradand dgrad could likewise be saved (in fact, our current implementation essentiallyrecomputes the gradient every time the covariant Hessian is applied).* In-lining Although we have separated the tangent and dtangent functions from thegrad and dgrad functions, one often finds that when these functions are in-lined, thatsome terms cancel out and therefore do not need to be computed. One might alsoconsider in-lining the move function when performing line minimizations so that onecould reuse data to make function evaluations and point updates faster.110* Sphere or O(n) Geometry Some problems are always set on the unit sphere or onO(n), both of which have simpler geometries than the Stiefel manifold. Geodesics onboth are easier to compute, as are tangent projections. To take advantage of this, oneshould modify the Geometric Implementation group.* Globals In order to have adaptable and readable code, we chose to use a globalstructures, SGParameters and FParameters, to hold data which could be put in theargument lists of all the functions. The user may wish to modify his/her code so thatthis data is passed explicitly to each function as individual arguments.* Preconditioning The Hessian solve by conjugate gradient routine, invdgrad, does nothave any preconditioning step. Conjugate gradient can perform very poorly withoutpreconditioning, and the preconditioning of black box linear operators like dgrad is avery difficult problem. The user might try to find a good preconditioner for dgrad,but this may be difficult.4.6 Geometric technicalitiesThis section is a brief explanation of the geometric concepts used to produce the methodsemployed by sg_min and its subroutines.4.6.1 ManifoldsA manifold, M, is a collection of points which have a differential structure. In plain terms,this means that one is able to take derivatives of some reasonable class of functions, theC" functions, defined on M. What this class of differentiable functions is can be somewhatarbitrary, though some technical consistency conditions must be satisfied. We will not bediscussing those conditions in this section.We consider the manifold StiefR(n, k) (for purposes of clearer explanation we restrict ourdiscussion to the real version of the Stiefel manifold) of points which are written as n x kmatrices satisfying the constraint Y*Y = I. We will select for our set of C' functions thosereal valued functions which are restrictions to Stief(n, k) of functions of nk variables which111are C"(Rnxk) about Stief(n, k) in the Rnxk sense. It should not be difficult to convinceoneself that the set of such functions must satisfy any technical consistency condition onecould hope for.Additionally, M = Stief(n, k) is a smooth manifold. This means that about every pointof Stief(n, k) we can construct a local coordinate system which is C'. For example, considerthe pointY=and a point in the neighborhood of Y,A)For small B, we can solve for A in terms of the components of B and the components ofan arbitrary small antisymmetric matrix A by solving S*S = I - B*B for symmetric S andletting A = eAS. This can always be done smoothly and in a locally 1-to-1 fashion for smallenough B and A. Since the components of B are all smooth functions (being restrictionsof the coordinate functions of the ambient Rnxk space), and since the solution for A isC"(R(n- k)xk D Rkxk) for small enough B and A, we have shown that any point of Stief(n, k)can be expressed smoothly as a function of (n - k)k + k(k - 1)/2 = nk - k(k + 1)/2 variables.The only difference between is a euclidean rigid motion; therefore, this statement holds(0for all points in Stief(n, k).Summarizing, a manifold is a set of points, with a set of differentiable functions definedon it as well as a sense of local coordinatization in which the dependent coordinates canbe represented as smooth functions of the independent coordinates. Specifically, we seethat there are always nk - k(k + 1)/2 independent coordinates required to coordinatizeneighborhoods of points of Stief(n, k). This number is called the dimension of the manifold,and it should be the same for every point of a manifold (if not, then the manifold is eitherdisconnected or not actually smooth).1124.6.2 The difference between a tangent and a differentialGiven a smooth point x on any manifold and a smooth path p(s) in the manifold such thatp(O) = x, one can construct a differential operator on the C" functions defined near x bythe ruleDpf = df(p(s)) S=o0dsNow, it turns out that these differential operators at x form a finite dimensional vector space,and the dimension of this vector space is equal to the dimension of the manifold. This vectorspace is called the tangent space of M at x, and is often denoted Tx(M).For the Stiefel manifold, and, generally, all smooth constraint manifolds, the tangentspace at any point is easy to characterize. The constraint equation Y*Y = I can be thoughtof as k(k + 1)/2 independent functions on Stief(n, k) which must all have constant value. IfH is a tangent vector of Stief(n, k) at Y, it must then be thatH*Y + Y*H = 0,which is obtained by taking dY(t)*Y(t) = 0 t=o (where Y(0O) = H). In a constraint manifold,the tangents at Y can equivalently be identified with those infinitesimal displacements of Ywhich preserve the constraint equations to first order.For H to be a tangent, its components must satisfy k(k + 1)/2 independent equations.These equations are all linear in the components of H, and thus the set of all tangents is ank - k(k + 1)/2 dimensional subspace of the vector space Rnxk.A differential operator DH may be associated with an n x k matrix, H, given bydDHf = f (Y + tH) t=o.dtDH is the directional derivative operator in the H direction. One may observe that for Hto be a tangent vector, the tangency condition is equivalent to DH(Y'Y) = 0.While tangents are n x k matrices and, therefore, have associated differential operators,differentials are something else. Given a C o function f and a point Y, one can consider the113equation DHf for H E Rnxk (H is not necessarily tangent). This expression is linear in Hand takes on some real value. It is thus possible to represent it as a linear function on thevector space Rnxk,DHf = tr(H*Z),for some appropriate n x k matrix Z whose values depend on the first order behavior of fnear Y. We identify this Z matrix with df, the differential of f at Y. This is the samedifferential which is computed by the dF functions in the sample problems.For any constraint manifold the differential of a smooth function f can be computedwithout having to know anything about the manifold itself. One can simply use the differ-entials as computed in the ambient space (the unconstrained R n xk derivatives in our case).If one then restricts one's differential operators to be only those in tangent directions, thenone can still use the unconstrained df in tr(H*df) to compute DHf for H E Ty(Stief(n, k)).This is why it requires no geometric knowledge to produce the dF functions.4.6.3 Inner products, gradients, and differentialsIn flat spaces, we often identify the differential of a function with its gradient. However,when dealing with a more general setting, one can run into problems making sense out ofsuch a definition.For example, the gradient is a vector, and it should be possible to think of vectors asinfinitesimal displacements of points. In Stief(n, k), any infinitesimal displacement 6Y mustsatisfy 6Y*Y + Y*6Y = 0. Thus, df may not always be a vector, since it does not necessarilysatisfy this equation. A gradient should be an infinitesimal displacement that points in thedirection of the displacement which will give the greatest increase in f.If the tangent space has an inner product, though, one can find a useful way to identifythe df uniquely with a tangent vector. Let ip(Hi, H2) be a symmetric nondegenerate bilinearform on the tangent space of Stief(n, k) at Y. Then one can define the gradient, G, implicitlyby,ip(G, H) = tr(H*df) = DHf.114dFtangent spaceilmii i Mimit m IIFigure 4-8: The unconstrained differential of F(Y) can be projected to the tangent space toobtain the covariant gradient, G, of F.Since ip is a nondegenerate form, this is sufficient to define G. The function tangent carriesout this projection from differentials to tangents (shown in Figure 4-8). This operation isperformed by grad to produce the gradient of the objective function.4.6.4 Getting around Stief(n, k)We've now laid the groundwork for a meaningful definition of the gradient in the setting of aconstraint manifold. At this point, one could run off and try to do a steepest descent searchto maximize one's objective functions. Trouble will arise, however, when one discovers thatthere is no sensible way to combine a point and a displacement to produce a new pointbecause, for finite s, Y + sH violates the constraint equations, and thus does not give apoint on the manifold.For any manifold, one updates Y by solving a set of differential equations of motion ofthe formd Y=Hdtd H = -F(H,H).dtThe F term is crafted to ensure that H remains a tangent vector for all times t, thus keepingthe path on the manifold. It is called the connection. (The connection term also depends onY.)To see how these equations could be satisfied, we take the infinitesimal constraint equation115IimilmmimIem iMimi iMfor H,H*Y + Y*H = 0,and differentiate it with respect to t, to getF(H, H)*Y + Y*F(H, H) = 2H*H.This can be satisfied generally by F(H, H) = Y(H*H) + T(H, H) where T(H, H) is somearbitrary tangent vector.In the next subsection, we will have reason to consider F(HI, H 2) for H 1 : H 2. Fortechnical reasons, one usually requests that F(H 1, H2) = F(H2 , H 1) (this is called the torsionfree property), and that ip(F(H1, H3 ), H2) + ip(Hi, F(H 2, H3 )) = 0 (this is called the metriccompatability property). These two properties uniquely determine the T(H1, H 2) term, andthereby uniquely specify the connection. On any manifold, the unique connection with theseproperties is called the Levi-Civita connection.The connection for the Stiefel manifold using two different inner products (the Euclideanand the canonical) can be found in the work by Edelman, Arias, and Smith (see [1, 31, 56]).The function connection computes F(H 1, H2 ) in the template software.Usually the solution of the equations of motions on a manifold are very difficult to carryout. For the Stiefel manifold, analytic solutions exist and can be found in the aforementionedliterature, though we have found very little performance degradation between moving alongpaths via the equations of motion and simply performing a QR factorization on Y + sH,as long as the displacements are small. Since QR factorization is cheaper, the softwarehas this as its default (motion= 'approximate'). We have provided the two equations ofmotion solvers (for the Euclidean and canonical connections) as well as movement by QRfactorization in move to satisfy both the purists and the pragmatists respectively.4.6.5 Covariant differentiationWith the notion of a gradient and a notion of movement that respects the constraints of themanifold, one might wish to begin with an optimization routine of some sort. A steepest116Figure 4-9: In a flat space, comparing vectors at nearby points is not problematic.Figure 4-10: In a curved manifold, comparing vectors at nearby points can result in vectorswhich do not lie in the tangent space.descent could be implemented from these alone. However, in order to carry out sophisti-cated optimizations, one usually wants some sort of second derivative information about thefunction.In particular, one might wish to know by how much one can expect the gradient to changeif one moves from Y to Y + EH. This can actually be a difficult question to answer on amanifold. Technically, the gradient at Y is a member of Ty(Stief(n, k)), while the gradientat Y + EH is a member of Ty+eH(Stief(n, k)). While taking their difference would work finein a flat space (see Figure 4-9), if this were done on a curved space, it could give a vectorwhich is not a member of the tangent space of either point. (see Figure 4-10).A more sophisticated means of taking this difference is to first move the gradient atY + EH to Y in some manner which translates it in a parallel fashion from Y + EH to Y,and then compare the two gradients within the same tangent space. One can check that forV E Ty+EH(Stief(n, k)) the ruleV -4 V + eFr(V, H),117-- -------- - -'- -""IIIC ------~ ~II - IIIYliilliiwhere F(V, H) is the Levi-Civita connection, takes V to an element of Ty(Stief(n, k)) andpreserves inner product information (to first order in c). This is the standard rule for paralleltransport which can found in the usual literature ([14, 69, 44, 41], and others).Using this rule to compare nearby vectors to each other, one then has the following rulefor taking derivatives of vector fields:dDHG = G(Y + sH)s=0o + r(G, H),dswhere G is any vector field (but we are only interested in derivatives of the gradient field).This is the function implemented by dgrad in the software.In an unconstrained minimization, the second derivative of the gradient g = Vf along avector h is the Hessian [ f ] times h. Covariantly, we then have the analogy,&2[2f - ]h = ( -. V)g - DHG.4.6.6 Inverting the covariant Hessian (technical considerations)Since the tangent space of Stief(n, k) is a subspace of Rnxk, the covariant Hessian must beinverted stably on this subspace. This requires any algorithms designed to solve DHG = Vfor H to really be pseudoinverters in a least squares or some other sense.A second consideration is that many useful functions f(Y) have the property thatf(YQ) = f(Y) for all block diagonal orthogonal Q (i.e.Q 0 ... 0Q= 0 Q2 ... 00 0 ... Qpwhere the Qi are orthogonal matrices). In this case, all tangent vectors of the form(A, 0 ... 0H=Y 0 A 2 -.. 00 0 --- AP118where the Hi are antisymmetric, have the property DHf = 0. These extra symmetry vectorsare then null vectors of the linear system DHG = V. Because of these null directions, theeffective dimension of the tangent space is reduced by the dimension of this null space (thisis the dimension returned by the dimension function).Thus, we see that in order to invert the covariant Hessian, we must take care to usea stable inversion scheme which will project out components of H which do not satify theinfinitesimal constraint equation and those which are in the direction of the additional sym-metries of f. The invdgrad function carries out a stable inversion of the covariant Hessian bya conjugate gradient routine, with the dgrad function calling the function nosym to projectout any extra symmetry components.119------- - ~~~ ~----~'Y ~ imilmIikW1,1161' ii rYIIIY I J15'' IIYII I 1II I lrll120Chapter 5Multiscale Computation withInterpolating Wavelets5.1 IntroductionWavelets offer a means of approximating functions that allows selective refinement. If regionsof an image or a signal have exceptionally large variations, one need only store a set ofcoefficients, determined by function values in neighborhoods of those regions, in order toreproduce these variations accurately. In this way, one can have approximations of functionsin terms of a basis that has spatially varying resolution. This approach reduces the memorystorage required to represent functions and may be used for data compression.Physical applications often involve multiple physical fields which interact in space withnon-linear couplings. Efficient implementations must minimize not only storage to representthe these fields but also the processing required to describe their interactions. It is highlydesirable to perform the needed operations with a fixed, limited number of floating pointoperations for each expansion coefficient and to minimize the number of points in space atwhich physical interactions must be evaluated.As a concrete example of a realistic application, consider the computation of the quantummechanical electronic structure of a collection of atoms in three dimensions. For otherexamples of physical applications, the reader may wish to consult [7], [9], [35], et al.. Arias121-- -rr~- ^- --- l --r~- -r--r-irrr~'l"lYirr~I-^.-~,~yl~~------ ---^ IP-l-lll - - __ YI.I~. ..IIYI~- 1 C-~ I aand coworkers [15],[5], and other works which have appeared after the original submission ofthis manuscript nearly one year ago [65],[64], have studied the use of multiresolution bases inquantum mechanical computations. (For a review, see [3].) It is a consequence of quantumphysics that, near atomic nuclei, electronic wave functions vary rapidly and, that far fromthe nuclei, the wave functions tend to be much more smooth. From this observation, one cananticipate that the fine scale coefficients in a multiresolution analysis of the electronic wavefunctions and associated physical fields will be significant only near the atomic cores, allowingfor truncation. This problem is thus an ideal candidate for multiresolution techniques.Within density functional theory [42], the quantum physics of electrons and nuclei in-volves two types of fields, the Schr6dinger wave function {4'(r)} for each electron i and theelectrostatic potential O(r) arising from the average electronic density n(r) - Ei i (r) 2.Within the local density approximation (LDA) [45], the solution for the correct values ofthese fields is obtained at the saddle-point of lowest energy of the Lagrangian functionalL'LDA({4ij}, ) = f3 r I3VV (r)1 2 + d3r Vnuc(r)n(r) + d3 r c,(n(r))n(r)- Jd3r (r)n(r) - Jd3r IV(r)1 2 . (5.1)Here, we work in units h = m = e = 1, Vnuc(r) is the total potential which each electronfeels due to the presence of the atomic nuclei, and e,,(n) is a non-algebraic function knownonly through tabulated values. (For a review of density functional theory, see [51].)In practice, one finds the fields {4(r)}, O(r) by" expanding the fields in terms of unknown coefficients within some basis setO' (x) = Zcib,(x) (5.2)a(x)= db,);* evaluating Equation (5.1) in terms of the unknown coefficients c and d;* determining analytically the gradients of the resulting LLDA (c, d) with respect to thosecoefficients; and122* proceeding with conjugate gradients to locate the saddle point.All follows directly once one has expressed the Lagrangian as a function of the expansioncoefficients.In doing this, we note that each term represents a local coupling in space, but that onecoupling, 0(r)n(r), is cubic in the field coefficients c and d, and another, 'xc(n(r))n(r), isonly known in terms of tabulated values. Expanding the product of two wavelets in terms ofwavelets on finer levels would make possible the treatment of the cubic coupling to some levelof approximation. (See, for example, [10].) However, this route becomes increasingly difficultfor higher order interactions and is hopeless for non-algebraic or tabulated interactions, suchas exc(n(r)). For higher order interactions it is natural, and for non-algebraic and tabulatedinteractions necessary, to evaluate the interactions at some set of points in space and thenrecover expansion coefficients for the result. One then relies upon the basis set to provideinterpolation for the behavior at non-sample points.The benefits of both truncated wavelet bases and interpolation on dyadically refined gridsare given by the use of interpolating scaling functions [27], [11], [12], [24], [26] (or interpolets[70], [3]), which are functions with the following properties (from [26], pp. 6-7).Let O(x) be an interpolet, then(INT1) cardinality: 0(k) = 6 0,k for all k E Zn(INT2) two-scale relation: 0(x/2) = ZEycz cy,(x - y)(INT3) polynomial span: For some integer m > 0, any polynomial p(x) of degree mcan be represented as a formal sum EYC, a(y)o(x - y).Cardinality allows the fast conversion between uniform samples and interpolating scalingfunctions and has subtle yet profound consequences for the resulting multiresolution basis.In particular, as is evident from our algorithms below, the expansion coefficient for a ba-sis function on a particular scale is independent of the samples of the function for pointsassociated with finer scales. Consequently, the expansion coefficients which we obtain forfunctions maintained in our basis are identical to what would be obtained were function sam-ples available on a complete grid of arbitrarily fine resolution. This eliminates all error in123the evaluation of non-linear, non-algebraic and tabulated interactions beyond the expansionof the result in terms of a finite set of basis functions.The cy in the two-scale relation are referred to as scaling coefficients, and cardinalityactually implies that c, = 0(y/2). The two-scale relation allows the resolution to varylocally in a mathematically consistent manner.The polynomial span condition captures, in some sense, the accuracy of the approxi-mation. By cardinality, we actually have a(y) = p(y). We shall call m the polynomialorder.Interpolets thus can be thought of as a bridge between computations with samples ondyadically refined grids and computations in a multiresolution analysis. The former pointof view is useful for performing local nonlinear operations, while the latter is useful for theapplication of local linear operators.This manuscript explores O(N) algorithms that calculate transforms and linear operatorsfor grids of variable resolution but return, for the coefficients considered, exactly the sameresults as would be obtained using a full, uniform grid at the highest resolution without theneed to introduce artificial temporary augmentation points to the grid during processing.We thus show that with relatively mild conditions on the variability of the resolution pro-vided by the grid, interpolet bases provided the ultimate economy in the introduction of gridpoints: only as many samples in space need be considered as functions used in the basis. Thefour transforms (forward, inverse, and the dual to each) mapping between coefficients andfunctions samples which we discuss here are particular to interpolet bases. For the applica-tion of operators in such bases, we show that the familiar non-standard multiply of Beylkin,Coifman and Rokhlin[8] shares with the transforms the property of correctness without theneed to introduce additional grid points. Furthermore, we weaken the condition on grid vari-ability by using a modification of the non-standard multiply. We generalize the non-standardmultiply so that communication may proceed between nearby but non-adjacent levels andthereby obtain less stringent conditions on the variability of the grid. All of theoretical re-sults in this manuscript are presented in a general d-dimensional space. Illustrative examplesfor the purpose of discussion will be given in d = 1 and d = 2 dimensions. The examplesof applications in the final section will be in d = 3 dimensions. Our focus is entirely on124Ck 2k - Zn for k > O, Z n for k < 0.Dk Ck-1 -Ck for k > 0, for k < 0.Ok min(k, m) where m is the largest integer such that 2m divides all thecomponents of y.Fk (S) the space of functions over S C Z n .Ik( , S) the space of linear combinations of q((x - y)/2 0 k(Y)) for y E S c Zn.L¢ the mapping from S - ZIk(, S) which takes y -+((x - y)/2 0 k(Y)) andlinearly extended to a map from Fk(S) - Ik(0, S).SFkl(Z")/Fk(Zn - S).4(() Ik(0, Zn)/I_(0, Zn _ S).'i the zero-lift representative of v E FS , i.e. ' E Fk, such that b(y) = v(y),y S and i(y) = O, y V S.V* the dual space of the vector space V.Jk1 ,k 2 the map, k, - k 2, given by Jkl,k2 = O Lk1Jk short for Jo,k.J-k short for Jk,0.(f101g) the matrix element f f(x)Og(x)d"x.Table 5.1: Notation for multiscale computations.interpolet bases, and so it remains an open question whether these results hold true or canbe adapted to more general wavelet systems.Our organization is as follows. In Sections 2 and 3, we explain how to manipulate andconstruct interpolet expansions and some aspects of how well they perform. These sectionswill present nothing new to the experienced wavelet user, but will explain our notationalconventions (which are summarized in Table 5.1) and recapitulate common theorems ([19],[53], [20], [58], et al.) for wavelet novices. In Section 4, we describe how nonuniform basescan be conceptualized in the framework of interpolet expansions and then use our resultsto develop algorithms for the transforms. Section 5 details the algorithm for V 2 and otheroperators. Section 6 gives some practical details for the reader interested in implement-ing these algorithms. Finally, Section 7 compares, in three dimensions, timings of theseimplementations with the timings of naive algorithms. This final section also explores theconvergence of a preconditioned conjugate gradients algorithm in the solution of Poisson'sequation for the full three dimensional electrostatic potential arising from the nuclei in thenitrogen molecule.125~ ..~_~....~~.^~. U131Clx..- .LIIL~~I~LL~-- -- C .-~l~--XI~D1^1^-~-- ll1 -~-II C~~~~~~~~--~-~l~--~~IIXIYL_~~__ ^I -- -C--~ ---~- 1_ ~111_ ~QI~-5.2 Introduction to interpoletsThere is a unique set of interpolets on R having symmetry and minimal support for agiven polynomial order m = 21 - 1 (the Deslauriers-Dubuc functions [24]). These are thefunctions with which this article is primarily concerned (our results carry over to more generalinterpolets, and no use will actually be made of minimal support or symmetry).To determine the cy's, one sets c2j = 6m0 and c, = c-_. One may solve the Vandermondesystem,S1 1 .. 1 / c1 32 --- (21- 1)2 c31 34 -.. (21- 1)41 321-2 ... (21 - 1)21-2 C2 1- 1to obtain the remaining cy's. These coefficients satisfy the21- 1.The scaling coefficients for m = 1 areC_1 = C1 = 0.5, co = 1,and for m = 3 (the example used for Figure 5-1.) they are) 0)conditions for polynomial order3 = 16 _1 = C1 = 16, co 1.One may then take tensor products of q's and cy's to form interpolets in higher dimen-sions.5.2.1 Interpolet Multiresolution AnalysisWe are concerned with recursive representations of functions from samples at integer pointson both uniform and refined grids. There are many definitions which make the expositionmore clear.126a10.50-8 -7-6-5-4-3-2-1 0 1 2 3 4 5 6 7 8b10. 5 '8 -7-6-5-4-3-2-1 0 1 2 3 4 5 6 7 8Figure 5-1: (a) shows q(x/2) = .q5(x - 3) + '2( - 1) + O(x) + -2-(x + 1) + =10(x + 3).(b) shows the functions 1( - 3), -q(x - 1), O(x), 'q(x + 1), 0(x + 3).Definition 5.2.1 For k > 0, let Ck = 2kZ n , and let Dk = Ck-1 - Ck. For k < 0, letCk = Zn, and Dk = .We consider Ck to be the set of coarse lattice points on the lattice 2k-lZn and Dk, thedetail lattice points, to be those points on 2k-lZn which are not coarse. Note: Dk U CkCk-1, and Z = Ck U Dk U Dk-1 U ... U D 1 is a partition of Z.Definition 5.2.2 We let Ok(y) = min(k, m) where m is the largest integer such that 2mdivides all of the components of y. We call Ok(y) the level of the point y.Given the partition Zn = Ck U Dk U Dk-1 U ... U D 1, we haveOk (Y) k, yE C k1-1, y EDi.Definition 5.2.3 Let S C Z n . Let O(x) be an interpolet. Let Ik(o,S) be the space offunctions given by formal sums of the form Eyes a(y)O(- ).2 k(Y)127...,.~I-~---- ---~-----*------ -- ~- --- -^ --U-------- ' ~U~-~-- -- - --- ----------~I"-- I _aIr*-;t- U~- ~- -- ^Where q and S are understood, we may simply write k.Definition 5.2.4 Let S C Zn . Let Fk (S) be the vector space of R- or C- valued functionson Z n which are zero at any point not in S (i.e. with support contained in S).Where S is understood, we may simply write Fk. Note: TYk(S) = Fk(S n Ck) ( .Fk(S nDk) Fk(S n Dk-1) ED Fk(S n D).The meaning of the k subscript will be established by the next definition, which will linkvectors in Fk with functions in 1k. It is for this reason that while the F's are technicallyidentical, they are semantically different. In practice, the Y's are the actual data structuresbeing stored on the computer.Definition 5.2.5 Let O(x) be an interpolet. Let to : S - ZIk(/, S) be defined byx- y2k (Y)This definition extends linearly to the mapping to : Fk(S) - Ik(, S) defined by:yeSZ. e.(LkV)(X) = E v(y)q((x- y)/2k)± E v(y)((x - y)/2k-1).+ + v(y)(x- y).yEsnCk yESnDk yESnD1The set S can be thought of as the set of points in a refined grid. The to identificationsallow one to think of S as a set of functions, { (tY)y E S}, which form a basis of IZk (, S). Wewill sometimes refer to S as a refined grid and sometimes as a basis with this identificationunderstood.One should think of the Fk as spaces of coefficients for function expansions in the corre-sponding Ik spaces, in the basis S. The o simply associate a set of coefficients in Fk witha function in Ik. When 0 is understood, we may write just tk.128We are now in a position to state the basic theorems of interpolet expansions on uniformgrids.Theorem 5.2.6 Let O(x) be an interpolet on R n . Then each mapping Lk : Yk(S) -+ Zk(q, S)(k = 1, 2,...) is an isomorphism.Proof: Since the map Lk is surjective, it is only necessary to show that tk is injective.By the definition, tkV = 0 if and only if there exist v E .Fk such thatv(y)((x - y)/2k) +yESnDkv(y)O(x - y).v(y)((x - y)/2k - 1) + . + E:yESnD1Let z E SnCk. By (INT1), we have 0((z-y)/2k) = 6(z-y)/2k,o = 6z,y, for y e SnCk, andalso ((z - y)/2') = 0, for y E SnD,. So, (tkV)(Z) = V(Z),Z E Sn Ck, therefore v(z) = 0, z ESn Ck. This being so, one then has (Lkv)(z) = v(z), z ES n Dk, SO v(z) = 0, z c S n Dk.Once again, (tkV)(z) = v(z),z c Sn Dk-l, thus we must have v(z) = O,z E Sn Dk-1.Continuing in this manner, we deduce that v(y) = 0, y c S, thus v = 0.Corollary 5.2.7Ik(S) = k(S n Ck) Ik (Sn Dk) E '.' " ZEk(S n Di)Since the sum is direct, the expansion is unique.This corollary is a consequence of observations in the above proof.Theorem 5.2.8 Let O(x) be an interpolet on R . ThenVk, k (0, Z n ) = Ik-1(0, Z"n).Consequently,Vkl, k2,i kl (, Z n ) = Ik 2 (, Zn).1290= EyESnCkProof: To prove that 'k C 'k-1 we note that Ik C k-1 U k(Ck). Thus we just needto show Ik(Ck) C $k-1.Translating by z E Ck and inserting powers of 2 where appropriate, one can rewrite(INT2) for 0 as((x - z)/2k) = ((x - z)/2k-1) + Cy/2k-1 0(( - z - y)/2k-1).yEDkThe terms in the right hand side are elements Ik-1. Thus Ik C -k-1-To prove Ik-1 C Zk note that any element of Ik-1 can be expressed as:f(x)= EYECk-1a(y)q((x - y)/2 k - 1 ) + EyGDk-1a(y)o((x - y)/2k- 2 ) +---yCD 1All the terms in this expansion but the first are elements of Zk. Since Ck-1 = Ck U Dk wemay split the first sum up as,a(y)¢((x - y)/2 k - 1) = E a(y)cb((x- y)/2±k - 1 ) + : a(y)o((x - y)/2 k - 1)yEDkThe second term is also an element of k-.Rewriting (INT2) one has (y E Ck):((x - y)/2 k - 1 ) = ¢((x - y)/2k) - Cz/2k-14((x - z - y)/2k-1).zEDky E Ck, z E Dk, so y + z E Dk, thus the right hand side is made up of elements of Ik. Thus,Ik-1 C Zk130a(y)(x - y)EyECk-15.2.2 Interpolet transformsCorollary 5.2.9 Interpolet DecompositionThe set of isomorphisms Lk induces a set of isomorphismsJkl,k2 Yki -4 k 2Jkl,k2 = L Lki-We refer to these isomorphisms as interpolet transforms. It is our convention to let Jk = JO,kand J-k = Jk,o. The reader will note that the Ji are linear transformations on the coefficientspaces, and are thus the primary object of computation.We now turn to a study of the J's. It is clear from the definition that for kl < k2 ,Jkl,k2 = Jk 2-1,k 2 O Jk 2-2,k 2-1 O ... Jkl,kl+l, and similarly for kl > k2 . Thus, we need onlystudy the Jk,k+l and Jk+l,k mappings.Theorem 5.2.10 Computation TheoremLet v E Fk (Zn).(Jk,kV)(y) = v(y), y V Dk+1where v'(y) = v(y) - EECk+1 C(y-z)/ 2 k(Z).(Jk+l,kv)() = v (y), y B Dk+1v'(y), y Dk+1where v'(y) = v(y) + EzECk+1 c(y+z)/2kv(z).Proof: For v E Fk we haveLkV = E v(y)o((x - y)/2k) + v(y)((x - y)/2k-1) + . + v(y)( - y)yECk yEDk yED1131expanding the first term,tkV = E v(y)q((x - y)/2k) + E v(y)q((x - y)/2k) + E v(y)o((x - y)/2 k - 1 ) + .yECk+1 yEDk+1 yEDkusing (INT2), 0((x - z)/2k) = 0((x - z)/2 k +l ) - ZyEDk+1 C(y-z)/ 2 k ((X - y)/2k),tkV = v(y)q((x - y)/2k+1) + Z v'(y)q((x - y)/2k) + Z v(y)q((x - y)/2k-1) +...YECCk+1 yDk+1 yCDk- v(y), y C Dk+1v'(y), yE Dk+1where v'(y) = v(y) - ZZECk+1 C(y-z)/2kv(Z).The proof for Jk+l,k is similar.OSimilar to what one might get with wavelets, we see that we can compute the coefficientsof interpolet expansions on uniform lattices by a pyramid algorithm. Computationally,this procedure can be carried out by first computing the D 1 coefficients with Jo,1, then bycomputing the D2 coefficients from the C 1 data with J 1,2, and so on. In this sense, it is nodifferent from standard multiresolution decompositions.A feature of the interpolet decomposition is that the transformations all have a particularlower triangular form. That is, if we write v E Fk as a vector with its Ck+1 components first,its Dk+1 components second, and the rest of its components third, then the transformationtakes the form,I 0 0 0 0 1VCk+M I 0 0 0 VDk+lJk,k+lv = 0 0 I 0 0 VDk0 0 0 I 0The0 0 0 0 inverse, i obtained by replacing with -The inverse, Jk+l,k, is obtained by replacing M with -M.1325.3 Accuracy of interpolet approximationGiven a function f(x) on R", one can form an interpolet approximation to f by the formula:f (x) ~ f (y)(x - y) = tof,yECowhere f on the right hand side is thought of as a function restricted to Co = Z (a more cum-bersome but more precise notation is Loff(Y)})vezE). This approximation has the propertythat (tof)(z) = f(z),z E Zn.Starting from the expansion tof E T o (, Zn) one can construct equivalent expansions,tk(Jkf) E Zk(¢, Zn). The coefficients Jkf are referred to as the interpolet transform of thefunction f.If f(x) is sufficiently smooth, then we can expect that the coefficients (Jkf)(y), y E D1, 1 Ok(Z) + p), or Figure 5-3 (for p = 1).touch? z ES Z--SyES ok okyE Zn - S not ok okOne example of a 1-level touching good basis for one dimensional 3rd order interpolets137--r^r~-r-- ^ - - --I*X~1--- -L-----~' ---L~---l-L"-L-~I~----LI~ - 1_ _1111 -. _1- ~~~-~-~~~-~-~--pl__~__Ll____sP__ ~----I-- -1 ~Ia3=S=Z-So o o o o o O *0- o 0 o I0 1 2 3 4 5 6 7 8 10 12 14xFigure 5-4: Our example of the 1-level touching good basis in one dimension. Note that thetwo functions plotted do not touch. Tick marks denote the set S C Z.(the ones being used as an example) is the set of interpolets centered at the pointsS= cou U {-7-21-5- **, - 7 ..., ,.21(see Figure 5-4). The support of the interpolet on level 0 at y is [-3 + y, 3 + y], the union ofall the supports of the interpolets in S on level 0 is [-10, 10]. The support of an interpoleton level 1 at y is [-6 + y, 6 + y], thus the only interpolets on level 1 which touch theinterpolets on level 0 are precisely those ones at points -14,...,14, which are precisely theones included in S. No interpolet not included on level 1 touches an interpolet included onlevel 0 so the definition is satisfied. The argument proceeds similarly on higher levels. Inthree dimensions, this example corresponds to nested concentric cubes of size 15 ... 21 ateach level 1 < n. Figure 5-5 shows a randomly generated generic example.An example of a 2-level touching good basis for 3rd order is the set of interpolets centeredat the pointsnU -5 -2' -3- 27,...3 21, 2 '138o o 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 000000 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 000000 00 000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 000 0 0 0 000000o 0 0 0 0 0 0 0 0000000 0 0 0 000000000 0 0 0 0 000000000 0 0 0 00000000o 0 0 0 0 0 0 000000000 0 0 0 0 000000000 0 0 0 0000000 0 0 0 0 0 0000000 0 0 0 0 0 0 000 0 0 0 0 0 0 0 000000 0 0 0 0 0000000 0 0 0 0 0000000o 0 0 0 0 0 0 0 0000000 0 0 0 000000000 0 0 0 0 0 0000000 0 0 0 0000000S 0 0 0 0 0 0 0 0 0000000 0 0 00000 0 0 0 0 0 0 0 0 0 0 00 0 0 00000000000000 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0o O 0 0 0o O 0 0 00 0 0 0 00 0 0 0 0Figure 5-5: A generic example of a truncation which meets our definitions of good and 1-leveltouching. In black are the points of S c Z n .1390 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 0 0 0 0 0 0 0 0 0 00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0 0 0 00 0 0 0 0 0 000 0 0 0 0 0 000 0 0 0 0 0 000 0 0 0 0 0 000 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0Figure 5-6: An example of a 2-level touching good basis which can be used for a diatomicmolecule (two atomic cores). The pointssignificant.in S C Z n are where residual values could beNote (for n > 2) that this set is not 1-level touching since the level 1 interpolet centered aty = 12 is not included, while an interpolet it touches, namely the level 0 interpolet centeredat y = 5, is included. Figure 5-6 shows an example of a 2-level touching good basis.Finally, the set of pointsS= C U {-3- 2' -1. 2' 1- 2',3 21)l=0forms a good basis, but is not 1-level touching or 2-level touching, but it is 3-level touching.The above examples are meant to suggest that the good basis and the p-level touchingdefinitions can be thought of, informally, as conditions telling one how quickly one can changefrom one resolution to another. Essentially, any nested set of refined regions can satisfy theseconditions so long as the margins around the set of points at a given resolution are wide1400 000 000 00O0 000 0enough.From a computational point of view, what these conditions do is ensure that data-pathswhich carry coefficient information between different resolutions are not broken by zeroedcoefficients at intermediate levels.It is clear from the preceding discussion, in a good basis, one has Jkl,k2 = Jk 2 -1,k2 OJk2-2,k 2 -1 ... o Jkl,ki+l1 We may now generalize the computation theorem to a truncatedbasis.Theorem 5.4.4 Good Basis Computation TheoremLet S be a good basis. v E .Fk(S), y e Sand , is any member of the equivalence class of v(Jk,k+lf)(Y) f v(y), y E S - Dk+1Sv'(y), y e Dk+1where v'(y) = v(y) - EzEznCk+1 C(y-z)/2kv().(l,(Y) v(y), y E S - Dk+1Sv'(y), yE Dk+1where v'(y) = v(y) + EzEsnck+l C(y+z)/ 2 kv(Z).Proof:In a good basis, the computations of all Jkl,k2 are independent of the representative.Thus, this algorithm computed on v gives a member of the same class as would be computedon v.O]Thus, the pyramid algorithm of the uniform basis, Z, has a counterpart in a good basisS, allowing the computation of the expansion coefficients in .Fk(S) from the values of thefunction in To(S) (and also has the lower triangular structure).In a good basis, one thus has the ability to perfectly reconstruct the multiscale coefficientsof a function for the basis functions associated with the points of the refined grid S by simplyapplying the pyramid algorithm on zero-lifts at each stage of the algorithm. The above141I~__~, -.-...r------~-ir-~r~L-~-li -Il---~~lrx ~U' -~--L-...-~.-~~-. --- - I -theorem establishes this as true, even though we do not necessarily expect the data zeroedduring the lift to be small. (The function may have significant sample values throughoutthe domain of the representation). Also, with exact recovery of sample values, it is easy toperform local nonlinear point-wise operations of the form f(x) -4 G(f(x)) (e.g. ef(x)), orpoint-wise multiplication (i.e. f(x), g(x) - f(x)g(x)), despite the truncation.The reader may note that this result is "analysis-free" in that we have sparsified the com-putation, not by proving that the coefficients vanish outside of the truncation for some classof functions, but by showing the the coefficients we wish to compute have no dependency onthe coefficients we omitted. Computationally, this means the data-structure in the computerrequires no intermediate augmentations (contrast with [38]).The advantages conferred by the additional the p-level properties are seen in the contextof operator evaluation, and will be the subject of the next two sections.5.5 Multilevel algorithms for V2 and other operatorsGiven v, w E Fk (Zn), we may compute the stiffness matrix of a model operator, V2,KtkVJV2 tkW) = J(tkv)(X)V 2(tkw) (x)d Xby changing the expansions,(toJ-kV V2 LOJ-kW) Z (J-kv)(y)(Jkw)(Z) - Y)V 2 0(x - z)dnz.y,zE ZnThis reduces the computation to the computation of the matrix elements ( (x - (y - z)) V2 10(x))which can be done by solving the associated eigen-problem obtained by applying (INT2).In particular, let L - (,( - y) V2 1¢(x)), thenL= cz cZ2 ((2x- 2y- z1 -Z 2) + V2 (2x)zi,z 2 EZnL = 5 2 2 -nCziCz 2 Ly+z -Z2zi,z 2 EZn142which we solve by standard methods (found in [58], for example). In subsequent sectionsof this article, we will also define L + = (¢(x - y) V 2 1¢(x/2)), L ++ = (¢(x - y)lV21(x/4)),L- = (0((x - y)/2)|V 2 1(x)), and L = (0((x - y)/4) V 2 1q(x)), which can be computedfrom L by employing (INT2). Although it is true by Hermiticity that L- = L+ andL-- = L+ + , we will make no use of this fact.One can write a matrix expression (LkVIV 21LkW) = vttkLLJ-kW, where the L is theToeplitz matrix with coefficients L'. In practice, one typically formulates the preceding asa computation of u = JtkDJkw for some w E .Fk(Z). Then (tkV V 2 1kw)= vtu. u may bethought of as an element of Fk(Zn)*, the dual space of Yk(Zn). The task is then to computethe coefficients u(y), y C Z n = Ck E Dk E ... D 1. Any algorithm for computing xtAy can beadapted to an algorithm for Ay, and for purposes of making proofs, it is somewhat easier tokeep thinking of the computation as (Lkv V 2 LkW), which is the point of view we shall take.However, computing (tk' V72 LkW)= v t kLJkCV, by just applying the transforms, andthe Toeplitz matrix is problematic, since this process makes it necessary to either represent vand w on a uniform grid, or to compute a matrix element between each pair of functions in thetruncated expansion which touch. In the first case, one ends up with an O(N) computationfor (typically) a very large N. In the second case, one chooses between extremes which areO(N) and quite complicated or simple and O(N2 ).The following sections outline the design of multilevel algorithms for V 2 for both 1-leveland 2-level touching bases. Both algorithms are derived according to the following format:break up the expansion of (tkV V72 I kW) into a decomposition over elements at the samelevel and adjacent levels.rewrite the expansions in terms of the matrix elements between elements of those levelsand the transforms of higher/lower level elements.implement the algorithm by computing those terms separately.establish correctness in a p-level truncated basis.Although only the 1-level and 2-level algorithms have been explored in any detail, thissame process will generally work to produce O(N) p-level algorithms for any p.1435.5.1 V 2 in 1-level decompositionThe 1-level decomposition of (tkVIV 2~LkW) isLkV IV2 LkW = kIVv2 ItkW)± + KLkv IV2 LkW) + tkV V2 ItkW)where( kVj'V2kW)0tkVU2 LkW)+tLkVIVkW= SOk ()=k (z)= SOk()Ok (z)tLkjv2tLkz) v(Y)W(Z)KLkyv2 LkZ) y(Y)W(Z)KLkjV2 Lkz) V(Y)W(Z).That is, we express the product as contributions from levels to the same level, higher levels,and lower levels. We will now investigate each of these terms individually. In the languageof matrices, these respectively correspond to diagonal, above diagonal, and below diagonalblocks of the V2 matrix.L(kV2 1 kW)= V t() 0(I)-(I)-()+ ()+ (1)+()0 ()- (I)()- () (()-( (Diagonal blocks: (I)OFor some fixed 1, 0 < 1 < k a term of the formEOk(z)=Ok(y)=l(((x - y)/2 ) 2|¢((x - z)/2')) v(y)w(z)contributes to (1)OHowever, (((x - y)/2) IV2l((X - z)/2')) = 2(n-2) (q - y) (- z)) = 21(n-2)Loz-y"144Thus we have(tk v2k w)O = 2'(n--2) Lzv()y)w(z)Super-diagonal blocks: (J)'For some fixed 1, 0 < 1 < k a term of the formS (0((x - y)/21) I2((X - z)/2") v(y)w(z)0 k(z)=l'>I=Ok(y)contributes to ()+ .Applying the inverse transform (Jk,l+l) to the w coefficients in the sum allows us to writethis term as0((x - y)/2 ) V 2 $((x - z)/2'+1)) v(y)(Jk,l+lw)(z).1+1 (z)=1+1>=0k(Y)Thus we haveKkUlV2 kw)+ = 21(n-2) LE yv(y)(Jk,l+w)(z).Sub-diagonal blocks: (I)-For some fixed 1, 0 < 1 < k a term of the formEOk(Z)=1<1'=Ok(Y) (¢((x - y)/21') 21 ((x- z)/21)) v(y)w(z)contributes to ()-.This time applying the inverse transform (Jk,I+) to the v coefficients in the sum allowsus to write this term asOkE(z)=Ok(Z)Which corresponds to the matrix decomposition:Lk VIV2 LkW)The key idea is to evaluate the diagonal and first off diagonal blocks of the V 2 matrix andthen to compute the other blocks above and below the tridiagonal through the transforms.() o, (+1, (-, I++, and (I)--The definitions for the contributions in the decomposition proceed just as they did for the1-level case.(tkV V 2 LkW) ° = 21(n-2) Lo ,v(y)w(z)SOk()=Ok(y)= Z152wo (1)+ ()++ (1)++= V (r (I)° (I)+ (I)++(I)- (D- 0 (I +(I) (I)-- ( )- (DOK LkVIVjkkW )+K kVIVILkW)KLkVIV2ILkW)k tkVIV2 kW)-21(n-2)1= 21(n-2)1= 21(n-2)= 21(n-2)1S L-y,v(y)w(z)k (z)- = k (y)=lLv (y) (z)Ok(z)=Ok(y)-=1S Lzfv(y) (Jk,1+2W) (z)Ok(Z)-l>l= k (Y )S L-,(Jk,1+2v) (y)w(z)Ok(z)=1k (), +2(z)=+2 L_++Iwtmp(z)wtmp +- J+li,twtmpend forfor I = 0 to k - 2utmp +- Jt+l,iutmputmp(y) -- utmp(y) + 21(n- 2) Oz+2(y)=t+2,0k(2)=1 L -w(z)end foru -- u + utmp5.6 Efficient implementationWe have produced a very successful 3D implementation of all of the above algorithms for theinterpolet used as this chapter's example (m = 3). Implementation details are given in thissection. The ideas used to make this implementation efficient for all of the above algorithms.The purpose of this section is to give additional information to readers who wish toimplement these algorithms themselves.1545.6.1 Data structuresThe interpolet data and function samples are kept in a sequence of blocks at various levels.Each block at level k contains the points of a rectangular subset of Ck-1. Since Dk =Ck-1 - Ck, we use the collection of blocks at level k < p (p being the top level) to represent arectangular subset Ok, ignoring the Ck points of each of these blocks. In our implementation,these extra Ck points hold the value 0 in between operations and take on useful intermediatevalues during operations. Since we are working in 3 dimensions, this multiplies the storagerequired by a factor of about 7, which we found an acceptable cost for its advantages.The coefficients for the transforms and operators are kept in various 3D arrays. Al-though it is possible to build the coefficients upon demand from a set of 1D arrays ofcoefficients, we have found that the arithmetic cost of doing this is much greater than thecost of storing them (about 10 flops are required for each V2 coefficient, while the 3D ar-rays are still small enough to be stored in cache). We have (in Fortran notation) the filterscs(0:3,0:3,0:3), SAMELEVEL(0:5,0:5,0:5), ONELEVEL(0:8,0:8,0:8), and (for 2-level algo-rithms) TWOLEVEL(0:14,0:14,0:14) (note: we have made use of the fact that our operatorsare symmetric to cut the size of these arrays by a factor of 1 and use ONELEVEL andTWOLEVEL for both upward and downward inter-level communication).5.6.2 Implementation overviewThe blocks described in the previous section are used as the fundamental objects for manip-ulation. The computation proceeds by employing block-to-block subroutines for the variousoperations, having every block at a level send data to every block at the same level, one levelup or down, or (for 2-level algorithms) two levels up or down.The number of blocks at each level is not very large, and if a subroutine determinesthat the intersection of two blocks is empty (which it does by examining the boundingrectangles), then it returns immediately. Thus, while this algorithm is to be O(B 2) where Bis the number of blocks, it remains O(N) where N is the number of actual points, becauseB is much smaller than N.155-- --I -^L--^I-'^"I-I~"--".---LLIII 1111 - lsllllli5.6.3 Block-to-block subroutinesThe block-to-block subroutines are all designed to take two blocks (source and destination)and a set of filter coefficients and place the result of convolving the filter with the sourceblock in the overlapping points of the destination block. There is a block-to-block subroutinefor the interpolet transform, its transpose, its inverse, and its inverse transpose, as well asoperator application routines for the same-level operator, up-one-level operator, down-one-level operator, up-two-level operator, and the down-two-level operator.All of these routines precompute the bounding box for the points in the destinationblock which are in the range of influence of the source block and, for each point in thissub-block, the bounding box for the points in the source block in the domain of influenceof the destination point. The result of this precomputation is that the only data values ofthe source (destination) which are accessed are the ones which are read (modified). Thisdecreases the number of data accesses in our test problems by a factor of 7.Additionally, blocking the computation generally increases the locality of access for thedata. More data requests hit the cache than would occur in a more arbitrarily arrangedconstruction.5.7 Results and conclusionsWith interpolets, it is possible to carry out O(N) computations in truncated bases, whereN is the number of elements retained in the truncation, without having to augment thegrid of points associated with the functions maintained in the basis. Along with allowingone to compute common linear physical operations, interpolet algorithms also allow one totransfer between function values and multiscale expansion coefficients on grids of variableresolution recovering the same results as one would obtain working with data on a full gridof arbitrary resolution but without introducing additional grid points into the calculation.This allows local nonlinear couplings to be computed quickly without the introduction ofadditional sample points and without the introduction of additional approximations whichmust be controlled.156These algorithms have been implemented in Fortran90 and have subsequently beenadopted for use in electronic structure computations as described in the introduction. Priorto this, we had been using very simple, naive O(N 2 ) algorithms which implement each trans-form and operator as multiplication by the multiscale representation of the correspondingmatrix. These multiplies check all points for being within interaction range and then con-struct the appropriate matrix element as needed. This is required in the na'be approachbecause the variety of inter-scale matrix elements is too wide to store in table of reasonablesize. This algorithm ultimately scales quadratically with the number of refinement levels forour application. This is because, as described in the introduction, basis functions are keptin the basis whenever they contain an atomic nucleus within their support. All functionsin this subset of significant size of the basis functions associated with each atomic centertherefore touch one another, and the multiscale matrices contain dense blocks connectingall of these elements of a given center with one another. Because the number of functionsassociated with a given center grows linearly with the number of refinement scales k, thenumber of operations required in the naive approach of multiplying directly by these densematrices scales quadratically with the number of functions in the basis. For reference, a typ-ical number of refinement levels in electronic structure calculations of the lighter elementswould be k = 5, as employed in the carbon atom [5] and the nitrogen molecule [3].A comparison with the previous implementation in Fortran90 on the same processor(Superscalar SPARC Version 9, UltraSPARC) demonstrates the speed improvements andscaling which can be achieved with the new approach. The "time" axis is the CPU timetaken by one application of the V2 operator. The "k" axis represents the number of levelsof refinement made in the basis and is proportional to the number of points in S.Figure 5-7 compares the runtimes of V2 in three dimensions on a 1-level touching goodbasis with 3rd order interpolets consisting of concentric cubes of size 153 centered aboutone atomic nucleus, as would be appropriate for the calculation of the electronic structureof a single atom. Although there is initially a significant O(N) contribution, as a functionof the number of refinement levels k, the times for the naive approach show the constantincrements in slope characteristic of a quadratic function. The new approach compares veryfavorably and is about 40 times faster for typical values of k. (Note the difference in vertical157previously used laplacian 1-level laplaciano o4]E E3100501 2 3 4 5 6 1 2 3 4 5 6k kFigure 5-7: The previously used implementation is on the left, and the implementationemploying a 1-level touching algorithm is on the right. (Note the difference in scale on thevertical axes.)scale between the two figures.)Although the comparison in Figure 5-7 is quite favorable for the new algorithm, one mustbear in mind that given the typical decay in the interpolet expansion coefficients about anatom[5], [3], the functions which are appropriate to maintain in the expansions tend to havethe 2-level touching property, not the 1-level touching property. Figure 5-8 compares theruntimes of V 2 in three dimensions on a 2-level touching good basis of concentric cubes ofsize 93, where the speed up just as dramatic as before, now by approximately a factor of 30.Figure 5-9 compares the runtimes of V 2 in three dimensions on a 2-level touching goodbasis of two refinement centers, with refinements now consisting of cubes of size 93 (similarto figure 6). This situation arises in the the calculation of the electronic ground state of thenitrogen molecule, N2. Note that with the introduction now of two atomic centers the timesare again consistent with the scalings described above: The runs times only double in themultilevel algorithm but quadruple in the naive algorithm.158previously used laplacian 2-level laplacian3 4k3 4kFigure 5-8: The previouslyemploying a 2-level touchingvertical axes.)used implementation is on the left,algorithm is on the right. (Note theand the implementationdifference in scale on the159previously used laplacian 2-level laplacian150 43.53100"2.5E E/ I '"1.55010.50 1 C - -01 2 3 4 5 2 3 4 5 6k kFigure 5-9: The previously used implementation is on the left, and the implementationemploying a multilevel algorithm is on the right. (Note the difference in scale on the verticalaxes.)160Having considered the efficiency of the algorithms, we next turn to the use of thesealgorithms in the solution of Poisson's equation to determine electrostatic fields, which wasthe rate limiting step in the calculations carried out in [5] and [3]. From those calculations,we were aware that the combination of conjugate gradients with the simple preconditioningconsisting of just applying the inverse of the diagonal elements of the Laplacian matrix leadsto an algorithm requiring very few iterations. It was the application of the operator withineach conjugate gradient iteration which limited the efficiency of those earlier calculations.Figures 5-10,5-11,5-12 illustrate results for varying levels of refinement in the two differentsystems. The first system consists of refinements centered about a single atomic center withina cubic super cell of side 15 Bohr radii with periodic boundary conditions. (One Bohr radiusis approximately 0.529 Angstroms.) The second system contains two refinement centersseparated at a distance of 2 Bohr radii, approximately the inter-nuclear separation in thenitrogen molecule. This latter system resides within a rectangular supercell of dimensions(15 Bohr) 2 x (17 Bohr). In both cases, the spacing of the grid at the coarsest scale is 1 Bohr,and the finest spacing is 2 - k Bohr. At k = 22, the greatest refinement considered in ournumerical experiments, the finest grid spacing is approximately 0.24 x 10- 6A. A full gridat this resolution would contain 2.8 x 1023 points. Our truncated basis contains only about60,000 functions in this case.Figure 5-10 compares, as a function of the number of refinement levels k, the conditionnumber of the Laplacian represented in a truncated interpolet basis (the "stiffness matrix"for the basis) and in an untruncated orthogonal basis at the corresponding resolution. Thefigure also shows the effect on the condition number of the interpolet stiffness matrix of thesimple diagonal preconditioner described above. The condition numbers for the truncatedinterpolet bases were determined numerically using the operators implemented as describedabove. The curves indicate results for the system with a single atomic center, and thesymbols indicate results for the two atom system. Comparing the results for the one andtwo atom cases suggests that apart from some transient behavior for small k, the conditionnumber is not sensitive to the number of atoms and depends primarily on the number ofrefinement levels k.Although finite basis representations of the Laplacian constructed from orthogonal func-161~......~_r. ...l......~.~.--~P----. -C ^.III~I--_-~-~ClrI- P(PI^IIIII~ 1L-~- ~~~~-~~--11- I- C S~-IPII -~~- ~I~~~ _I 1 -AI\E: 4C1010100Condition of Stiffness Matrix-.- Orthogonal Laplacian- -,o Interpolet Laplacian-,x Precondtloned Interpolet Laplacian ,r AcAr.WoXox1 2 3 4 5Number of refinement levels ->6 7Figure 5-10: The condition number of the Laplacian operator represented in truncated mul-tiresolution interpolet basis as a function of the number of refinement levels k with andwithout simple diagonal preconditioning and compared with the condition number in an or-thogonal basis with the same resolution. Lines indicate results for bases with a single atomiccenter of refinement, and points represent results for two atomic centers corresponding tothe nitrogen molecule.162Convergence of Preconditioned Conjugate Gradients for N2A 10 -10a)1010-100 20 40 60 80 100 120Number of iterations ->Figure 5-11: Convergence of the solution to Poisson's equation for the nuclear potential in anitrogen molecule in an interpolet basis with k = 8 levels of refinement about each nucleus.Performance of Preconditioned Conjugate Gradients for N24.54-3 .5 .........-oE.021.5i0 5 10 15 20 25Number of refinement levels ->Figure 5-12: Effective condition number of Poisson's equation for the nuclear potential ina nitrogen molecule with simple diagonal preconditioning as a function of k, the number oflevels of refinement about each nucleus.163...... ~ ..~.x,..~ --------- -----~-- ------u --- -CY -)-----4a------'-, ~c~pY.~----~~x..l .--~ --C-"rPT -~~I"~II~--1~1.---., .CIIIWs~----~CIC- II LI~~-~~-- I-I --tions at a given resolution should all have similar condition numbers, the fact that the inter-polet basis is not orthogonal allows the condition numbers of multiscale interpolet operatorsto be quite different that their single-scale counter parts. Compared to an orthogonal basis,the condition number in the interpolet representation is already over two orders of magnitudesuperior at the typical k = 5 levels of refinement. This comparison continues to improvewith increasing scale. The orthogonal condition number scales inversely as the square of thespacing on the grid of maximum resolution whereas the interpolet condition number scalesinversely with approximately the 5/4 power of the resolution, as determined from the slope inthe figure. The interpolet basis itself therefore provides an intrinsic form of preconditioning.Figure 5-10 shows that our simple explicit, diagonal preconditioner improves the scaling ofthe condition number, which now scales merely as the inverse of the resolution. (Note thelower slope of the lower curve.) At k = 5 levels of refinement the improvement is only afactor of three but becomes more significant as the number of refinements increases.Figure 5-11 shows the convergence of the preconditioned conjugate gradient algorithmin the solution of Poisson's equation. As a simple example, we solve for the electrostaticpotential which arises from the two nuclei in a nitrogen molecule. In this calculation we usek = 8 levels of refinement, a somewhat higher level of resolution than would be employed incalculations of the N2 molecule. For this calculation, the charge density of each nucleus ismodeled as a three dimensional Gaussian of root mean square width a along each directionequal to the spacing on the finest scale. After an initial phase of about twenty iterations, theconvergence becomes nearly perfectly exponential. This procedure reduces the magnitude ofthe residual vector by ten orders of magnitude in one hundred iterations. This is very goodperformance for a system consisting of 14,000 degrees of freedom with a Laplacian operatorwith a nominal single-scale condition number of about 65,000 at this level of resolution.The slope of this exponential portion of the convergence curve corresponds to a reductionin error at each iteration by 25%. One would obtain the same error reduction in a simpleweighted iterative Jacobi algorithm (with the inverse of the maximum eigenvalue as theweight) applied to an operator with condition number c e 4. The quantity c, the inverseof the fractional improvement in the magnitude of the residual, we define as the effectivecondition number for the algorithm.164Figure 5-12 shows this effective condition number c for the conjugate gradient algorithmwith simple diagonal preconditioning as a function of the number of refinement levels k forthe solution of Poisson's equation for the nuclei in the nitrogen molecule. In all cases theextent of the nuclei a is again set to the spacing of the finest grid. We note that after aboutsix refinements, the effective condition number is essentially constant. The example fromFigure 5-11 is therefore representative of the typical rate of convergence attained. Theseresults indicate that, regardless of the amount of refinement, a constant number of iterationswill suffice to produce a result of a given accuracy, even as the nominal condition numberfor an orthogonal stiffness at the corresponding resolution approaches 1.8 x 1013 at k = 22.Because the computational work involved in each iteration is linear in the number of pointsin the basis, this approach appears to produce the solution to Poisson's equation in O(N)time for these multiresolution bases.165166Appendix AMatlab source for Stiefel-GrassmannoptimizationThis appendix presents the source code for the sgmin routines that were dicussed in chapter4. We have included the sources for the basic routines followed by the sources for the variousexamples.Sources from the top-level sg_min directory:function [fn,Yn] = sgmin(YO,varargin)% SGMIN Stiefel/ Grassmann minimization meant for template use.% [fn, Yn] = SG MIN(YO) minimizes F(Y) where F is% a matlab function defined in F.m (dF.m and ddF.m),% where Y satisfies Y'*Y=L% [fn, Yn] = SGMIN(YO, rc,mode,metric,verbose,gradtol,ftol,partition)% Required Argument: YO expected orthonormal, YO'*YO = I% Optional Arguments: (may be specified in nearly any order)167% rc={ 'real', 'complex'} specifies real vs complex computation.% mode={'frcg', 'newton'} picks search direction:% Fletcher-Reeves Nonlinear Conjugate Gradient% Newton's method (Linear CG Hessian inverter)% metric={ 'flat', 'euclidean', 'canonical'}% motion={ 'approximate', 'exact'}% verbose={ 'verbose', 'quiet'}% gradtol={ first of any scalar arguments } convergence tolerance% ftol ={ second of any scalar arguments } convergence tolerance 20% partition = cell array describing symmetries in F% Defaults:% rc: 'real' if isreal(YO), 'complex' otherwise% partition: automatically determined% SGMIN(YO,rc, 'newton', 'euclidean', 'approximate', 'verbose', e- 6, l e- 8, partition)% Output:% fn = function minimum% Yn = minimizing argument will satisfy Yn'*Yn=I.% role parses the arguments, sets the global parameters and calls the% minimizers 30if - exist('F', 'file'), error('F.m must be in matlab''s path'), endif - exist('dF', 'file'), error('dF.m must be in matlab' 's path'), endif - exist('ddF','file'),error('ddF.m must be in matlab''s path'), endglobal SGParameters;SGParameters = [];[YO,r] = qr(YO,O);SGParameters.verbose = 1; 40nas = length(varargin);metarg = 0; rcarg = 0; partarg = 0; ftolarg = 0; gradtolarg = 0;mdarg = 0; motarg = 0;for j=l:nas,if (ischar(varargin{j }))168if (strcmp(lower(varargin{j}), ' approximate'))SGParameters.motion = 0; motarg=1;elseif (strcmp(lower(varargin {j}), 'exact '))SGParameters.motion = 1; motarg=1; 50elseif (strcmp(lower(vararginf{j}), 'flat '))SGParameters.metric = 0; metarg=l;elseif (strcmp(lower(varargin{j }), 'euclidean'))SGParameters.metric = 1; metarg=1;elseif (strcmp(lower (varargin {j}), ' canonical'))SGParameters.metric = 2; metarg=l;elseif (strcmp(lower(varargin {j}), 'real'))SGParameters.complex = 0; rcarg=1;elseif (strcmp(lower(varargin{j }), 'complex'))SGParameters.complex = 1; rcarg=1; 60elseif (strcmp(lower(varargin{j }), 'quiet'))SGParameters.verbose = 0; verbarg=1;elseif (strcmp(lower (varargin {j }), 'verbose'))SGParameters.verbose = 1; verbarg=1;elseif (strcmp(lower(varargin{j}), 'frcg'))SGParameters.Mode = 1; mdarg=1;elseif (strcmp(lower(varargin{j }), 'newton'))SGParameters.Mode = 0; mdarg=1;endelseif (iscell(varargin{j})) 70part = varargin{j}; partarg=1;elseif (isnumeric(varargin{j }))if (gradtolarg)SGParameters.ftol = varargin{j}(1);ftolarg=1;elseSGParameters.gradtol = varargin{j} (1);gradtolarg=1;endend 80end169% SGParamters.complex = 0 for real and 1 for complex.if (rcarg)if (-isreal(YO)) & (1 -SGParameters.complex)warning('YO has imaginary part, but real computation has been declared.endelseSGParameters.complex = ~isreal(YO);end% SGParameters.metric = 0 for flat, 1 for euclidean, and 2 for canonical 90if (-metarg)SGParameters.metric = 1;end% SGParameters.motion = 0 for approximate, 1 for exactif (-motarg)SGParameters.motion = 0;endif (-gradtolarg)SGParameters.gradtol = le-6;end 100if (-ftolarg)SGParameters.ftol = le-8;endif (-mdarg)SGParameters.Mode = 0;end% Make a partition using a possible given one given oneif (partarg)SGParameters.partition = part;else 110SGParameters.partition = partition(YO);endSGParameters.dimension = dimension(YO);if (SGParameters.Mode == 1)170[fn,Yn] = sgfrcg(YO);else[fn,Yn] = sg_newton(YO);end 120function [fn,Yn]= sg_newton(Y)% SG_NEWTON(Y) Optimize the objective function, F(Y) over all% Y such that Y'*Y = L Employs a local iterative search with% initial point Y and terminates if the magnitude of the gradient% falls to gradtol*(initial gradient magnitude) or if the relative% decrease in F after some iteration is less than ftol.% [fn, Yn]= SG NEWTON(Y)% Y is expected to satisfy Y'*Y = I.% Yn will satisfy Yn'*Yn = L to% role high level algorithm, Newton's Methodsglobal SGParameters;gradtol = SGParameters.gradtol;ftol = SGParameters.ftol;if (SGParameters.verbose)global SGdata;flps = flops;SGdata=[];end 20g = grad(Y); mag = sqrt(ip(Y,g,g)); oldmag = mag;f = F(Y); oldf = 2*f; N = 0;if (SGParameters.verbose)SGdata = [];disp(sprintf('%s\ts\t\ts\t\ts\t\t%s', 'iter', 'grad', 'F(Y)',' flops' , 'step type'));SGdata(N+1,:) = [N mag f flops-flps];disp(sprintf(' %d\tYe\t%e\t%9d\ts ',N,mag,f,flops-flps, 'none'));endwhile (mag>eps) & (mag/oldmag>gradtol) & (abs(oldf/f-1)>ftol)N= N+1; 30171sdir = -g;gsdir = ip(Y,sdir,-g); sdir=sdir*sign(gsdir); gsdir=abs(gsdir);sa = fmin( 'Fline',-abs(f/gsdir),abs(f/gsdir),[O le-4],Y,sdir);fsa = F(move(Y,sdir,sa));Hsdir = dgrad(Y,sdir); sdirHsdir = ip(Y,sdir,Hsdir);sb = fmin ('Fline',-abs(gsdir/sdirHsdir) ,abs(gsdir/sdirHsdir),[O le-4],Y,sdir);fsb = F(move(Y,sdir,sb));ndir = invdgrad(Y,-g,gradtol*oldmag); 40gndir = ip(Y,ndir,-g); ndir = ndir*sign(gndir);na = fmin('Fline' ,-abs(f/gndir),abs(f/gndir),[O le-4],Y,ndir);fna = F(move(Y,ndir,na));Hndir = dgrad(Y,ndir); ndirHndir = ip(Y,ndir,Hndir);nb = fmin (' Fline' ,-abs(gndir/ndirHndir),abs(gndir/ndirHndir),[O le-4],Y,ndir);fnb = F(move(Y,ndir,nb));if (fsaf)if (SGParameters.verbose) disp(' sg_newton: fmin overshoot!'), endfprime = ip(Y,g,dir); 60while (newf>f)t = t^2*fprime/ (t*fprime+f-newf) /2;newf = F(move(Y,dir,t));endend172Y = move(Y,dir,t); oldf= f; f = newf;g=grad(Y); mag=sqrt(ip(Y,g,g));if (SGParameters.verbose)SGdata(N+1,:) = [N mag f flops-flps];disp(sprintf(' .d\te\te\t%9d\ts ' ,N,mag,f,flops-flps,steptype)); 70endendfn = f;Yn = Y;function [fn,Yn]= sgfrcg(Y)% SG_FRCG(Y) Optimize the objective function, F(Y) over all% Y such that Y'*Y = I. Fletcher-Reeves CG iterative search with% initial point Y and terminates if the magnitude of the gradient% falls to gradtol*(initial gradient magnitude) or if the relative% decrease in F after some iteration is less than ftol.% fn, Yn]= SG CG(Y)% Y is expected to satisfy Y'*Y = L% Yn will satisfy Yn'*Yn = I. 10% role high level algorithm, Fletcher-Reeves Methodglobal SGParameters;gradtol = SGParameters.gradtol;ftol = SGParameters.ftol;if (SGParameters.verbose)global SGdata;flips = flops;SGdata= ;end 20g = grad(Y); mag = sqrt(ip(Y,g,g)); oldmag = mag;f = F(Y); oldf = 2*f; N = 0;olddir = 0;if (SGParameters.verbose)SGdata = H;173disp(sprintf('7s\ts\t\ts\t\t%s', 'iter' 'grad', 'F(Y) ',' flops'));SGdata(N+1,:) = [N mag f flops-flps];disp(sprintf( '%d\te\tZe\t%9d' ,N,mag,f,flops-flps));endwhile (mag>eps) & (mag/oldmag>gradtol) & (abs(oldf/f-1)>ftol) 30N= N+1;dir = -g;if (N>1)Hessolddir = dgrad(Y,olddir);alpha = ip(Y,dir,Hessolddir)/ip(Y,olddir,Hessolddir);dir = dir-alpha*olddir;endgdir = ip(Y,dir,-g); dir = dir*sign(gdir);cga = fmin ( 'Fline',-abs(f/gdir),abs(f/gdir),[O le-4],Y,dir); 40fcga = F(move(Y,dir,cga));Hessdir = dgrad(Y,dir); dirHdir = ip(Y,dir,Hessdir);cgb = fmin('Fline',-abs(gdir/dirHdir),abs(gdir/dirHdir),[O le-4],Y,dir);fcgb = F(move(Y,dir,cgb));if (fcgaf)if (SGParameters.verbose) disp(' sg_frcg: fmin overshoot!'), endfprime = ip(Y,g,dir);while (newf>f)t = t^2*fprime/(t*fprime+f-newf)/2; 50newf = F(move(Y,dir,t));endend[Y,olddir] = move(Y,dir,t); oldf= f; f = newf;g=grad(Y); mag=sqrt(ip(Y,g,g));if (SGParameters.verbose)SGdata(N+1,:) = [N mag f flops-flps];disp(sprintf(' %d\tYe\te\t%9d' ,N,mag,f,flops-flps));endend 60174fn = f;Yn = Y;function G = grad(Y)% GRAD computes the gradient of the energy F at the point Y.% G = GRAD(Y)% Y is expected to satisfy Y'*Y=L.% G will satisfy G = tangent(Y,G)% role geometrized objective function, this is the routine called to produce%5 the geometric gradient of the unconstrained differential of F.% The analog is G = A*Y-b. 10G = tangent(Y,dF(Y));function W = dgrad(Y,H)% DGRAD Computes the tangent vector, W, which results from applying the% geometrically correct hessian at the stiefel point, Y, to the% tangent vector H.% W = DGRAD(Y,H)% Y is expected to satisfy Y'*Y = I% H is expected to satisfy H = tangent(Y,H)%W will satisfy W = tangent(Y, W)10% role geometrized objective function, the is the routine called to apply% the covariant hessian of F to a tangent vector. The analog is% W = A *H, here A is the hessian of F.global SGParameters;met = SGParameters.metric;if (met==O)W=tangent(Y,ddF(Y,H));elsedf = dF(Y); g = tangent(Y,df);175W = connection(Y,g,H)+... 20dtangent(Y,df,H)+...tangent(Y,ddF(Y,H));end% this line removes the grassmann degrees of freedom from YW = nosym(Y,W);function H = invdgradcg(Y,W,tol)% INVDGRAD_CG Inverts the operator dgrad. i.e. solves for H satisfyingdgrad(Y,H) = W. Uses a conjugate gradient algorithm witha tolerance of gep*norm(W), or a tolerance of tol if given.H = INVDGRAD CG(Y, W)% Y is expected to satisfy Y'*Y=I% W is expected to satisfy W = tangent(Y, W)% H will satisfy H = tangent(Y,H)% 10% role geometrized objective function, this is the function called when one% wishes to invert the geometric hessian. The analog is H = Hess\ W.global SGParameters;dim = SGParameters.dimension;gep = SGParameters.gradtol;b = W;% p = reshape(dgrad(Y, reshape(q,n,k)),n*k,1) is the black box which% applies the hessian to a column vector, q, to produce a column vector p. 20% Any iterative inversion algorithm that applies the matrix as a black% box can be used here.oldd = 0;x = 0*b;r = -b;rho = ip(Y,r,r);if (nargin==3)176gepr = tol;elsegepr = sqrt(ip(Y,r,r))*gep;endposdef=1;cn = 0;while (posdef & sqrt(ip(Y,r,r))>gepr & cn<2*dim)% application of the hessian.cn = cn+l;if (cn==l)d = -r;elsebeta = rho/oldrho;d = -r+beta*oldd;% It is important% off the tangentendto make sure long term roundoff does not blow ussurface. So we reproject.d = tangent(Y,d);Ad = dgrad(Y,d);dAd = ip(Y,d,Ad);if (dAd<0)posdef=0;elsedist = rho/dAd;x = x+dist*d;r = r+dist*Ad;oldd = d;oldrho = rho;rho = ip(Y,r,r);endendif (SGParameters.verbose & posdef==0)disp(' invdgrad: Hessian not positive definite, CG terminating early');if (cn==l) x = -b; endend177if (SGParameters.verbose & cn==2*dim)disp(' invdgrad: max iterations reached inverting the hessian by CG'),if (cn==l) x = -b; endendH = tangent(Y,x);function H = tangent(Y,D)% TANGENT produces tangent H from unconstrained differential D.% H will be the unique tangent vector such that for any% tangent W, real(trace(D'*W)) = ip(Y,W,H).H = TANGENT(Y,D)% Y is expected to satisfy Y'*Y=I% D is expected to be the same size as Y% H will satisfy H = tangent(Y,H)% 10% role geometric implementation, this function helps produce the gradient and%covariant hessian.global SGParameters;met = SGParameters.metric;vert = Y'*D;verts = (vert+vert')/2;if (met==O)H = D - Y*verts;elseif (met==l)H = D - Y*verts; 20elseif (met==2)H = (D - Y*vert')/2;endfunction T = dtangent(Y,H,dY)% DTANGENT computes the differential of the tangent map. That is% T = d/dt tangent(params, Y+t*dY,H).178% T = DTANGENT(Y,H,dY)% Y is expected to satisfy Y'*Y = I% H is expected to satisfy H = tangent(Y,H)% dY is expected to be the same size as Y% T will be the same size as Y10o% role geometric implementation, helps to produce a geometricallycorrect covariant hessian.global SGParameters;met = SGParameters.metric;vert = Y'*H; verts = (vert+vert')/2;dvert = dY'*H; dverts = (dvert+dvert')/2;if (met==O)T = -dY*verts-Y*dverts;elseif (met==l)T = -dY*verts-Y*dverts; 20elseif (met==2)T = (-dY*vert'-Y*dvert')/2;endfunction [Yo,Ho] = move(Yi,Hi,t)% MOVE moves the point Yi and direction Hi along a geodesic of length t% in the metric to a point Yo and a direction Ho.% [Yo,Ho] = MOVE(Yi,Hi,t)% Yi is expected to satisfy Yi'*Yi = I% Hi is expected to satisfy Hi = tangent(Yi,Hi)% Yo will satisfy Yo'*Yo = I% Ho will satisfy Ho = tangent(Yo,Ho)% o10% role geometric implementation, the analog to Yo = Yi+Hi*t, Ho = Hi.global SGParameters;met = SGParameters.metric;mot = SGParameters.motion;[n,k] = size(Yi);179--I- --~^- - --~~"--PI"4;minXlimII IIIIII MNI INN 11 0, 11iiif (t==O)Yo = Yi;if (nargout==2)Ho = Hi;end 20return;endif (met==O)% Move by straight lines with a qr projection back to the manifoldif (nargout==2) magl = sqrt(ip(Yi,Hi,Hi)); end[Yo,r] = qr(Yi+t*Hi,O);if (nargout==2)Ho = tangent(Yo,Hi);mag2 = sqrt(ip(Yo,Ho,Ho));Ho = Ho*magl/mag2; 30endelseif (met==1)if (mot==0)% This section computes approximate euclidean geodesics% using polar decompositionif (nargout==2) magl = sqrt(ip(Yi,Hi,Hi)); end[U,S,V] = svd(Yi+t*Hi,O); Yo = U*V';if (nargout==2)Ho = tangent(Yo,Hi);mag2 = sqrt(ip(Yo,Ho,Ho)); 40Ho = Ho*magl/mag2;endelse% This section computes exact euclidean geodesicsif (nargout==2) magi = sqrt(ip(Yi,Hi,Hi)); enda = Yi'*Hi;[q,r] = qr(Hi-Yi*a,O);mn = expm(t*[2*a, -r'; r, zeros(k)]);nm = expm(-t*a);Yo = (Yi*mn(l:k,l:k) + q*mn(k+l1:2*k,l:k))*nm; so180[Yo,r] = qr(Yo,O);if (nargout==2)q =(Yi*mn(1:k,k+1:2*k)+q*mn(k+1:2*k,k+1:2*k));Ho = Yi*a+q*r*nm;Ho = tangent(Yo,Hi);mag2 = sqrt(ip(Yo,Ho,Ho));Ho = Ho*magl/mag2;endendelseif (met==2) 60if (mot==0)% this section computes approximate canonical geodesics using approximate% matrix inverse.if (nargout==2) magl = sqrt(ip(Yi,Hi,Hi)); enda = Yi'*Hi;[q,r] = qr(Hi-Yi*a,O);geo = t*[a, -r'; r, zeros(k)];mn = (eye(2*k)+geo/2)/(eye(2*k)-geo/2);mn = mn(:,l:k);Yo = Yi*mn(l:k,:) + q*mn(k+1:2*k,:); 70if (nargout==2)Ho = Hi*mn(l:k,:) - Yi*(r*'mn(k+1:2*k,:));Ho = tangent(Yo,Hi);mag2 = sqrt(ip(Yo,Ho,Ho));Ho = Ho*magl/mag2;endelse% This section computes exact canonical geodesicsif (nargout==2) magl = sqrt(ip(Yi,Hi,Hi)); enda - Yi'*Hi; so[q,r] = qr(Hi-Yi*a,O);geo = t*[a, -r'; r, zeros(k)];mn = expm(geo);mn = mn(:,l:k);Yo = Yi*mn(l:k,:) + q*mn(k+1:2*k,:);181YIIIY litif (nargout==2)Ho = Hi*mn(l:k,:) - Yi*(r'*mn(k+1:2*k,:));Ho = tangent(Yo,Hi);mag2 = sqrt(ip(Yo,Ho,Ho));Ho = Ho*magl/mag2; 90endendendfunction f = Fline(t,Y,H)% FLINE the energy, e, along the geodesic passing through Y% in the direction H a distance t.% f= FLINE(t, Y,H)% Y is expected to satisfy Y'*Y = I.% H is expected to satisfy H = tangent(Y,H)% role high level algorithm, basic window dressing for the fmin functionf = F(move(Y,H,t)); 10function C = connection(Y,H1,H2)% CONNECTION Produces the christoffel symbol, C, for either the canonicalor euclidean connections at the point Y.% C = CONNECTION(Y,H1,H2)% Y is expected to satisfy Y'*Y = I% Hi and H2 satisfy H1=tangent(Y,H1) and H2=... (similarly)% C will be a matrix of the same size as Y% role geometric implementation, an important term in computing the hessian toand in defining the geodesics.global SGParameters;met = SGParameters.metric;if (met==O)182% the unconstrained connectionC = zeros(size(Y));elseif (met==l)% the euclidean connection for stiefelC = Y*(H1'*H2+H2'*H1)/2;elseif (met==2) 20% the canonical connection for the stiefelb = H1'*H2-H1'*Y**Y'*H2;C = (H1*H2'*Y+H2*H1'*Y)/2+Y*(b+b')/2;endfunction i = ip(Y,H1,H2)% IP computes the inner produce of H1,H2 which are tangents at the% stiefel point Y.% i = IP(Y,H1,H2)% Y is expected to satisfy Y'*Y=I% H1,y are expected to satisfy H1 = tangent(Y,H1), H2 = tangent(Y,H2)% role geometric implementation, the analog of real(H1'*H2)global SGParameters; 1omet = SGParameters.metric;if (met==O)% unconstrained metrici = sum(sum(real(conj(H1).*H2)));elseif (met==1)% euclidean metrici = sum(sum(real(conj (H1).*H2)));elseif (met==2)% canonical metrici = sum(sum(real(conj(H1).*H2)))-sum(sum(real(conj(Y'*H1).*(Y'*H2))))/2; 20endfunction H = nosym(Y,H)183-- -- -- '' -~l - " Ift'' s il I jiIIIIIY IIIIYII~r% NOSYM Orthogonal projection to remove the block diagonal components% of the tangent vector H corresponding to the block diagonal% right symmetries of F(Y).% H = NOSYM(Y,H)% Y is expected to satisfy Y'*Y=I% H is expected to and sill satisfy H = tangent(Y,H)% role geometrized objective function, necessary to reduce the number of 10% dimensions of the problem and to have a well conditioned hessian.% Somewhat analogous to projecting H to a feasible set of search% directions.global SGParameters;part = SGParameters.partition;vert = Y'*H;for j = 1:length(part),H(:,part {j})=H(:,part{j}) -Y(:,part{j})*vert(part{j} ,part{j});endfunction part = partition(YO)% PARTITION makes a guess at the partition of the function F% based on the properties of dF at a randomly selected point Y.% part = partition(YO);% YO is expected to satisfy YO'*YO = I% P = size(YO,2)% role geometrized objective function, determines the feasible subspace% of possible search directions for the function F. 10global SGParameters;cm = SGParameters.complex;[N,P] = size(YO);[Y,r] = qr(randn(N,P)+i*cm*randn(N,P),O);A = Y'*dF(Y);As = (A+A')/2;184Aa = (A-A')/2;M = (abs(Aa)<1e-7*abs(As));scorebd = ones(1,P); np=O; 20for j=l:P,if (scorebd(j))np=np+l;part{np} = find(M(j,:));scorebd(part np}) = O*part{np};endendMcomp = O*M;for j=l1:length(part), 30Mcomp(part {j },part {j }) = ones(length(part {j}));endgoodpart = (sum(sum(abs(Mcomp-M)))==0O);if ~goodpartwarning('unable to find consistent partition of F.');part = num2cell(l:P);endfunction dim = dimension(YO)% DIMENSION correctly counts the dimension, dim (stiefel dimension minus% the grassmann degrees of freedom).% dim = DIMENSION% role geometrized objective function, used to set the parameters describing%o the particular submanifold of the stiefel manifold on which the% computation will occur.global SGParameters; 10part = SGParameters.partition;cm = SGParameters.complex;185[N,P] = size(YO);np = length(part);if (-cm)dim = N*P;dim = dim - P*(P+1)/2;for i=l:np,k = length(part{i}); 20dim = dim - k*(k-1)/2;endelsedim = 2*N*P;dim = dim - P ^2 ;for i=l:np,k = length(part{i});dim = dim - k^2;endend 30186Sources from the sg_min/finitediff subdirectory:function df = dF(Y)% dF Computes the differential of F, that is,% df satisfies real(trace(H'*df)) = d/dx (F(Y+x*H)).% df = DF(Y)% Y is expected to satisfy Y'*Y = I% df is the same size as Y% role objective function, this is the routine called to compute the% differential of F. 10[N,P] = size(Y);ep = le-6;for k=1l:P, for j=1:N,Yp = Y; Yp(j,k) = Yp(j,k)+ep;Ym = Y; Ym(j,k) = Ym(j,k)-ep;df(j,k) = (F(Yp)-F(Ym))/ep/2;end, endfunction ddf = ddF(Y,H)% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Y+x*H).% ddf = DDF(Y,H)% Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y% role objective function, the function is called to apply the 10/0 unconstrained hessian of F to a vector.ep = le-6;ddf = (dF(Y+ep*H)-dF(Y-ep*H))/ep/2;187Sources from the sg_min/examples/procrustes subdirectory:function f = F(Y)% F Computes the energy associated with the stiefel point Y.% In this case, by the equation f = l|AY- YB|I 12/2.f = F(Y)% Y is expected to satisfy Y'*Y=I% role objective function, the routine called to compute the objective%function.global FParameters; 10A = FParameters.A;B = FParameters.B;q = A*Y-Y*B;f = sum(sum(real(conj(q).*(q))))/2;function df = dF(Y)% dF Computes the differential of F, that is,% df satisfies real(trace(H'*df)) = d/dx (F(Y+x*H)).df = DF(Y)% Y is expected to satisfy Y'*Y = I% df is the same size as Y% role objective function, this is the routine called to compute the% differential of F. 10global FParameters;A = FParameters.A;B = FParameters.B;q = A*Y-Y*B;df = A'*q-q*B';function ddf = ddF(Y,H)188% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Y+x*H).% ddf = DDF(Y,H)% Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y% role objective function, the function is called to apply the to% unconstrained hessian of f to a vector.global FParameters;A = FParameters.A;B = FParameters.B;dq = A*H-H*B;ddf = A'*dq-dq*B';function Y = guess()% GUESS provides the initial starting guess Y for energy minimization.% In this case it makes a somewhat random guess.% Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF.global FParameters; 10inA = FParameters.A;inB = FParameters.B;[n,n] = size(inA);[p,p] = size(inB);Y = randn(n,p);Y = (inA*Y)/inB;[Y,r] = qr(Y,O);189----- W MM-N-AMN A111 h, I U Ifunction parameters(A,B)% PARAMETERS initializes the parameters for an instance of aminimization problem.% PARAMETERS(A,B)% A is expected to be a square matrix% B is expected to be a square matrix% role sets up the global parameters used at all levels of computation.10global FParameters;FParameters = [];FParameters.A = A;FParameters.B = B;function [A,B] = randprob(A = randn(10);B = randn(4);190Sources from the sg_min/examples/jordan subdirectory:function f = F(Y)% F Computes the objective value associated with the stiefel point Y.% In this case, by the equation f = IAY- YB(Y)II ^2/2.% f = F(Y)Y is expected to satisfy Y'*Y=I% role objective function, the routine called to compute the objective% function.global FParameters; 10A = FParameters.A;Ab = Block(Y);f = A*Y - Y*Ab; f = sum(sum(real(conj(f).*f)))/2;function df = dF(Y)% dF Computes the differential of F, that is,df satisfies real(trace(H'*df)) = d/dx (F(Y+x*H)).% df = DF(Y)Y is expected to satisfy Y'*Y = I% df is the same size as Y% role objective function, this is the routine called to compute the% differential of F. 10global FParameters;A = FParameters.A;Ab = Block(Y);f = A*Y - Y*Ab;df = A'*f - f*Ab';function ddf = ddF(Y,H)% ddF Computes the second derivative of F, that is,191-- --ddf = d/dx dF(Y+x*H).ddf = DDF(Y,H)Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y%o role objective function, the function is called to apply the 10% unconstrained hessian of F to a vector.global FParameters;A = FParameters.A;Ab = Block(Y);dAb = dBlock(Y,H);f = A*Y - Y*Ab;df = A*H - H*Ab - Y*dAb;ddf = A'*df - df*Ab' - f*dAb';function Y = guess()% GUESS provides the initial starting guess Y for energy minimization.Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF.global FParameters;inA = FParameters.A; 10blocks = FParameters.blocks;eigs = FParameters.eigs;[n,k] = size(inA);[ne,ns] = size(blocks);off = 0;Q = eye(n); A = inA;192for i=l:ne, for j=l:ns,if (blocks(ij) -= 0)[u,s,v] = svd(A(off+1:n,off+ l1:n)-eye(n-off) *eigs(i)); 20v = v(:,n-off:-1:1);q = [eye(off) zeros(off,n-off); zeros(n-off,off) v];A = q'*A*q;Q = Q*q;s = diag(s); s = flipud(s);off = off + blocks(i,j);endend, endY = Q(:,l:off);function [B,eigvs] = Block(Y)% BLOCK computes the matrix B appearing in the energy function% F = IIAY- YB(Y)II ^2 for jordan block problems.% [B,eigvs] = Block(Y)% Y is expected to be a stiefel point (Y'*Y= I)% B will be an upper triangular matrix in staircase form the%size of Y'*Y% eigvs will be the eigenvalues along the diagonal of B. If% FParameters.type == 'orbit' then eigvs == FParameters.eigs. o10% role objective function, auxiliary routine used by F, dF, ddFglobal FParameters;A = FParameters.A;eigvs = FParameters.eigs;blocks = FParameters.blocks;bundle = FParameters.bundle;[n,k] = size(Y);B = Y'*A*Y; 20[ne,ns] = size(blocks);off = 0;193for i=l:ne,if (bundle)eigvs(i)= trace(B (off+1 :off+sum(blocks(i,:)),off+1 :off+sum(blocks(i,:))))/sum(blocks(i,:),endfor j=l:ns,if (blocks(i,j)>O)B(off+ 1:k,off+1 :off+blocks(i,j))= eigvs(i)*eye(k-off,blocks(i,j));off = off + blocks(i,j); 30endendendfunction B = dBlock(Y,H)% dBLOCK This is an auxiliary procedure used by ddF which computes% B = d/dx Block(Y+x*H)% B = DBLOCK(Y,H)% Y is expected to be a stiefel point (Y'*Y= I)% H is expected to satisfy H = tangent(Y,H)% B will be a staircase matrix the size of Y'*Y% role objective function, auxiliary function for ddF 10global FParameters;A = FParameters.A;eigvs = FParameters.eigs;blocks = FParameters.blocks;bundle = FParameters.bundle;[n,k] = size(Y);B =H'*A*Y+Y'*A*H;[ne,ns] = size(blocks);off = 0; 20for i=l:ne,if (bundle)deigvs(i)= trace(B (off+1 :off+sum(blocks(i,:)),off+1 :off+sum(blocks(i,:))))/sum(blocks(i,:)194elsedeigvs = O*eigvs;endfor j=l:ns,if (blocks(ij)>O)B(off+1 :k,off+1:off+blocks(i,j))= deigvs(i)*eye(k-off,blocks(i,j));off = off + blocks(i,j); 30endendendfunction parameters(A,eigvs,serge,type)% PARAMETERS initializes the parameters for an instance of aminimization problem. This functions must be executedbefore any attempts to do any jordan structure minimzations.% PARA ME TERS(A, eigvs, serge, type)% A is expected to be a square matrix% eigvs is expected to be a lxp matrix of eigenvalues, assumed to% be distinct.% serge is expected to be a pxq matrix of serge characteristics. o10% structures (zero padded).% for example, if one wishes to describe a matrix that% has a substructure ofSJ_4(1 .0)xJ 2(1.0)xJ 2(1.0)xJ3(2.0)xJ 3(2.0)%o then one has% eigvs = [1.0 2.0]serge = [4 2 2; 3 3 0]f% type is expected to be the string 'bundle' or 'orbit'% role sets up the global parameters used at all levels of computation. 20global FParameters;FParameters = 0;FParameters.A = A;195010111JUNIX1P IIIIYYII W INII iFParameters.eigs = eigvs;serge = fix(max(serge,0*serge));% convert the serge characteristics into weyr characteristics.blocks = [];for j = 1:size(serge,1),for k=l1:max(serge(j,:)), 30blocks(j,k) = sum(serge(j,:)>=k);endendFParameters.blocks = blocks;if (strcmp(lower(type), 'orbit'))FParameters.bundle = 0;elseif (strcmp(lower(type), 'bundle'))FParameters.bundle = 1;else'Either bundle or orbit must be specified.', 40endfunction [M,eigvs] = post(Y)% POST Computes the matrix which corresponds to the one nearest to% FParameters.A with the jordan structure FParameters.blocks% and eigenspace Y.% [M,eigvs] = post(Y)% Y is expected to be a stiefel point (Y'*Y=I)% M will be the same size as FParameters.A with jordan structure% encoded by FParameters.blocks, eigenvalues of eigvs, and eigenspace% of Y. 10% eigvs will be the vector of jordan eigenvalues of M.global FParameters;A = FParameters.A;[B,eigvs] = Block(Y);M = A - (A*Y-Y*B)*Y';196Sources from the sg_min/examples/tracemin subdirectory:function f = F(Y)% F Computes the energy associated with the stieffel point Y.% In this case, by the equation f = trace(Y'AY)/2.% f = F(Y)% Y is expected to satisfy Y'*Y=I% role objective function, the routine called to compute the objective%o function.global FParameters; 10A = FParameters.A;f = sum(sum(real(conj(Y).*(A*Y))))/2;function df = dF(Y)% dF Computes the differential of F, that is,de satisfies real(trace(H'*df)) = d/dx (F(Y+x*H)).% df = DF(Y)% Y is expected to satisfy Y'*Y = I% df is the same size as Y% role objective function, this is the routine called to compute thedifferential of F. 10global FParameters;A = FParameters.A;df = A*Y;function ddf = ddF(Y,H)% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Y+x*H).ddf = DDF(Y,H)197% Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y% role objective function, the function is called to apply the 10% unconstrained hessian of F to a vector.global FParameters;A = FParameters.A;ddf = A*H;function Y = guess(P)% GUESS provides the initial starting guess Y for energy minimization.% Y will be nxP where FParameters.A is an nxn matrix. Uses a% random matrix.% Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF. 10global FParameters;inA = FParameters.A;[n,n] = size(inA);[Y,r] = qr(randn(n,P),0);function parameters(A)% PARAMETERS initializes the parameters for an instance of a% minimization problem.% PARAMETERS(A)% A is expected to be a square matrix% role sets up the global parameters used at all levels of computation.198global FParameters; 1oFParameters = U;FParameters.A = A;function K = Kinetic(n)K = 2*eye(n) - diag(ones(1,n-1),-1) - diag(ones(1,n-1),l);199magn owlMgiollll lumu llati I1Sources from the sg_min/examples/ldatoy subdirectory:function f= F(Y)% F(Y) computes the energy of a toy Ida particle in a box configuration.% F(Y) = 1/2*trace(Y'*A *Y)+c*1/4 *sum(sum(Y. ^ 2,2). ^2,1).% f = F(Y)% Y is expected to satisfy Y'*Y=I% role objective function, the routine called to compute the objective% function.global FParameters; 10A = FParameters.A;c = FParameters.c;f = trace(Y' *(A*Y))/2+c*sum(sum(Y. ^2,2).^2,1)/4;function df = dF(Y)% dF Computes the differential of F, that is,% df satisfies real(trace(H'*df)) = d/dx (F(Y+x*H)).% df = DF(Y)% Y is expected to satisfy Y'*Y = I% df is the same size as Y% role objective function, this is the routine called to compute the% differential of F. 10global FParameters;A = FParameters.A;c = FParameters.c;df = (A*Y)+c*( (sum(Y.^ 2,2)*ones(1,size(Y,2))).*Y);function ddf = ddF(Y,H)% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Y+x*H).200% ddf = DDF(Y,H)% Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y% role objective function, the function is called to apply the to% unconstrained hessian of F to a vector.global FParameters;A = FParameters.A;c = FParameters.c;ddf= (A*H);ddf = ddf+c*( (sum(Y.^ 2,2)*ones(1,size(Y,2))).*H );ddf = ddf+c*( (sum(2*Y.*H,2)*ones(1,size(Y,2))).*Y);function Y = guess(p)% GUESS provides the initial starting guess Y for energy minimization.% Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF.global FParameters;A=FParameters.A; 10n = length(A);[V,D]=eig(A); D=diag(D);[D,I] = sort(D); V = V(:,I);Y = V(:,l:p);function parameters(A,c)% PARAMETERS initializes the parameters for an instance of a% toy lda problem with laplacian A and self coupling%constant c.201% PARAMETERS(A,c)% A is expected to be a square symmetric matrix% role sets up the global parameters used at all levels of computation.10global FParameters;FParameters = [;FParameters.A = A;FParameters.c = c;function K = Kinetic(n)K = 2*eye(n) - diag(ones(1,n-1),-1) - diag(ones(1,n-1),1);202Sources from the sg_min/examples/simschur subdirectory:function f = F(Q)% F Computes the objective value associated with the stiefel point Q.% In this case, by the equationS f = sum(sum(qAql. ^ 2))/2+sumsmsum(qBql. ^2))/2;%0 where qAql = masked out part of Q'*A *Q and qBql similarly.% f = F(Q)% Q is expected to satisfy Q'*Q=I% role objective function, the routine called to compute the objective 10% function.global FParameters;A = FParameters.A;B = FParameters.B;mask = FParameters.Mask;qAq = Q'*A*Q;qBq = Q'*B*Q;qAql = (1-mask).*qAq;qBql= (1-mask).*qBq;f = sum(sum(qAql.^2))/2+sum(sum(qBql.^2))/2; 20function df = dF(Q)% dF Computes the differential of F, that is,% df satisfies real(trace(H'*df)) = d/dx (F(Q+x*H)).% df = DF(Q)% Q is expected to satisfy Q'*Q = I% df is the same size as Q% role objective function, this is the routine called to compute the% differential of F. 10global FParameters;A = FParameters.A;2031,, 1 ' l m 1lulluB = FParameters.B;mask = FParameters.Mask;qAq = Q'*A*Q;qBq = Q'*B*Q;qAql = (1-mask).*qAq;qBql = (1-mask).*qBq;df = Q*(qAq' *qAql+qAq*qAql' + qBq' *qBql+qBq*qBql');function ddf = ddF(Q,dQ)% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Q+x*dQ).% ddf = DDF(Q,dQ)% Q is expected to satisfy Q'*Q = I% dQ is expected to be the same size as Q% ddf will be the same size as Q% role objective function, the function is called to apply the to% unconstrained hessian of F to a vector.global FParameters;A = FParameters.A;B = FParameters.B;mask = FParameters.Mask;qAq = Q'*A*Q;qBq = Q'*B*Q;qAql = (1-mask).*qAq;qBql= (1-mask).*qBq;dqAq = dQ'*A*Q+Q'*A*dQ; 20dqBq = dQ'*B*Q+Q'*B*dQ;dqAql = (1-mask).*dqAq;dqBql = (1-mask).*dqBq;ddf = dqAq' *qAql+dqAq*qAql'+ dqBq' *qBql+dqBq*qBql';ddf = ddf+qAq' *dqAql+qAq*dqAql' + qBq' *dqBql+qBq*dqBql';ddf = Q*ddf;ddf = ddf+dQ*(qAq' *qAql+qAq*qAql' + qBq' *qBql+qBq*qBql');204function Y = guess()% GUESS provides the initial starting guess Y for energy minimization.% Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF.global FParameters;A = FParameters.A;B = FParameters.B;% use average[Y,T] = schur(A+B);% using random guessY = randn(size(Y));% use identity% Y = eye(size(A));function parameters(A,B)% PARAMETERS initializes th% a Schilders problem.% PARAMETERS(A,B)% A,B are expected to be real%o% rolee parameters for an instance of amatrices of the same size.sets up the global parameters used at all levels of computation.global FParameters;FParameters.A = A;FParameters.B = B;n = length(A);Mask = triu(ones(n));for 1=1:2:n-1,Mask(l+l,1)=l;205endFParameters.Mask = Mask;206Sources from the sg_min/examples/indscal subdirectory:function f = F(Y)% F Computes the energy associated with the stiefel point Y.% In this case, by the equation f = sumi IISi- Y D i Y1 ^2/2.% f = F(Y)% Y is expected to satisfy Y'*Y=I% role objective function, the routine called to compute the objective% function.global FParameters;Ss = FParameters.Ss;df = 0;for k=1:size(Ss,3),S = Y'*Ss(:,:,k)*Y; D = max(diag(S),0);f = sum(sum( (S - diag(D)).^2 ))/2;endfunction df = dF(Y)% dF Computes the differential of% df satisfies real(trace(H'*df))% df = DF(Y)Y is expected to satisfy Y'*Y% df is the same size as Y% rolthat is,d/dx (F(Y+x*H)).e objective function, this is the routine called to compdifferential of F.global FParameters;Ss = FParameters.Ss;df = 0;for k=1:size(Ss,3),S = Y'*Ss(:,:,k)*Y; D = max(diag(S),0);df = df+2*Y*S*(S - diag(D));207ute the= Iendfunction ddf = ddF(Y,H)% ddF Computes the second derivative of F, that is,% ddf = d/dx dF(Y+x*H).% ddf = DDF(Y,H)Y is expected to satisfy Y'*Y = I% H is expected to be the same size as Y% ddf will be the same size as Y% role objective function, the function is called to apply the% unconstrained hessian of F to a vector.global FParameters;Ss = FParameters.Ss;ddf = 0;for k=1:size(Ss,3),S = Y'*Ss(:,:,k)*Y; D =max(diag(S),0);dS = H'*Ss(:,:,k)*Y+Y'*Ss(:,:,k)*H;dD = (D>0).*diag(dS);ddf = ddf+2*H*S*(S - diag(D));ddf = ddf+2*Y*dS*(S - diag(D));ddf = ddf+2*Y*S*(dS - diag(dD));endfunction YO = guess()% GUESS provides the initial starting guess Y for an F minimization.% Y = GUESS% Y will satisfy Y'*Y = I% role objective function, produces an initial guess at a minimizer% ofF.global FParameters;208Ss = FParameters.Ss; toQ = 0;% take an average% av = Ss(:,:,l);% for k=2:size(Ss,3)% av = av +Ss(:,:,k);% end% [YO,D] = eig(av);% using random guess% YO = randn(size(av));% using the identity 20YO = eye(size(Ss(:,:,1)))+ le- 1*randn(size(Sz(:,:,1)));function parameters(varargin)% PARAMETERS initializes the parameters for an instance of a% INDSCAL minimization problem.% PARAMETERS(S1,S2,...,Sn)% Si are expected to be symmetric matrices of equal size.% role sets up the global parameters used at all levels of computation.global FParameters; 10FParameters = U;FParameters.Ss = cat(3,varargin{:});function [A,B] = exact()A = diag(1:10);B = diag((1:10).^2);function [A,B] = noisy()[A,B] = exact;randn('state',0);209A = A + le-4*randn(size(A));B = B + le-4*randn(size(B));210Bibliography[1] T.A. Arias A. Edelman and S.T. Smith. The geometry of algorithms with orthogonalityconstraints. SIAM Journal on Matrix Analysis and Applications, to appear.[2] I. Andersson. Experiments with the conjugate gradient algorithm for the determinationof eigenvalues of symmetric matrices. Technical Report UMINF-4.71, University ofUmea, Sweden, 1971.[3] T. A. Arias. Multiresolution analysis of electronic structure: semicardinal and orthog-onal wavelet bases. Submitted to Reviews of Modern Physics, 1997.[4] T. A. Arias, A. Edelman, and S. T. Smith. Curvature in conjugate gradient eigenvaluecomputation with applications to materials and chemistry calculations. In J. G. Lewis,editor, Proceedings of the 1994 SIAM Applied Linear Algebra Conference, pages 233-238, Philadephia, 1994. SIAM.[5] T.A. Arias, K.J. Cho, P. Lam, and M.P. Teter. Wavelet transform representation ofthe electronic structure of materials. In Rajiv K. Kalia and Priya Vashishta, editors,Proceeedings of the '94 Mardi Gras Conference: Toward Teraflop Computing and NewGrand Challenge Applications, page 23, Commack, New York, 1995. Nova Science Pub-lishers.[6] M. Arioli, I. S. Duff, and D. Ruiz. Stopping criteria for iterative solvers. Report RAL-91-057, Central Computing Center, Rutherford Appleton Laboratory, 1992.211[7] S. Bertoluzza and G. Naldi. A wavelet collocation method for the numerical solutionof partial differential equations. Applied Computational and Harmonic Analysis, 3:1-9,1996.[8] G. Beylkin, R.R. Coifman, and V. Rokhlin. Fast wavelet transforms and numericalalgorithms. Commun. Pure and Appl. Math, 44:141-183, 1991.[9] G. Beylkin and J. M. Keiser. On the adaptive numerical solution of nonlinear partialdifferential equations in wavelet bases. PAM Report, 262, 1995.[10] G. Beylkin and J. M. Keiser. On the adaptive numerical solution of nonlinear partialdifferential equations in wavelet bases. Journal of Computational Physics, 132:233-259,1997.[11] G. Beylkin and N. Saito. Wavelets, their autocorrelation functions, and multiresolutionrepresentation of signals. Expanded abstract in Proceedings ICASSP-92, 4:381-384, 1992.[12] G. Beylkin and N. Saito. Multiresolution representations using the auto-correlationfunctions of compactly supported wavelets. IEEE Transactions on Signal Processing,41:3584-3590, 1993.[13] J. W. Bruce and P. J. Giblin. Curves and Singularities (Second Edition). CambridgeUniversity Press, Cambridge, England, 1992.[14] I. Chavel. Riemannian Geometry-A Modern Introduction. The Cambridge UniversityPress, 1993.[15] K. Cho, T.A. Arias, J.D. Joannopoulos, and Pui K. Lam. Wavelets in electronic struc-ture calculations. Physical Review Letters, 71:1808, 1993.[16] S.-N. Chow, C. Li, and D. Wang. Normal Forms and Bifurcation of Planar VectorFields. Cambridge, University Press, 1994.[17] A. Cohen and R. Danchin. Multiscale approximation of vortex patches. preprint, 1997.212[18] S. Dahlke, W. Dahmen, R. Hochmuth, and R. Schneider. Stable multiscale bases andlocal error estimation for elliptic problems. Applied Numerical Mathematics, 23(1):21-47, 1997.[19] I. Daubechies. Orthonormal bases of compactly supported wavelets. J. Comm. PureAppl. Math., 41:909-996, Oct 1988.[20] I. Daubechies. Ten Lectures on Wavelets. Number 61 in CBMS/NSF Series in AppliedMath. SIAM, 1992.[21] J. de Leeuw and W. Heiser. Theory of multidimensional scaling. In P. R. Krishnaiahand L. N. Kanal, editors, Handbook of Statistics, Vol 2, pages 285-316. North-HollandPublishing Co., 1982.[22] J. W. Demmel. A Numerical Analyst's Jordan Canonical Form. PhD thesis, Univ. ofCalifornia at Berkeley, Berkeley, CA, USA, 1983.[23] J. W. Demmel. Applied Numerical Linear Algebra. SIAM, 1997.[24] G. Deslauriers and S. Dubuc. Symmetric iterative interpolation processes.Approx., 5(1):49-68, 1989.[25] M. P. do Carmo. Differential geometry of curves and surfaces. Prentice-Hall, 1976.[26] D. L. Donoho. Interpolating wavelet tranforms. Preprint, Department of Statistics,Stanford University, 1992.[27] S. Dubuc. Interpolation through and iterative scheme. J. Math. Anal. Appl., 114:185-204, 1986.[28] A. Edelman and E. Elmroth. Personal communication, January 1997.[29] A. Edelman and Y. Ma. Non-generic eigenvalue perturbations of jordan blocks. LinearAlgebra Appl., 273:45-63, 1998.[30] A. Edelman and Y. Ma. Staircase failures explains by orthogonal versal forms. Submittedto SIAM J. Matrix Ana. and Appl., 1998.213Constr.101=19 1. -"IIII lr 101 .1[31] A. Edelman and S. T. Smith. On conjugate gradient-like methods for eigen-like prob-lems. BIT, to appear. See also Proc. Linear and Nonlinear Conjugate Gradient-RelatedMethods, eds. Loyce Adams and J. L. Nazareth. SIAM, Philadelphia, PA, 1996.[32] L. Elden. Algorithms for the regularization of ill-conditioned least-squares problems.BIT, 17:134-145, 1977.[33] T. F. Fairgrieve. The application of singularity theory to the computation of JordanCanonical Form. Master Thesis at Computer Science in Univ. of Toronto, 1986.[34] R. Fletcher. Practical Methods of Optimization. John Wiley & Sons, New York, 2dedition, 1987.[35] J. Frohlich and K. Schneider. An adaptive wavelet galerkin algorithm for one- andtwo-dimensional flame. Computations. Eur. J. Mech, 13:439-471, 1994.[36] Z. Fu and E. M. Dowling. Conjugate gradient eigenstructure tracking for adaptivespectral estimation. IEEE Trans. Signal Processing, 43(5):1151-1160, 1995.[37] P. E. Gill, W. Murray, and M. H. Wright. Practical Optimization. Academic Press,New York, 2d edition, 1981.[38] S. Goedecker and O.V. Ivanov. Linear scaling solution of the coulomb problem usingwavelets. Preprint (see http://xxx.lanl. gov/abs/physics/9701024), Submitted Jan1997.[39] G. H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins UniversityPress, Baltimore, Maryland, 2d edition, 1989.[40] G. H. Golub and J. H. Wilkinson. Ill-conditioned eigensystems and the computation ofthe Jordan canonical form. SIAM Review, 18:578-619, 1976.[41] S. Helgason. Differential Geometry, Lie Groups, and Symmetric Spaces. AcademicPress, New York, 1978.214[42] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Physical Review, 136:B864-B871, 1964.[43] T. Kato. Perturbation Theory for Linear Operators. Springer-Verlag, New York, NY,USA, 1966.[44] S. Kobayashi and K. Nomizu. Foundations of Differential Geometry, volume 1 and 2.Wiley, New York, 1969.[45] W. Kohn and L.J. Sham. Self-consistent equations including exchange and correlationeffects. Physical Review, 140:A1133, 1965.[46] J. Liandrat, V. Pierrier, and P. Tchamitchian. Numerical resolution of the regularizedBurgers equation using the wavelet transform. Tech. Report CPT-89/P. 2320 - Centerof Theoretical Physics, Marseille, France, 1989.[47] J. Liandrat and P. Tchamitchian. Resolution of the ld regularized burgers equationusing a spatial wavelet approximation. report, Technical Report 90-83, NASA ICASE,December 1990.[48] C. F. Van Loan. On estimating the condition of eigenvalues and eigenvectors. LinearAlgebra Appl., 88/89:715-732, 1987.[49] A. N. Malyshev. On wilkinson's problem. Technical report, Report 140, Dept. ofInformatics, University of Bergen, Norway, November 1997.[50] S. G. Nash and Ariela Sofer. Linear and Nonlinear Programming. McGraw-Hill, NewYork, 1995.[51] M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, and J.D. Joannopoulos. Iterativeminimization techniques for ab initio total energy calculations: molecular dynamicsand conjugate gradients. Reviews of Modern Physics, 64:1045, 1992.[52] E. Polak. Computational Methods in Optimization. Academic Press, New York, 1971.215[53] N. Saito and G. Beylkin. Multiresolution representations. Applied and ComputationalHarmonic Analysis, 41:3584, 1993.[54] A. H. Sameh and J. A. Wisniewski. A trace minimization algorithm for the generalizedeigenvalue problem. SIAM J. Numerical Analysis, 19:1243-1259, 1982.[55] W. H. A. Schilders. Personal communication, September 1997.[56] S. T. Smith. Optimization techniques on Riemannian manifolds. Fields Institute Com-munications, 3:113-146, 1994.[57] I. Stich, R. Car, M. Parrinello, and S. Baroni. Conjugate gradient minimization of theenergy functional: A new method for electronic structure calculation. Phys. Rev. B.,39:4997-5004, 1989.[58] G. Strang and T. Nguyen. Wavelets and Filter Banks. Wellesley-Cambridge Press,Cambridge, Massachusetts, 1996.[59] S. H. Strogatz. Nonlinear dynamics and Chaos: with applications to physics, biology,chemistry, and engineering. Addison-Wesley, 1994.[60] J.-G. Sun. A note on simple non-zero singular values. Journal of Computational Math-ematics, 6:258-266, 1988.[61] J.-G. Sun. Sensitivity analysis of multiple eigenvalues. I. Journal of ComputationalMathematics, 6:28-38, 1988. Cited in [62].[62] J.-G. Sun. A note on local behavior of multiple eigenvalues. SIAM Journal on MatrixAnalysis and Applications, 10:533-541, 1989.[63] L. N. Trefethen. Pseudospectra of matrices. In D. F. Griffiths and G. A. Watson, editors,Numerical Analysis 1991. Longman Scientific and Technical, Harlow, Essex, 1992.[64] C.J. Tymczak and X.-Q. Wang. Orthonormal wavelet bases for quantum moleculardynamics. Physical Review Letters, 78:3654, 1997.216[65] S. Wei and M.Y Chou. Wavelets in self-consistent electronic structure calculations.Physical Review Letters, 76:2650, 1996.[66] J. H. Wilkinson. The Algebraic Eigenvalue Problem. Clarendon Press, Oxford, England,1965.[67] J. H. Wilkinson. Sensitivity of eigenvalues. Util. Math., 25:5-76, 1984.[68] J. H. Wilkinson. Sensitivity of eigenvalues II. Util. Math., 30:243-286, 1986.[69] Y.-C. Wong. Differential geometry of Grassmann manifolds. Proc. Natl. Acad. Sci.USA, 57:589-594, 1967.[70] D. Yegilleten and T. A. Arias. Synthesis of wavelet theory with the multigrid algorithm.Condensed Matter Physics, MIT, 1996.217

Nonlinear eigenvalue problems

Abstract

Similar works

Full text

Available Versions

DSpace@MIT