A computationally efficient approach for the removal of the phase shift singularity in harmonic resolvent analysis
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b A computationally efficient approach for the removalof the phase shift singularity in harmonic resolventanalysis
Alberto Padovan · Clarence W. Rowley
Received: date / Accepted: date
Abstract
The recently introduced harmonic resolvent framework [1] is concernedwith the study of the input-output dynamics of nonlinear flows in the proximity ofa known time-periodic orbit. These dynamics are governed by the harmonic resol-vent operator, which is a linear operator in the frequency domain whose singularvalue decomposition sheds light on the dominant input-output structures of theflow. Although the harmonic resolvent is a mathematically well-defined operator,the numerical computation of its singular value decomposition requires inverting amatrix that becomes exactly singular as the periodic orbit approaches an exact so-lution of the nonlinear governing equations. The very poor condition properties ofthis matrix hinder the convergence of classical Krylov solvers, even in the presenceof preconditioners, thereby increasing the computational cost required to performthe harmonic resolvent analysis. In this paper we show that a suitable augmen-tation of the (nearly) singular matrix removes the singularity, and we provide alower bound for the smallest singular value of the augmented matrix. We also showthat the desired decomposition of the harmonic resolvent can be computed usingthe augmented matrix, whose improved condition properties lead to a significantspeedup in the convergence of classical iterative solvers. We demonstrate this sim-ple, yet effective, computational procedure on the Kuramoto-Sivashinsky equationin the proximity of an unstable time-periodic orbit.
Keywords
Harmonic resolvent · Time-periodic systems · Singular linear systems
Throughout the years some fluid flows of interest have been studied using a split-ting approach, whereby the nonlinear terms in the Navier-Stokes equation are
Alberto PadovanDepartment of Mechanical and Aerospace Engineering, Princeton University, NJ 08544, USAE-mail: [email protected] W. RowleyDepartment of Mechanical and Aerospace Engineering, Princeton University, NJ 08544, USAE-mail: [email protected] Alberto Padovan, Clarence W. Rowley treated as forcing that acts on the linear terms. This approach is particularly con-venient since the linear dynamics lend themselves to classical linear analyses thatcan help uncover some of the fundamental mechanisms behind complicated phys-ical phenomena. For instance, the input-output analysis of the linearized Navier-Stokes operator about a steady base flow has helped shed light on some of theenergy amplification mechanisms in shear flows of interest. In particular, resolventanalysis was used to study the response of perturbations to spatio-temporal forcingin linearized channel flow [2]. A similar approach was implemented in [3] to inden-tify the most amplified velocity structures at selected temporal frequency-spatialwavenumber pairs of interest in turbulent pipe flow.Recently, this linear input-output framework was extended to analyze the dom-inant dynamics of perturbations about periodically time-varying base flows [1].This extended framework, known as harmonic resolvent analysis , is based on thesingular value decomposition of the harmonic resolvent operator, a frequency-domain linear input-output operator that governs the dynamics of time-periodicperturbations about a time-periodic base flow. Much like the singular value de-composition of the resolvent operator discussed in [2, 3], the singular value decom-position of the harmonic resolvent provides insight into the dominant input-outputstructures of the flow.While the harmonic resolvent operator is mathematically well-defined, the nu-merical computation of its singular values requires some care. In particular, we willsee in section 2.2 that computing the desired decomposition of the harmonic resol-vent requires inverting the linearized Navier-Stokes operator evaluated about theperiodic base flow. It is well known, however, that linearized periodic dynamics areneutrally stable in the direction of a phase shift along the periodic orbit [4]. Thelinearized Navier-Stokes operator will thus have a one-dimensional nullspace alongthe direction of the phase shift, and this singularity hinders the performance ofclassical Krylov-based solvers. This is especially problematic in large-scale applica-tions, where iterative solvers may be the only computationally feasible algorithmsfor the solution of linear systems.We will show in section 3 that a suitable augmentation of the linearized Navier-Stokes operator allows for the removal of the singularity, and we provide a lowerbound for the smallest singular value of the augmented operator. We then showthat the desired singular value decomposition of the harmonic resolvent operatorcan be obtained by working with the better conditioned augmented matrix. Insection 4 we apply this computational procedure to the Kuramoto-Sivashinskyequation in the proximity of an unstable periodic orbit, and we demonstrate thespeedup in convergence that is obtained by properly removing the singularity.Before moving forward, it is worth mentioning that the procedure we proposeshares some similarities with algorithms for the solution of singular linear systems,available in linear algebra packages such as PETSc [5]. These usually rely on theknowledge of the nullspace of the linear operator and of its complex conjugatetranspose to compute the least squares solution for the linear system at hand. emoval of the singularity in the harmonic resolvent analysis 3 q ( t ), and dynamics given bydd t q ( t ) = f (cid:0) q ( t ) (cid:1) . (1)We then decompose the state about a periodic base flow Q ( t ) with period T : q ( t ) = Q ( t ) + q ′ ( t ) = X ω ∈ Ω b ˆ Q ω e iωt + X ω ∈ Ω ˆ q ′ ω e iωt , (2)where Ω b ⊆ Ω ⊂ πT Z . Here, Ω b is the set of frequencies associated with the baseflow, while Ω is the set of frequencies associated with the perturbations that wewish to resolve. Upon substituting (2) into (1) one obtainsdd t q ′ ( t ) = D q f ( Q ( t )) | {z } A ( t ) q ′ ( t ) + h ′ ( t ) (3)where h ′ ( t ) contains higher-order terms. Formula (3) can be written in the fre-quency domain as iω ˆ q ′ ω = X α ∈ Ω ˆ A ω − α ˆ q ′ α + ˆ h ′ ω ∀ ω ∈ Ω. (4)For ease of notation, let ˆ q ′ be the vector of ˆ q ′ ω for all frequencies ω ∈ Ω , and letˆ h ′ be defined similarly. We then define the operator T by (cid:2) T ˆ q ′ (cid:3) ω = iω ˆ q ′ ω − X α ∈ Ω ˆ A ω − α ˆ q ′ α . (5)If the base flow Q ( t ) is an exact solution of (1), then it is straightforward to seethat T is singular. In particular, Q ( t + τ ) is also an exact solution, for any shift τ ,and by differentiating with respect to τ , one can show that (d / d t ) Q ( t ) exactlysatisfies (3) with h ′ = 0, and thus the Fourier coefficients of d Q / d t lie in thenullspace of T . That is, the nullspace of T is along the direction of a phase shiftof the base flow.On the other hand, if Q ( t ) is an approximate solution of (1), then T is nearly singular along the direction of the phase shift. As we are not usually interestedin the trivial phase shift, the harmonic resolvent was defined in [1] in a way thatremoves it. Specifically, letting v be the unit-norm vector in the direction of phaseshift given by d Q / d t V , and letting Σ = v ⊥ denote its orthogonal complement, wedefine the restricted operator T | Σ : Σ → W Σ , (6)where W Σ (the range of T | Σ ) is a codimension-1 subspace orthogonal to a unit-norm vector u . We will further discuss u and its efficient computation in the Alberto Padovan, Clarence W. Rowley upcoming sections. Notice that restricting the range and domain of the operator T is analogous to constructing a Poincar´e map by reducing the dynamics ontoa codimension-1 subspace pierced by the limit cycle [4]. Finally, the harmonicresolvent on W Σ is defined as H = ( T | Σ ) − . (7)Upon the removal of the phase shift direction, formula (4) may be written as T ˆ q ′ = ˆ h ′ ⇐⇒ ˆ q ′ = H ˆ h ′ . (8)2.2 Computational procedure and challengesAs mentioned in the introduction, we are interested in computing the singularvalue decomposition (SVD) of the harmonic resolvent in order to shed light onthe dominant input-output structures of the flow in the proximity of the time-periodic base flow. However, given the high dimensionality of systems arising fromthe discretization of partial differential equations, H cannot be computed andstored explicitly. Instead, given the sparse operator T , one can use a randomizedSVD algorithm to compute the leading singular values and singular vectors of theharmonic resolvent. We refer the reader to [6] for a detailed description of the algo-rithm. For the upcoming discussion, it suffices to point out that the computationof the randomized SVD requires evaluating matrix-vector products of the form H ˆ h ′ and H ∗ ˆ q ′ , where H ∗ is the complex conjugate transpose of H . Specifically,in order to evaluate H ˆ h ′ we solve the linear system T ˆ q ′ = (cid:0) I − uu ∗ (cid:1) ˆ h ′ (9)where ( I − uu ∗ ) is the orthogonal projection onto W Σ = u ⊥ . Likewise, to evaluate H ∗ ˆ q ′ we solve the linear system T ∗ ˆ w ′ = ˆ q ′ , ˆ h ′ = (cid:0) I − uu ∗ (cid:1) ˆ w ′ . (10)Since T is often singular or poorly conditioned, solving (9) or (10) may be problem-atic. If T is exactly singular, both direct solution algorithms and iterative solverswill fail. If T is very poorly conditioned, exact solvers based on matrix decomposi-tions (e.g., LU decomposition) may be a viable option in small- and moderate-sizedproblems. As the size of the problem increases, the cost associated with perform-ing a matrix decomposition grows polynomially and exact solution algorithms maybecome inaccessible. This limit can be quickly approached in two-dimensional fluidflows, where T may have dimensions on the order O (cid:0) (cid:1) − O (cid:0) (cid:1) . Finally, it-erative solvers may suffer even in the presence of carefully chosen preconditionersif T is poorly conditioned. Therefore, if the (near) singularity in T could be re-moved, this could significantly reduce the computational cost required to performthe harmonic resolvent analysis. emoval of the singularity in the harmonic resolvent analysis 5 We now present a computationally efficient way to remove the singularity from T .Throughout this section, we will work with an augmented linear operator e T = (cid:20) T uv ∗ (cid:21) , (11)where u and v are defined as in the previous section: that is, v is the unit-normvector in the direction of the phase shift, and u is a unit-norm vector orthogonalto the range of T | Σ (where Σ = v ⊥ ). While v is easily computed from the timederivative of the base flow, details on the computation of u will be provided insection 3.3. We note in passing that operators of the form of (11) arise in Newton-based harmonic balancing methods for the solution of nonlinear systems, where T would be the Jacobian matrix at the k th iteration and v a phase constrainton the k th update. In presenting the main results of this paper, we consider twoscenarios: when T is exactly singular and when T is nearly singular.3.1 Singular T We first consider the case when T is exactly singular, with its one-dimensionalnullspace spanned by the direction of phase shift about the base flow. This scenariois likely to arise when the base flow is computed with accuracy close to machineprecision using harmonic balancing methods. These methods are ubiquitous inmost branches of physics and a similar approach has recently been adopted in[7] to compute a time-periodic and spanwise-periodic solution for the transitionto turbulence in a forced boundary layer. The main result of this subsection ispresented in the proposition below, which states that augmenting T as in (11)removes the singularity, and that the harmonic resolvent operator can be definedin terms of the augmented matrix. Proposition 1
Consider the singular value decomposition of T ∈ C N × N given by T = N X j =1 σ j u j v ∗ j , (12) where σ N = 0 , and σ j > for all j < N . Let e T ∈ C ( N +1) × ( N +1) be defined as in(11), with v = v N and u = u N . Then the following hold:1. e T is invertible and its singular value decomposition is given by e T = N − X j =1 σ j e u j e v ∗ j + N +1 X j = N e u j e v ∗ j (13) where e v j = (cid:20) v j (cid:21) , e u j = (cid:20) u j (cid:21) , j ∈ { , , · · · , N − } (14) and e v N = (cid:20) v N (cid:21) , e u N = (cid:20) (cid:21) , e v N +1 = (cid:20) (cid:21) , e u N +1 = (cid:20) u N (cid:21) . (15) Alberto Padovan, Clarence W. Rowley
2. The harmonic resolvent, defined in (7), is given by H = N − X j =1 σ j v j u ∗ j . (16) Therefore, the singular values and singular vectors of H can be found from theSVD of e T , according to (14). The first part of the proposition states that if v and u in (11) are properlychosen, then the augmented operator is invertible and well-conditioned. More pre-cisely, the non-zero singular values of T agree with N − e T , and the zero singular value σ N of T is replaced by two singular values withvalue one. The second statement says that one can easily compute the singularvalues and singular vectors of the harmonic resolvent from the SVD of the better-conditioned augmented matrix e T .In order to prove the first statement it suffices to check that e T e v j = σ j e u j , ∀ j ∈ { , , · · · , N + 1 } . The proof of the second statement follows immediatelyfrom (12) and the definition of the harmonic resolvent in (7).3.2 Nearly singular T We now consider the case when T is not exactly singular. This is usually thecase when the base flow is computed via numerical integration of the governingequations, as in [1], and numerical errors and truncation errors slightly perturbthe null singular value. The perturbed matrix is then invertible, but it may stillbe poorly conditioned because of this small singular value, and it may thereforebecome necessary to remove the near singularity in order to improve the perfor-mance of iterative solvers. As in the previous section, we would like to show thatthe augmented matrix e T is better conditioned than T , and that the SVD of theharmonic resolvent operator can be computed using e T .We start by showing that the SVD of the harmonic resolvent can be computedusing e T . First and foremost, we recall that the SVD of the harmonic resolventis performed numerically by computing matrix-vector products of the form H ˆ h ′ and H ∗ ˆ q ′ . As mentioned in section 2.2, these are usually computed using T , bysolving the linear systems (9) and (10). The proposition below states that we cancompute these matrix-vector products using the better conditioned matrix e T . Proposition 2
Let e T , u , and v be defined as in (11), with T full rank. Then thefollowing two statements hold:1. If H denotes the harmonic resolvent defined by (7), then ˆ q ′ = H ˆ h ′ solveseither of the following systems: (cid:20) T uv ∗ (cid:21) (cid:20) ˆ q ′ λ (cid:21) = (cid:20) ˆ h ′ (cid:21) ⇐⇒ T ˆ q ′ = (cid:0) I − uu ∗ (cid:1) ˆ h ′ . (17)
2. If H ∗ denotes the adjoint of H , then ˆ h ′ = H ∗ ˆ q ′ solves either of the followingsystems: (cid:20) T ∗ vu ∗ (cid:21) (cid:20) ˆ h ′ λ (cid:21) = (cid:20) ˆ q ′ (cid:21) ⇐⇒ T ∗ ˆ w ′ = ˆ q ′ , ˆ h ′ = (cid:0) I − uu ∗ (cid:1) ˆ w ′ . (18) emoval of the singularity in the harmonic resolvent analysis 7 This means that the action of H or H ∗ on a vector (and hence the singularvalue decomposition of H ) may be computed using either T or e T . The proof ofthis proposition is given in the appendix. The next proposition establishes thatthe augmented matrix e T is better conditioned than the original matrix T , and weprovide a lower bound for the smallest singular value of e T . Proposition 3
Let σ be the smallest singular value of the operator T | Σ definedin (6), and suppose that k T v k = ε < . Then we have (cid:13)(cid:13)(cid:13) e T ˆ z ′ (cid:13)(cid:13)(cid:13) ≥ γ (cid:13)(cid:13) ˆ z ′ (cid:13)(cid:13) , ∀ ˆ z ′ (19) where γ = min { − ε, σ (1 − ε ) / } . (20) That is, the minimum singular value of e T is at least γ . In most examples we have encountered, the smallest singular value of T isin the direction of the phase shift; that is, ε is smaller than σ . In this case, theproposition states that the smallest singular value of e T is either close to one, orclose to σ . Thus, e T is better conditioned than T since we have removed its smallestsingular value. The proof of this proposition is also available in the appendix.3.3 Computing the vector u In the previous subsections, we have defined an augmented operator e T relying onthe knowledge of the appropriate vectors v and u that would remove the singularityfrom T . It is straightforward to compute the vector v , the unit-norm vector in thedirection of the phase shift, which is given by the Fourier coefficients of the timederivative of the base flow. However, computing u requires some care.Recall that u is the orthogonal complement of the range of T | Σ , where Σ = v ⊥ ,and consider the system (cid:20) T ∗ vw ∗ (cid:21) (cid:20) ˆ z ′ λ (cid:21) = (cid:20) (cid:21) , (21)where w is an arbitrarily chosen vector. We readily see that the system above canbe written as T ∗ ˆ z ′ = − λ v (22) w ∗ ˆ z ′ = 1 . (23)Any vector ˆ z ′ that satisfies the above equations must be orthogonal to the rangeof T | Σ , since for any ˆ q ′ that lies in Σ = v ⊥ we have h T ˆ q ′ , ˆ z ′ i = h ˆ q ′ , T ∗ ˆ z ′ i = h ˆ q ′ , − λ v i = 0 . (24)We therefore solve (21) for ˆ z ′ and set u = ˆ z ′ / k ˆ z ′ k . This approach may fail if wemistakenly choose w to lie entirely in the range of T | Σ , in which case w ∗ ˆ z ′ = 0and the system would have no solution. This risk, however small it may be, can beavoided by generating w in such a way that it is close to the orthogonal complement Alberto Padovan, Clarence W. Rowley of the range of T | Σ . For instance, letting T D be the block diagonal componentsof T , one can do so by letting w be the solution of the system T ∗ D w = v , (25)which can be understood as the block-Jacobi solution of (22) with λ = −
1. Evenif T is exactly singular, T D is full-rank and the reader may recognize in thisblock-diagonal operator the linear operator that governs the linearized dynamicsabout the temporal mean (i.e., its inverse would contain the well-known resolventoperators discussed in [2, 3]). We now illustrate the main results on the Kuramoto-Sivashinsky equation, whichis a one-dimensional partial differential equation that arises in the description ofinstabilities on interfaces and flame fronts. This equation was chosen for three rea-sons. First, it exhibits complex spatio-temporal dynamics similar to those that canarise from the Navier-Stokes equation [8]. Second, its spatio-temporal discretiza-tion is low-dimensional enough that we can compute the entire spectrum of theoperators T and e T , thereby providing empirical evidence for the theoretical re-sults developed in sections 3.1 and 3.2. Third, its spatio-temporal discretizationis high-dimensional enough that we can demonstrate the faster convergence ofKrylov solvers when the phase-shift singularity is removed.We consider the equation in the form ∂u∂t = − u ∂u∂x − ∂ u∂x − ∂ u∂x , (26)where the state u ( x, t ) is defined over the periodic spatial domain X = [ − L/ , L/ L = 39, for which thereexist a chaotic attractor and a number of unstable time-periodic orbits [9]. Wespecify, for the sake of completeness, that the spatial discretization is performedusing a Fourier-spectral method, and we retained 32 spatial wavenumbers.We henceforth omit the spatial dependence of the state variable for notationalsimplicity. Given a periodic orbit, denoted U ( t ), we linearize the dynamics byperforming the following expansion of the state u ( t ) = U ( t ) + u ′ ( t ) = X ω ∈ Ω b ˆ U ω e iωt + X ω ∈ Ω ˆ u ′ ω e iωt , (27)where u ′ ( t ) are time periodic perturbations about the periodic orbit. As discussedin section 2.1, Ω b is the set of frequencies associated with the periodic base flow,while Ω is the set of frequencies associated with the perturbations. We henceforthtake Ω = {− , − , · · · , } ω f , where ω f is the fundamental frequency of os-cillation. Upon substituting (27) into the nonlinear dynamics given by (26), weobtain iω ˆ u ′ ω = − (cid:18) ∂ ∂x + ∂ ∂x (cid:19) ˆ u ′ ω − X α ∈ Ω " ˆ U ω − α ∂∂x + ∂ ˆ U ω − α ∂x ˆ u ′ α +ˆ h ′ ω , ∀ ω ∈ Ω (28) emoval of the singularity in the harmonic resolvent analysis 920406080100120 t −
10 0 10 x − . − . . . U ( x , t ) (a) − − − − h ˆ U k ω f , ˆ U k ω f i × ω f kω f (b) Fig. 1: We show (a) one of the unstable periodic orbits, with period T = 2 π/ω f ≈ .
9, computed using an harmonic balancing approach, and (b) its energy spec-trum.where ˆ h ′ ω contains all the nonlinear terms at frequency ω . Letting ˆ u ′ denote thecollection of all Fourier modes, formula (28) can be written compactly as T ˆ u ′ = ˆ h ′ . (29)We can now verify the theoretical results obtained in section 3. Throughout theremainder of this section, the augmented operator e T is defined as in (11), where v is the unit norm vector in the direction of a phase shift given by the time derivativeof U ( t ), and u is computed following the procedure described in section 3.3.4.1 Singular T As discussed in the previous sections, when the base flow U ( t ) satisfies the dynam-ics exactly, the matrix T will be singular with a one-dimensional nullspace in thedirection of the phase shift about the base flow. We compute one of the unstable pe-riodic orbits that exist for the chosen configuration of the Kuramoto-Sivashinskyequation using a harmonic balancing method. In the Fourier expansion of thecandidate solution we consider 24 harmonics of the (yet unknown) fundamentalfrequency, and we therefore let Ω b = {− , − , · · · , } ω f . The spatio-temporalevolution and the energy spectrum of this orbit are shown in figures 1a and 1b,and the fundamental period of oscillation is found to be T = 2 π/ω f ≈ . T and e T . Specif-ically, from figure 2a, we observe that while T has a zero singular value ( σ N ≈ . × − ), e T does not, and its smallest singular value agrees with the smallestnon-zero singular value of T . Furthermore, we observe from figure 2b that two sin-gular values of e T have value one. Finally, we notice that except for the highlighteddifferences, the singular values of T and e T agree exactly. These observations areas we expect from P proposition 1. . . . . . .
998 Sing. values of T and e T near zeroSing. values of T and e T near unity σ N (a)(b) − − − − k R e s i du a l k (c) Fig. 2: We show (a) the singular values of T and e T close to zero (the smallestsingular value of T is σ N ∼ O (10 − ), indicating that T is singular for practicalpurposes), (b) the singular values of T and e T close to unity and (c) a convergenceplot obtained by solving the linear systems in (30) and (31). The dashed line in (b) is an extension of the 1-tick on the y axis. ( ◦ : data for T ; × : data for e T ).Figure 2c shows a representative convergence plot illustrating the performanceof the block-Jacobi-preconditioned GMRES solver on the linear systems T x = h (30) e T e x = ( h , , (31)where h is a random unit-norm vector. The benefits of removing the singularityfrom the linear operator T are clear. Specifically, we see that the solver usingequation (30) (blue curve) plateaus at a residual on the order of 10 − , while thesolver using equation (31) (red curve) converges to the desired tolerance in fewerthan 200 iterations.4.2 Nearly singular T We now consider the case when the base flow U ( t ) does not satisfy the govern-ing equation exactly. We introduce an error in the base flow described in theprevious section by truncating its highest frequency component, so that Ω b = {− , − , · · · , } ω f . The set of frequencies Ω associated with the perturbedstate u ′ ( t ) is kept unchanged.As expected, T is now non-singular, and its smallest singular value (shown infigure 3a) has order of magnitude O (10 − ). Constructing e T is beneficial nonethe-less, and we see from figure 3a that the smallest singular value is removed, andit is replaced with two singular values with value ≈ e T do not agreeexactly with those of T . We recall, however, that we are ultimately interested inthe singular values and singular vectors of the harmonic resolvent H , and we haveshown via proposition 2 that we are allowed to use e T to compute its singular valuedecomposition. emoval of the singularity in the harmonic resolvent analysis 1100 . . . . . .
998 Sing. values of T and e T near zeroSing. values of T and e T near unity σ N (a)(b) − − − − k R e s i du a l k (c) Fig. 3: Analog of figure 2, showing (a) the singular values of T and e T close to zero(the smallest singular value of T is σ N ∼ O (10 − )), (b) the singular values of T and e T close to unity and (c) a convergence plot obtained by solving the linearsystems in (30) and (31). The dashed line in (b) is an extension of the 1-tick onthe y axis. ( ◦ : data for T ; × : data for e T ).Finally, figure 3c shows a representative convergence plot illustrating the per-formance of the block-Jacobi-preconditioned GMRES solver on the linear systemsin equations (30) and (31). Although the speed-up in convergence is not as sub-stantial as the one shown in figure 2c, we observe that removing the smallestsingular value still leads to slightly faster convergence. It is well-known that linear time-periodic dynamics are neutrally stable in the di-rection of a phase shift about the time-periodic orbit, and this property manifestsitself in the form of a singularity in the linear operator that governs the dynam-ics. This singularity, in turn, leads to numerical difficulties in the context of theharmonic resolvent analysis, where the computation of the singular value decom-position of the harmonic resolvent operator requires the inversion of this singularoperator.We have proposed a computationally inexpensive solution to this problem,showing that a suitable augmentation of the singular matrix leads to the removalof the singularity and to a significant improvement in the condition propertiesof the resulting augmented linear operator. In our discussion, we have consid-ered the cases when the operator is exactly singular and when the operator isnearly singular. We then used the Kuramoto-Sivashinsky equation as an exampleto demonstrate that in both cases it is convenient to remove the (near) singular-ity in order to improve the convergence properties of classical Krylov-based linearsolvers.
Acknowledgements
The authors wish to thank Samuel Otto for providing useful comments that helpedimprove the presentation of this paper. This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-17-1-0084.
Conflict of interest
The authors declare that they have no conflict of interest.
A Proof of proposition 2
We start by verifying the first statement, which we copy here for clarity. (cid:20)
T uv ∗ (cid:21) (cid:20) ˆ q ′ λ (cid:21) = (cid:20) ˆ h ′ (cid:21) ⇐⇒ T ˆ q ′ = ( I − uu ∗ ) ˆ h ′ . (32)We first recall that v is the unit-norm orthogonal complement of Σ , while u is the unit-normorthogonal complement of the range of T | Σ . From the second row of the augmented linearsystem we see that h ˆ q ′ , v i = 0, which implies that ˆ q ′ ∈ Σ , and hence h T ˆ q ′ , u i = 0. We proceedby taking the inner product of the first row of the augmented system with u , to obtain h T ˆ q ′ , u i + λ h u , u i = h ˆ h ′ , u i ⇒ λ = h ˆ h ′ , u i , (33)where we have used h u , u i = 1. Substituting λ into the first row of the augmented systemestablishes the desired result.We proceed by verifying the second statement of the proposition, given by (cid:20) T ∗ vu ∗ (cid:21) (cid:20) ˆ h ′ λ (cid:21) = (cid:20) ˆ q ′ (cid:21) ⇐⇒ T ∗ ˆ w ′ = ˆ q ′ , ˆ h ′ = ( I − uu ∗ ) ˆ w ′ . (34)We first leverage the fact that T is invertible and we recall, as discussed in section 3.3, that u is given by z = ( T ∗ ) − v , u = z / k z k . (35)From the first row of the augmented system and using (35) we haveˆ h ′ = ( T ∗ ) − (cid:0) ˆ q ′ − λ v (cid:1) = ( T ∗ ) − ˆ q ′ − λ k z k u . (36)By the second row of the augmented system we have h ˆ h ′ , u i = 0, thus, using (36) we have h ( T ∗ ) − ˆ q ′ , u i − λ k z k h u , u i = 0 ⇒ λ = 1 k z k h ( T ∗ ) − ˆ q ′ , u i . (37)Substituting λ into (36) we obtainˆ h ′ = ( I − uu ∗ ) ( T ∗ ) − ˆ q ′ , (38)which is precisely what is given on the right of the equivalence sign in (34), with ˆ w ′ =( T ∗ ) − ˆ q ′ . B Proof of proposition 3
First, observe that since σ is the smallest singular value of T | Σ , it follows that k T ˆ q ′ k ≥ σ k ˆ q ′ k , for all ˆ q ′ ∈ Σ. (39)emoval of the singularity in the harmonic resolvent analysis 13Let us write ˆ z ′ = (ˆ q ′ + α v , λ ), where ˆ q ′ ∈ Σ and α, λ ∈ C . Note that any ˆ z ′ can be writtenthis way, since Σ is the orthogonal complement of v . Then (cid:13)(cid:13) ˆ z ′ (cid:13)(cid:13) = (cid:13)(cid:13) ˆ q ′ + α v (cid:13)(cid:13) + | λ | = (cid:13)(cid:13) ˆ q ′ (cid:13)(cid:13) + | α | + | λ | , since v ⊥ ˆ q ′ and k v k = 1. Furthermore, e T ˆ z ′ = (cid:20) T (ˆ q ′ + α v ) + λ uv ∗ (ˆ q ′ + α v ) (cid:21) = (cid:20) T (ˆ q ′ + α v ) + λ u α (cid:21) again since v ⊥ ˆ q ′ and k v k = 1. Now, (cid:13)(cid:13)(cid:13) e T ˆ z ′ (cid:13)(cid:13)(cid:13) = (cid:13)(cid:13) T ˆ q ′ + α T v + λ u (cid:13)(cid:13) + | α | == (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) + | α | k T v k + | λ | k u k ++ 2Re (cid:2) h α T v , T ˆ q ′ i + h λ u , T ˆ q ′ i + h α T v , λ u i (cid:3) + | α | == (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) + (1 + ε ) | α | + | λ | ++ 2Re (cid:2) h α T v , T ˆ q ′ i + h α T v , λ u i (cid:3) since u ⊥ T ˆ q ′ and k T v k = ε . Note that for any x, y , we have − Re h x, y i ≤ |h x, y i| ≤ k x k · k y k , thanks to the Cauchy-Schwarz inequality. Therefore, (cid:13)(cid:13)(cid:13) e T z (cid:13)(cid:13)(cid:13) ≥ (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) + (1 + ε ) | α | + | λ | − ε | α | (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) − ε | α | | λ | , since T v = ε . Since − ab ≥ − ( a + b ) for any a, b , we have (cid:13)(cid:13)(cid:13) e T ˆ z ′ (cid:13)(cid:13)(cid:13) ≥ (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) + (1 + ε ) | α | + | λ | − ( ε | α | + ε (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) ) − ( ε | α | + ε | λ | )= (1 − ε ) (cid:13)(cid:13) T ˆ q ′ (cid:13)(cid:13) + (1 − ε + ε ) | α | + (1 − ε ) | λ | ≥ (1 − ε ) σ (cid:13)(cid:13) ˆ q ′ (cid:13)(cid:13) + (1 − ε ) | α | + (1 − ε ) | λ | thanks to (39), since 0 < − ε <
1. So if we let γ = min { (1 − ε ) σ , (1 − ε ) , (1 − ε ) } = min { (1 − ε ) σ , (1 − ε ) } , we have (cid:13)(cid:13)(cid:13) e T ˆ z ′ (cid:13)(cid:13)(cid:13) ≥ γ ( (cid:13)(cid:13) ˆ q ′ (cid:13)(cid:13) + | α | + | λ | ) = γ (cid:13)(cid:13) ˆ z ′ (cid:13)(cid:13) . Taking square roots then establishes (19) for γ given by (20). References
1. A. Padovan, S.E. Otto, C.W. Rowley, Journal of Fluid Mechanics , A14 (2020)2. M.R. Jovanovi´c, B. Bamieh, Journal of Fluid Mechanics , 145 (2005)3. B.J. McKeon, A.S. Sharma, Journal of Fluid Mechanics , 336 (2010)4. J. Guckenheimer, P. Holmes,
Nonlinear Oscillations, Dynamical Systems, and Bifurcationsof Vector Fields (Springer, 2002)5. S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin,A. Dener, V. Eijkhout, W.D. Gropp, D. Karpeyev, D. Kaushik, M.G. Knepley, D.A.May, L.C. McInnes, R.T. Mills, T. Munson, K. Rupp, P. Sanan, B.F. Smith, S. Zampini,H. Zhang, H. Zhang, PETSc users manual. Tech. Rep. ANL-95/11 - Revision 3.12, ArgonneNational Laboratory (2019). URL
6. N. Halko, P.G. Martinsson, J.A. Tropp, SIAM Review , 217 (2011)7. G. Rigas, D. Sipp, T. Colonius, Journal of Fluid Mechanics , A15 (2021)8. P. Holmes, J.L. Lumley, G. Berkooz, C.W. Rowley, Turbulence, Coherent Structures, Dy-namical Systems and Symmetry , 2nd edn. (Cambridge, 2012)9. D. Lasagna, SIAM J. on Applied Dynamical Sysyems17