Optimal Solvers for Linear Systems with Fractional Powers of Sparse SPD Matrices
Stanislav Harizanov, Raytcho Lazarov, Pencho Marinov, Svetozar Margenov, Yavor Vutov
OOPTIMAL SOLVERS FOR LINEAR SYSTEMS WITH FRACTIONAL POWERSOF SPARSE SPD MATRICES
STANISLAV HARIZANOV, RAYTCHO LAZAROV, PENCHO MARINOV, SVETOZAR MARGENOV,AND YAVOR VUTOV
Abstract.
In this paper we consider efficient algorithms for solving the algebraic equation A α u = f , 0 < α <
1, where A is a properly scaled symmetric and positive definite matrix obtained fromfinite difference or finite element approximations of second order elliptic problems in R d , d = 1 , , u = A β − α F with F = A − β f with β positive integer. The approxi-mate solution method we propose and study is based on the best uniform rational approximationof the function t β − α for 0 < t ≤
1, and the assumption that one has at hand an efficient method(e.g. multigrid, multilevel, or other fast algorithm) for solving equations like ( A + c I ) u = F , c ≥ Introduction
Motivation for our study.
Let Ω be a bounded domain in R d , d = 1 , ,
3, with polygonalboundary Γ = ∂ Ω = ¯Γ D ∪ ¯Γ N , where Γ D has positive measure. Let q ( x ) ≥ a ( x ) ∈ R d × d be a symmetric and positive definite (SPD) matrix uniformly bounded in Ω, i.e.,(1) cξ T ξ ≤ ξ T a ( x ) ξ ≤ Cξ T ξ ∀ ξ ∈ R d , ∀ x ∈ Ω , for some positive constants c and C . Next, on V × V , V := { v ∈ H (Ω) : v ( x ) = 0 on Γ D } definethe bilinear form(2) A ( u, v ) := (cid:90) Ω (cid:0) a ( x ) ∇ u ( x ) · ∇ v ( x ) + q ( x ) u ( x ) v ( x ) (cid:1) dx. Under the assumptions on a ( x ), q , and Γ, the bilinear form is symmetric and coercive on V . Further,introduce T : L := L (Ω) → V , where for f ∈ L (Ω) the function u = T f ∈ V is the uniquesolution to A ( u, φ ) = ( f, φ ) , ∀ φ ∈ V, and ( v, u ), for u, v ∈ L (Ω) is the inner product in L (Ω).The goal of this paper is to study methods and algorithms for solving the finite element approx-imation of the operator equation(3) L α u = f, L = T − , L α u ( x ) = ∞ (cid:88) i =1 λ αi c i ψ i ( x ) , where u ( x ) = ∞ (cid:88) i =1 c i ψ i ( x ) , { ψ i ( x ) } ∞ i =1 are the eigenfunctions of L , orthonormal in L -inner product and { λ i } ∞ i =1 are the cor-responding eigenvalues that are real and positive.This definition is general, but different from the definition of the fractional powers of ellipticoperators with homogeneous Dirichlet data defined through Riesz potentials, which generalizes theconcept of equally weighted left and right Riemann-Liouville fractional derivative defined in onespace dimension to the multidimensional case, see, e.g. [3]. There is ongoing research about therelations of these two different definitions and their possible applications to problems in science andengineering, see, e.g. [2]. However, we shall focus on the current definition and withhold commentsand references on such works.Studying and numerically solving such problems is motivated by the recent development in thefractional calculus and its numerous applications to Hamiltonian chaos, anomalous transport, and Date : July 7, 2018. a r X i v : . [ m a t h . NA ] M a r S. HARIZANOV, R. LAZAROV, P. MARINOV, S. MARGENOV, Y. VUTOV super-diffusion, [31], anomalous diffusion in complex systems such as turbulent plasma, convec-tive rolls, and zonal flow system, [1], long-range interaction in elastic deformations, [25], nonlocalelectromagnetic fluid flows, [21], image processing, [10], nonlocal evolution equations arising in ma-terials science, [2]. A recent discussion about various anomalous diffusion models, their propertiesand applicability to chemistry and engineering one can find in [22]. These applications lead to var-ious types of fractional order partial differential equations that involve in general non-symmetricelliptic operators see, e.g. [17].An important subclass of such problems are the fractional powers of self-adjoint elliptic operatorsdescribed below, which are nonlocal but self adjoint. Assume that a finite element method has beenapplied to approximate the problem (3) and this resulted in a certain algebraic problem. The aimof this paper is to address the issue of solving such systems. The rigorous error analysis of suchapproximation is a difficult task that is outside the scope of this paper. Such error bounds arederived under certain assumptions that are interplay between the data regularity, regularity pick-up of L and the fractional power α . The needed justification is provided by the work of Bonito andPasciak, [3], for the problem (2) in the case when q = 0 and Γ D = ∂ Ω, see also earlier work [20].Following [3], one introduces ˙ H α := { v ∈ L : (cid:80) ∞ j =1 λ αj | ( v, ψ j ) | < ∞} and shows that ˙ H α is aHilbert space under the inner product A α ( v, w ) := ( L α/ v, L α/ w ), for all v, w ∈ ˙ H α . To set up afinite element approximation of L α u = f we first introduce its weak form: find u ∈ ˙ H α such that(4) A α ( u, v ) = ( f, v ) , ∀ v ∈ ˙ H α . This problem has a unique solution(5) u = T α f := ∞ (cid:88) j =1 λ − αj ( f, ψ j ) ψ j . Then for a finite elements space V h ⊂ ˙ H of continuous piece-wise polynomial functions definedon a quasi-uniform mesh with mesh size h one gets an approximate solution u h to (3) by setting(6) A αh u h = π h f, A − h = T h , where T h : V h → V h is the solution operator to the FEM of finding u h ∈ V h s. t. A ( u h , v ) = ( f, v ) ≡ ( π h f, v ) , ∀ v ∈ V h and π h : L (Ω) → V h is the orthogonal projection on V h . Here for T αh we use an expression similarto (5) but involving the eigenfunctions and eigenvalues of T h . As shown in [3], if the operator T satisfies the regularity pick up, i.e., there is s ∈ (0 ,
1] s.t. (cid:107) u (cid:107) ˙ H s (Ω) ≡ (cid:107)T f (cid:107) ˙ H s (Ω) ≤ c (cid:107) f (cid:107) ˙ H − s (Ω) and L is a bounded map of ˙ H s (Ω) into ˙ H − s (Ω), then for α > s one has(7) (cid:107) u − u h (cid:107) L = (cid:107)T α f − T αh π h f (cid:107) L ≤ Ch s (cid:107) f (cid:107) ˙ H δ , δ ≥ . The paper [3, see, Theorem 4.3] contains more refined results depending on the relationship betweensmoothness of the data δ , the regularity pick up s and the fractional order α . In the case of fullregularity, s = 1, the best possible rate for f ∈ L (Ω) is, cf. [3, Remark 4.1], (cid:107) u − u h (cid:107) L ≤ Ch α | ln h |(cid:107) f (cid:107) L . The bottomline of the error estimates from [3] is that, if f ∈ L (Ω) and s > α , then (cid:107) u − u h (cid:107) L is essentially O ( h α ). Therefore, the solution u h ∈ V h of (6) is an approximation to the solution u of (3). This fact makes our aim of solving the algebraic problem (6) justifiable. Finite elementapproximations of the elliptic problem (3) and also more general non-symmetric problems wererecently considered and studied by Bonito and Pasciak in [4].1.2. Algebraic problem under consideration.
Now let N be the dimension of V h and considerthat a standard nodal basis is used. Let A ∈ R N × N be a matrix representation of A h = T − h (defined in (3)) and (cid:101) u ∈ R N and (cid:101) f ∈ R N be vector representations through the nodal values of OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 3 u h ∈ V h , and π h f ∈ V h , respectively. Then we can recast the problem (6) in the following algebraicform:(8) find u ∈ R N such that A α u = (cid:101) f . The matrix A is symmetric and positive definite. The fractional power A α , 0 < α <
1, of a symmet-ric positive definite matrix A is expressed through the eigenvalues and eigenvectors { ( (cid:101) Λ i , Ψ i ) } Ni =1 of A . We assume that the eigenvectors are l -orthonormal, i.e. Ψ Ti Ψ j = δ ij and (cid:101) Λ ≤ (cid:101) Λ ≤ . . . (cid:101) Λ N .The spectral condition number k( A ) = (cid:101) Λ N / (cid:101) Λ = O ( h − ) for quasi-uniform meshes with mesh-size h .Then A = W DW T and A α = W D α W T , where W, D ∈ R N × N are defined as W = [ Ψ T , Ψ T , ..., Ψ TN ]and D = diag ( (cid:101) Λ , . . . , (cid:101) Λ N ). Then A − α = W D − α W T and the solution of A α u = (cid:101) f can be expressedas(9) u = A − α (cid:101) f = W D − α W T (cid:101) f . Obviously, we have the following standard equality for any β ∈ R : (cid:107) u (cid:107) A β + α = (cid:107) (cid:101) f (cid:107) A β − α with (cid:107) u (cid:107) A γ = u T A γ u , γ ∈ R . The formula (9) could be used in practical computations if the eigenvectors and eigenvalues areexplicitly known and the matrix vector multiplication with W is equivalent to a Fast FourierTransform when A is a circulant matrix. In such cases the computational complexity is almostlinear, O ( N log N ). However, this limits the applications to problems with constant coefficientsin simple domains and to the lowest order finite element approximations. More general is theapproach using an approximation of A with H -matrices combined with Kronecker tensor-productapproximation. This allows computations with almost linear complexity of the inverse of fractionalpower of a discrete elliptic operator in a hypercube (0 , d ∈ R d , for more details, see [9]. Firstattempt to apply H -matrices to solving fractional differential equations in one space variable isdone in [32].This work is related also to the more difficult problem of stable computations of the matrix squareroot and other functions of matrices, see, e.g. [7, 14, 16], where the stabilization of Newton methodis achieved by using suitable Pad´e iteration. However, in this paper we do not deal with evaluationof A α , instead we propose an efficient method for solving the algebraic system A α u = (cid:101) f , where A is an SPD matrix generated by approximation of second order elliptic operators. Our researchis also connected with the work done in [15], where numerical approximation of a fractional-in-space diffusion equation with non-homogeneous boundary conditions is considered. In [15], theproposed solver of the arising algebraic system relies on Lanczos method. First, the adaptivelypreconditioned thick restart Lanczos procedure is applied to a system with A . The gathered spectralinformation is then used to solve the system with A α . In [7] an extended Krylov subspace methodis proposed, originating by actions of the SPD matrix and its inverse. It is shown that for thesame approximation quality, the variant of the extended subspaces requires about the square rootof the dimension of the standard Krylov subspaces using only positive or negative matrix powers.A drawback of this method is the memory required to store the full dense matrix W needed toperform the reorthogonalization. Essentially, this approach and the method proposed and used in[12] rely on polynomial approximation of t − α , and the efficiency of the methods depends on thecondition number of A and deteriorates substantially for ill-conditioned matrices.1.3. Overview of existing methods.
The numerical solution of nonlocal problems is ratherexpensive. The following three approaches (A1 - A3) are based on transformation of the originalproblem to a local elliptic or pseudo-parabolic problem, or on integral representation of the solution,thus increasing the dimension of the original computational domain. The Poisson problem isconsidered in the related papers refereed bellow. A1 A Neumann to Dirichlet map is used in [5]. Then, the solution of fractional Laplacianproblem is obtained by u ( x ) = v ( x,
0) where v : Ω × R + → R is a solution of the equation − div (cid:0) y − α ∇ v ( x, y ) (cid:1) = 0 , ( x, y ) ∈ Ω × R + , S. HARIZANOV, R. LAZAROV, P. MARINOV, S. MARGENOV, Y. VUTOV where v ( · , y ) satisfies the boundary conditions of (2) ∀ y ∈ R + , lim y →∞ v ( x, y ) = 0 , x ∈ Ω , aswell as lim y → + (cid:0) − y − α v y ( x, y ) (cid:1) = f ( x ) , x ∈ Ω . The variational formulation of this equa-tion is well posed in the related weighted Sobolev space. The finite element approximationuses the rapid decay of the solution v ( x, y ) in the y direction, thus enabling truncation ofthe semi-infinite cylinder to a bounded domain of modest size. The proposed multilevelPCG method is based on the Xu-Zikatanov identity [30]. A2 A fractional Laplacian is considered in [27, 28] assuming the boundary condition a ( x ) ∂u∂n + µ ( x ) u = 0 , x ∈ ∂ Ω , which ensures L = L ∗ ≥ δ I , δ >
0. Then the solution of the nonlocal problem u canbe found as u ( x ) = w ( x, , w ( x,
0) = δ − α f, where w ( x, t ) , < t <
1, is the solution ofpseudo-parabolic equation ( t D + δ I ) dwdt + α D w = 0 , and D = L − δ I ≥
0. Stability conditions are obtained for the fully discrete schemes underconsideration. A further development of this approach is presented in [18] where the caseof fractional order boundary conditions is studied. A3 The following representation of the solution operator of (2) is used in [3]: L − α = 2 sin( πα ) π (cid:90) ∞ t α − (cid:0) I + t L (cid:1) − dt The authors introduce an exponentially convergent quadrature scheme. Then, (the approx-imation of u only involves evaluations of ( I + t i A ) − f , where t i ∈ (0 , ∞ ) is related to thecurrent quadrature node, and where I and A stand for the identity and the finite elementstiffness matrix corresponding to the Laplacian. A further development of this approach isavailable in [4], where the theoretical analysis is extended to the class of regularly accretiveoperators.There are various other problems leading to systems with fractional power of sparse symmetricand positive definite matrices. For an illustration we give the following examples. Consider anon-overlapping domain decomposition (DD) for the 2D model Laplacian problem on a regularmesh with an interface along a single mesh line. The elimination of all degrees of freedom from theinterior of the two subdomains reduces the solution to a system of equations on the interface, orthe Schur complement system. The matrix of this system is spectrally equivalent to A / , where A = h − tridiag ( − , , − A / , then the DD preconditioner has almost optimalcomplexity.In [12] the approach discussed in this paper was used to optimally solve the equation (8) with α = 0 .
5, when A is an SPD matrix and belongs to a class of particular weighted graph Laplacianmodels used in the volume constrained 2-phase segmentation of images. Then, after rescaling, sothat the spectrum of the rescaled matrix is in (Λ , t − / , t ∈ [Λ , well separated from 0, was used to construct an optimal solver.1.4. Our approach and contributions.
In Section 2 we introduce the mathematical problemand present the idea of the proposed algorithm. Let Λ be an upper bound for the spectrum of A ,namely, (cid:101) Λ j ≤ Λ, j = 1 , . . . , N . We rescale the system to the form(10) A α u = f , where A = A / Λ and f = (cid:101) f / Λ α , so that the spectrum of A , Λ j = (cid:101) Λ j / Λ, j = 1 , · · · , N is in (0 , A in the following assumption. Assumption 1.1. A is a symmetric, positive definite matrix and its spectrum is in the interval (0 , . OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 5
Next, we argue that instead of the system A α u = f one can solve the equivalent system A α − β u = A − β f with β ≥ A β − α f by P k ( A )( Q k ( A )) − f , for k integer, where P k ( t ) Q − k ( t ) := r βα ( t ) is the best uniform rational approxi-mation (BURA) of t β − α on the interval (0 , r βα ( t ), its approximation properties and questions regarding theimplementation of P k ( A )( Q k ( A )) − f .The properties of the best rational approximation have been an object of numerous studies. Inparticular, the distribution of the poles, zeros, and extreme points, and the asymptotic behavior ofthe error E α ( k, k ; β ) = max t ∈ [0 , (cid:12)(cid:12)(cid:12) t β − α − r βα ( t ) (cid:12)(cid:12)(cid:12) when k → ∞ are given in [26]. For example, it isknown that all poles lie on the negative real line and the error decays exponentially in k , namely,is O ( e − c √ k ), c >
0, see, relation (14).In Theorem 2.3 and Remark 2.4 we show that we can balance the finite element error (7) withthe error of the BURA (13) so that the total error is O ( h α ) when k ≈ β π ( β − α ) | ln h | . Thus,the feasibility of the method will depend on the possibility to address two key issues: (1) fora given 0 < α < k integer, compute the BURA P k ( t ) /Q k ( t ) and (2) implement P k ( A )( Q k ( A )) − f efficiently.To find the BURA for t β − α we apply the modified Remez algorithm, see, [19, 6]. The maindifficulty in implementing the algorithm is its instability for large k , outlined for example in [8].Our experience shows that for moderate k = 5 , , α < . k > P k ( A )( Q k ( A )) − will involve inversion of A − d j I , d j ≤ j = 0 , , . . . , k , see representation (15). The integer parameter k ≥ r α ( t ) of t − α on the interval (0 , β = 1 leads to (16). The nonpositivity of d j ensures thatthe systems with A − d j I , d j ≤ A . The positivity of c j means that the BURA approximation r α ( A ) is positive. Thebehavior of c j > k ≤
10 and small α we have developed, studied andexperimented with a new concept of multi-step BURA algorithm, outlined in Section 4.Finally, in Section 5 we present numerical experiments that illustrate the efficiency of the pro-posed algorithms. The first group of tests concerns scaled (normalized) matrices corresponding to1D Poisson equation where the exact solution is known and the BURA approximations are exactlycomputed. This setting allows numerically confirming the sharpness of theoretical estimates. Inparticular, some promising approximation properties are observed when different powers of A areinvolved in the multi-step BURA. The experiments with 2D fractional Laplacian illustrate the the-oretical results concerning balancing the rescaling effect in (13). Finally, we present 3D numericalexperiments involving jumping coefficients. The solution u is unknown, while u r is computed by apreconditioned conjugate gradient (PCG) solver that uses algebraic multigrid as a preconditioner.A multi-step setting of BURA is used to confirm the robustness with respect to the PCG accuracy.2. Solution strategy
The idea and theoretical justification of the method.
The goal of this study is topresent a new robust solver of optimal complexity for solving the system (8) for a large class ofsparse SPD matrices assuming that such a solver is available for α = 1. This assumption holds ina very general setting, when A is generated by finite element or finite difference approximation ofsecond order elliptic operators. Such matrices are used in the numerical tests presented in Section5. Note that A α is dense and in general not known. This means in particular that the standard S. HARIZANOV, R. LAZAROV, P. MARINOV, S. MARGENOV, Y. VUTOV iterative solution methods are not applicable since even in the case when A α v is computable justone such computation requires O ( N ) arithmetic operations.We consider the class of rational functions R ( m, k ) := { r ( t ) = P m ( t ) /Q k ( t ) : P m ∈ P m , Q k ∈ P k } , where P k is the set of all polynomials of degree k . For a given univariate function g ( t ), 0 ≤ t ≤ r ∗ ( t ) = P ∗ m ( t ) Q ∗ k ( t ) ∈ R ( m, k ) of the problem(11) min r ∈R ( m,k ) max t ∈ [0 , | g ( t ) − r ( t ) | = max t ∈ [0 , | g ( t ) − r ∗ ( t ) | is called Best Uniform Rational Approximation (BURA) of g ( t ). Deffinition 2.1.
The minimizer r βα ( t ) = P ∗ m ( t ) Q ∗ k ( t ) for g ( t ) = t β − α is called β -Best Uniform RationalApproximation ( β -BURA) and its error is denoted by E α ( m, k ; β ) := max t ∈ [0 , (cid:12)(cid:12)(cid:12) t β − α − r βα ( t ) (cid:12)(cid:12)(cid:12) . Our algorithm is based on the following Lemma:
Lemma 2.1.
Let A satisfy the assumption 1.1 and F = A − β f so that u = A β − α F . Let r βα ( t ) be the best uniform rational approximation of t β − α on [0 , and consider u r = r βα ( A ) F to be anapproximation to u . Then the following bound for the error holds true (cid:107) u r − u (cid:107) A γ ≤ E α ( m, k ; β ) (cid:107) f (cid:107) A γ − β ∀ γ ∈ R . (12) Proof.
Consider the representation of F with respect to the eigenvectors of A , F = (cid:80) Ni =1 F i Ψ i sothat (cid:107) F (cid:107) A γ = (cid:80) Ni =1 Λ γi F i , for any γ ∈ R . Since r βα is analytic in (0 , r βα ( A ) Ψ i = r βα (Λ i ) Ψ i , i = 1 , . . . , N . Using the orthonormal propertyof the eigenvectors Ψ Ti Ψ j := (cid:104) Ψ i , Ψ j (cid:105) = δ ij we get easily (cid:107) u r − u (cid:107) A γ = (cid:107) r βα ( A ) F − A β − α F (cid:107) A γ = (cid:42) N (cid:88) i =1 ( A γ ( r βα ( A ) − A β − α ) F i Ψ i , N (cid:88) i =1 ( r βα ( A ) − A β − α ) F i Ψ i (cid:43) = N (cid:88) i =1 F i Λ γi (cid:16) r βα (Λ i ) − Λ β − αi (cid:17) ≤ max t ∈ [0 , | r βα ( t ) − t β − α | N (cid:88) i =1 Λ γi F i ≤ E α ( m, k ; β ) (cid:107) F (cid:107) A γ . To complete the proof take into account that F = A − β f . (cid:3) It is important to keep in mind that the above estimate is for the scaled system (10), where f = (cid:101) f / Λ α and A = A / Λ. As a corollary we get the following bound for the solution through theoriginal (unscaled) data:
Corollary 2.2.
The following estimate holds true for the solution of (8) : (13) (cid:107) u r − u (cid:107) A γ ≤ E α ( m, k ; β )Λ β − α (cid:107) (cid:101) f (cid:107) A γ − β . Among various classes of best rational approximations, the diagonal sequences r ∈ R ( k, k ) ofthe Walsh table of t α , 0 < α < E α ( k, k ; β ) when k → ∞ are well known. For example, Theorem 1 from[26] shows that lim k →∞ e π √ ( β − α ) k E α ( k, k ; β ) = 4 β − α | sin π ( β − α ) | OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 7 holds for any 0 < α < ≤ β , β integer. This could be written as an approximate relation (used inpractice)(14) E α ( k, k ; β ) ≈ β − α | sin π ( β − α ) | e − π √ ( β − α ) k . In order to convince ourselves in the feasibility of practical use of (14) in Table 1 we present theobtained results for E α ( k, k ; β ) for various k and β when the BURA P ∗ k ( t ) /Q ∗ k ( t ) is computed usingthe Remez algorithm, see for more details Section 3. The results show that for relatively small k we can get good approximation E α ( k, k ; β ). It is quite clear from this table that E α ( k, k ; β ) is acouple of orders of magnitude smaller for β = 3 compared with β = 1. However, this comes at acost. From (13) we observe that: (1) the errors are measured in two different ways and (2) thescaling factor enters into the play with a negative impact on the accuracy for larger β . Thus, thevalues β = 2 and β = 3 have not been used in our computations, we are giving the approximationproperties of BURA for these values just for comparison. In Table 1 in parenthesis we show thecomputed values from the asymptotic formula (14). These and other computations, see e.g. [29,Tables 2.1 – 2.7], show that asymptotic formula is quite accurate and the relation (14) could beused for fairly low k . α E α (5 ,
5; 1) E α (6 ,
6; 1) E α (7 ,
7; 1) E α (5 ,
5; 2) E α (5 ,
5; 3)0.75 2.7348E-3 (3.60E-3) 1.4312E-3 7.8650E-4 (9.82E-4) 1.9015E-6 6.8813E-80.50 2.6896E-4 (3.88E-4) 1.0747E-4 4.6037E-5 (6.28E-5) 9.5789E-7 5.5837E-80.25 2.8676E-5 (4.16E-5) 9.2522E-6 3.2566E-6 (4.47E-6) 2.8067E-7 2.4665E-8
Table 1.
Errors E α ( k, k ; β ) of BURA P ∗ k ( t ) /Q ∗ k ( t ) of t β − α on [0 , Lemma 2.3.
For β ≥ integer there is a constant C α,β > and an integer k ≥ such that for k ≥ k the following error bound holds true (cid:107) u r − u (cid:107) A γ ≤ C α,β Λ β − α e − π √ ( β − α ) k (cid:107) (cid:101) f (cid:107) A γ − β . Remark 2.4.
The error in solving the algebraic problem (8) using the proposed method shouldbe balanced with the approximation error given by (7) . In the case of second order problems on aquasi-uniform mesh with size h we have Λ ≈ h − . Then in case of best possible convergence rate ofthe finite element solution, namely, O ( h α | ln h | ) , cf. [3, Remark 4.1] , we can take k ≈ β π ( β − α ) | ln h | , and get the total error O ( h α | ln h | ) (this includes the finite element approximation error and errorof approximately solving the algebraic problem). This represents the foundation of the method we propose and study in this paper. The feasibilityof such approach depends substantially on the possibility to efficiently compute r βα ( A ) f . Onepossible implementation is proposed in next subsection.2.2. Efficient implementation of the method.
We first bring some important facts about thebest uniform rational approximation r α ( t ), m = k , of t − α on [0 ,
1] for 0 < α <
1, see, e.g. [24, 26].
Lemma 2.5. ( [24, Lemma 2.1] ) Let m = k and < α < . Then the following statements arevalid: (1) The best rational approximation r α ( t ) has numerator and denominator of exact degree k ; (2) All k zeros ζ , . . . , ζ k and poles d , . . . , d k of r α are real and negative and are interlacing,i.e. with appropriate numbering one has > ζ > d > ζ > d > · · · > ζ k > d k > −∞ ; S. HARIZANOV, R. LAZAROV, P. MARINOV, S. MARGENOV, Y. VUTOV (3)
The function t − α − r α ( t ) has exactly k + 2 extreme points η , . . . , η k +2 on [0 , and withappropriate numbering we have η < η < · · · < η k +2 = 1 ,η − αj − r α ( η j ) = ( − j E α ( k, k ; 1) , j = 1 , . . . , k + 2 . Then we introduce d = 0 so that r α ( t ) is represented as a sum of partial fractions(15) t − r α ( t ) := 1 t P ∗ k ( t ) Q ∗ k ( t ) = 1 t (cid:80) kj =0 p j t j (cid:80) kj =0 q j t j = k (cid:88) j =0 c j t − d j . These notations are used in the tables below.This Lemma allows us to have the following implementation of the method:Step 1: Find all poles 0 = d > d > d > · · · > d k ;Step 2: Find the representation (15) of r k ( t ) as a sum of partial fractions;Step 3: Compute the approximate solution by(16) u r := A − r α ( A ) f = k (cid:88) j =0 c j ( A − d j I ) − f This shows that to find u r = r α ( A ) A − f we need to solve one system A v = f and k separateindependent systems ( A − d i I ) v = f for i = 1 , . . . k with SPD matrices A − d i I . Remark 2.6.
Our numerical tests show that often we achieve accuracy of E α ( k, k ; 1) ≈ − orbetter with k = 5 . For example, for α = 0 . and k = 5 we get E α ( k, k ; 1) = 2 . ∗ − . Thecoefficients of the 1-BURA are given in Table 3 and they result in the following partial fractionsrepresentation: r . ( t ) t = P ∗ ( t ) t Q ∗ ( t ) = 0 . t + 0 . t + 0 . . t + 0 . . t + 0 . . t + 0 . . t + 3 . . We note that the nonpositivity of d j ensure that the sytems ( A − d j I ) v = f can be solvedefficiently and the positivity of c j , shown in [11, Theorem 1], guarantees no loss of significant digitsdue to subtraction of large numbers. Remark 2.7.
Using representation of r α ( t ) by partial fractions is just one possible way to compute r α ( A ) A − f . Another possibility is to use the zeros and the poles to compute consecutively the factorsin the formula r α ( A ) A − f = c k (cid:89) j =1 ( A − ζ j I )( A − d j I ) − A − f . Due to the interlacing of the zeroes and the poles this will lead to stable computations. Moreover,it will preserve the monotonicity of the solution (see [11] ), which is a desired feature in someapplications. The only substantial difference is that the computations with partial fractions can bedone in parallel.
Now we present some examples of BURA r α within the class R ( k, k ) for α = 0 . , . , .
25. InTables 2 - 4 we show the computed coefficients (15) of r α for α = 0 . , . , . A − α via the application of theoperator A − r α ( A ) involves solving six systems of linear equations with SPD matrices; wehave assumed that each evaluation of ( A− d j I ) − f can be computed approximately by PCGmethod with optimal complexity.(2) Summing these six solutions is a stable process since the coefficients in the sum of fractions(17) are small and positive and there should not expect any loss of accuracy (or stability)that might come from subtracting large numbers. OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 9 j p j q j c j d j Table 2.
The coefficients in the representation (15) of the best rational ap-proximation P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , α = 0 .
75; from Table 1 we have E α (5 ,
5; 1) =2.7348E-3 j p j q j c j d j Table 3.
The coefficients in the representation (15) of the best rational ap-proximation P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , α = 0 .
5; from Table 1 we have E α (5 ,
5; 1) =2.6896E-4 j p j q j c j d j Table 4.
The coefficients in the representation (15) of the best rational ap-proximation P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , α = 0 .
25; from Table 1 we have E α (5 ,
5; 1) =2.8676E-53.
Best uniform rational approximation of t β − α Theoretical background and numerical methods.
Let r βα ( t ) = P ∗ m ( t ) Q ∗ k ( t ) be the β -BURA of t β − α for t ∈ [0 , k and m , the rational function has the following representation:(18) P ∗ m ( t ) Q ∗ k ( t ) = m (cid:80) j =0 p j t jk (cid:80) j =0 q j t j = m (cid:80) j =0 ¯ p j T j (2 t − k (cid:80) j =0 ¯ q j T j (2 t −
1) = ¯ P ∗ m ( s )¯ Q ∗ k ( s )where T ( s ) = 1 , T ( s ) = s, . . . , T j ( s ) = 2 sT j − ( s ) − T j − ( s ) , j = 2 , , . . . ; s ∈ [ − ,
1] areorthogonal base functions. These are the well-known Chebyshev polynomials and in our case s = 2 t −
1, since t ∈ [0 , ,
1] the element of best uniformapproximation exists. Due to the equioscillation theorem, there are at least ( m + k + 2) points { η i } m + k +21 , where the error r βα ( t ) − t β − α have extremes and the sign alternates. We use orthogonalbase functions, because there are numerical difficulties (instabilities) for finding r βα in the standardmonomial basis { t j } (see [8]). The benefits of working in Chebyshev basis are illustrated in Ta-ble 5, where the maximal values of m and k for which the element of best uniform approximationof t − α can be successfully computed (the algorithm converges) are documented. The left pairsin the table correspond to best polynomial approximation ( k = 0), while the right ones corre-spond to best ( k, k )-rational approximation ( m = k ). It is evident that apart from the choice ofbase functions, calculations heavily depend on the used precision for arithmetic operations (single,double, quadruple). This is due to the non-differentiability of t − α at zero. The function is only(1 − α )-H¨older continuous (i.e. in C , − α [0 , { η i } k +21 of r α are clustered in a neighborhood of zero to account for the steep slope there. For example,when k = 5 and α = 0 .
75 the first two points are η = 0 and η ≈ · − , while the ninth pointvalue is still just η ≈ .
05. Therefore, to accurately compute { η i } and capture the sign changesof r α ( t ) − t − α between them one must use high precision arithmetics. This fact has been knownand attempts to compute the BURA and the error E α ( k, k ; β ) for large k has required using highprecision arithmetic, e.g., see [29].Precision Base α = 0 . α = 0 . α = 0 . t j (8,0), (4,4) (8,0), (4,4) (8,0), (2,2)Single T j ( s ) (40,0), (3,3) (50,0), (2,2) (50,0), (2,2)Double t j (19,0), (6,6) (19,0), (4,4) (19,0), (2,2)Double T j ( s ) (60,0), (5,5) (60,0), (5,5) (50,0), (4,4)Quadro t j (35,0), (6,6) (35,0), (4,4) (35,0), (2,2)Quadro T j ( s ) (95,0), (11,11) (95,0), (11,11) (95,0), (7,7) Table 5.
Maximal values ( m, k ) for which Algorithm 3.1 converges.3.2.
Modified Remez algorithm for computing BURA.
We suggest the following (modifiedRemez) algorithm for finding the β -BURA of t β − α on [0 , f ( s ) = (cid:0) s (cid:1) β − α for s ∈ [ − ,
1] (see [19], [6]).
Algorithm 3.1. Input: ( α, β ) , ( m, k ) , N (maximal number of algorithm iterations), V (maximalnumber of inside iterations for solving the non-linear system in Step 3(ii)), δ > (accuracy). Initialization: (cid:96) , s (0) , and ¯ r , satisfying • (cid:96) = m + k + 2 . • (cid:110) s (0) i (cid:111) (cid:96)i =1 – strictly monotonically increasing sequence in [ − , . • ¯ r ( s ) = ¯ P m ( s )¯ Q k ( s ) : (cid:16) f ( s (0) i ) − ¯ r ( s (0) i ) (cid:17) / (cid:16) f ( s (0) i +1 ) − ¯ r ( s (0) i +1 ) (cid:17) < , ∀ i = 1 , . . . , (cid:96) − . FOR n = 1 , , . . . DO (1) Updating the equioscillation point set: FOR i = 1 , . . . , (cid:96) DO (i) ¯ τ ( n ) i := sup − ≤ τ ≤ s ( n − i { f ( τ ) = ¯ r n − ( τ ) } , ¯ τ ( n ) i := inf s ( n − i ≤ τ ≤ { f ( τ ) = ¯ r n − ( τ ) } . (ii) s ( n ) i = argmax ¯ τ ( n ) i ≤ s ≤ ¯ τ ( n ) i | f ( s ) − ¯ r n − ( s ) | , η ( n ) i = | f ( s ( n ) i ) − ¯ r n − ( s ( n ) i ) | . END FOR (iii) s ( n ) ∗ = argmax − ≤ s ≤ | f ( s ) − ¯ r n − ( s ) | , η ( n ) ∗ = | f ( s ( n ) ∗ ) − ¯ r n − ( s ( n ) ∗ ) | . OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 11 (iv) IF (cid:16) s ( n ) ∗ / ∈ { s ( n ) i } (cid:96) (cid:17) THEN FIND j s.t. s ( n ) j < s ( n ) ∗ < s ( n ) j +1 .IF (cid:16) sgn (cid:16) f ( s ( n ) ∗ ) − ¯ r n − ( s ( n ) ∗ ) (cid:17) = sgn (cid:16) f ( s ( n ) j ) − ¯ r n − ( s ( n ) j ) (cid:17)(cid:17) THEN s ( n ) j = s ( n ) ∗ .ELSE s ( n ) j +1 = s ( n ) ∗ . (2) Convergence check: IF (cid:16) max i η ( n ) i − min i η ( n ) i < δ (cid:17) OR ( n = N + 1) STOP . ELSE (3) Updating the rational approximation:
Solve iteratively the non-linear system (cid:16) f (cid:16) s ( n ) i (cid:17) − ¯ r n (cid:16) s ( n ) i (cid:17)(cid:17) = ( − i E n , i = 1 , . . . , (cid:96) for the unknown E n and the coefficients of ¯ r n : (i) ( E n , ¯ r n ) = ( E n − , ¯ r n − ) . (ii) FOR v = 1 , , . . . DO: Solve the (cid:96) × (cid:96) linear system of equations m (cid:88) j =0 ¯ p ( n,v ) j T j ( s ( n ) i ) − (cid:16) f ( s ( n ) i ) − ( − i E ( v − n (cid:17) k (cid:88) j =1 ¯ q ( n,v ) j T j ( s ( n ) i ) + ( − i E ( v ) n = f ( s ( n ) i ) . IF (cid:12)(cid:12)(cid:12) E ( v ) n − E ( v − n (cid:12)(cid:12)(cid:12) < (cid:15) , OR v > V GO to (iii). ELSE v = v + 1 and REPEAT. (iii) ¯ r n ( s ) = (cid:16)(cid:80) mj =0 ¯ p ( n,v ) j T j ( s ) (cid:17) / (cid:16) (cid:80) kj =1 ¯ q ( n,v ) j T j ( s ) (cid:17) . (iv) E n = E ( v ) n . (4) n = n + 1 . GO to Step 1. Output: n ; ¯ r βα ( s ) = ¯ r n ; E α ( m, k ; β ) = | E n | ; s ∗ = { s ( n ) i } (cid:96)i =1 . Then we take η = ( s ∗ + 1) / r βα ( t ) from ¯ r βα ( s ) using (18).Several remarks, concerning the computer implementation of the proposed algorithm are in order:(1) In the Initialization step, we usually take s (0) to be uniformly sampled on [ − , r we apply least-squares optimization techniques. All the rationalfunctions ¯ r n are normalized with respect to the constant term in the denominator, i.e.,¯ r n ( s ) = m (cid:88) j =0 ¯ p ( n ) j T j ( s ) / k (cid:88) j =1 ¯ q ( n ) j T j ( s ) ∀ n ≥ . (2) In order to increase the computational efficiency of the algorithm, we compute neither thesequences ¯ τ ( n ) and ¯ τ ( n ) in Step 1(i) nor the local extrema s ( n ) in Step 1(ii) . Instead, we searchfor the maximal value of | f ( s ) − ¯ r n − ( s ) | on a small, discretized interval around s ( n − i , decreasethe mesh size, and repeat the process several times around the current maximizer. Such simplelocalization techniques seem to work fine for our numerical examples.(3) In Step 3(ii) we apply Aitken-Steffensen acceleration but instead of E ( v − n for the systemsplitting we use a combination of the values { E ( v − i ) n } i =1 from the previous three steps.4. Numerical accuracy and Multi-step BURA method
In this section we investigate the numerical accuracy of the proposed algorithm and a multi-step generalization of the BURA-method. The analysis, presented here is theoretical in nature,so we consider the full generality of the proposed solution strategy, namely the ( m, k ) β -BURAapproximation.4.1. Properties of the fractional decomposition.
For given ( m, k, β ) the partial fraction de-composition of t − β r βα ( t ) has the general form(19) t − β r βα ( t ) = m − k − β (cid:88) j =0 b j t j + β (cid:88) j =1 c ,j t j + p (cid:88) j =1 c j t − d j + p (cid:88) j =1 B j t + C j ( t − F j ) + D j where k = p + 2 p . We always consider triples for which m < k + β . One reason for such aparameter constraint comes from the fact that t − β r βα has a leading term of degree t m − k − β , whileit approximates the power function t − α , α >
0. Another reason is the numerical simplificationof (19), where the index set for the first sum becomes empty. In all our numerical examples thedenominator of r βα has no complex roots, thus we concentrate on the case p = 0 from now on.Then β -BURA can be rewritten in the following way:(20) 1 t β P ∗ m ( t ) Q ∗ k ( t ) = m (cid:80) j =0 p j t j t β (cid:32) k (cid:80) j =0 q j t j (cid:33) = β (cid:88) j =1 c ,j t j + k (cid:88) j =1 c j t − d j . The first representation in (20) is the best approximation written as a standard rational function,while the second one is its partial fraction decomposition (19), the way this approximation is usedin the implementation of the method.Let P ∗ m ( t ) Q ∗ k ( t ) = β − (cid:88) j =0 b ∗ j t j + k (cid:88) j =1 c ∗ j t − d j . We have used that β > m − k and we set the extra coefficients { b ∗ j } β − m − k +1 to zero whenever m − k < β −
1. Then, straightforward computations give rise to(21) 1 t β P ∗ m ( t ) Q ∗ k ( t ) = β (cid:88) j =1 b ∗ β − j t j + k (cid:88) j =1 (cid:32) c ∗ j /d βj t − d j − β (cid:88) i =1 c ∗ j /d β − i +1 j t i (cid:33) Comparing the coefficients in front of the corresponding terms in (20) and (21), we derive(22) c ,j = b ∗ β − j − k (cid:88) i =1 c ∗ i /d β − j +1 i , c j = c ∗ j /d βj . Various useful identities follow from (22). We want to highlight a couple of them. Due to theChebyshev’s equioscillation theorem(23) c ,β = b ∗ − k (cid:88) i =1 c ∗ j /d j = P ∗ m (0) Q ∗ k (0) = p q = ± E α ( m, k ; β ) . In particular, for ( k, k ; 1) we have c = E α ( k, k ; 1) due to Lemma 2.5.The second one is(24) c , + k (cid:88) i =1 c i = b ∗ β − = (cid:26) p m /q k , m − k = β − , m − k < β − . Finally, (22) allows for stable numerical computations of the coefficients { c j } , as the fractionaldecomposition of r βα can be accurately derived in Chebyshev basis.4.2. Accuracy Analysis.
In this subsection we briefly discuss issues related to the numericalaccuracy of the developed framework. We do not go into details, since thorough analysis of thealgorithm is outside the scope of the paper. However, certain observations in this direction areworth mentioning, so that the reader can make conclusions for the full picture.Lemma 2.1 quantifies the error between u r = r βα ( A ) A − β f and u = A − α f . All the estimationsare under the assumption that u r can be exactly computed by an optimal numerical solver. Withinthe adopted setup m < k + β and p = 0, u r has the following representation(25) u r = β (cid:88) i =1 c ,i A − i f + k (cid:88) i =1 c i ( A − d i I ) − f := β (cid:88) i =1 c ,i v ,i + k (cid:88) i =1 c i v i , OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 13 which is the corresponding simplification of (16). The practical derivation of u r involves k + β applications of such a solver, that independently solves each of the involved large-scale linearsystems with a right-hand-side f . The numerical stability of each solution process depends on thecondition number of the underlined linear operator. Remark 4.1.
The matrix
A − d i I is better conditioned than A whenever d i < or d i > Λ + Λ N .If d i > Λ + Λ N , the condition numler k( −A + d i I ) is uniformly bounded dependinding only on d i . Since we are interested in operators A which spectrum is normalized to lie inside (0 , ≈ N ≈
1. The poles of r βα are outside of (a neighborhood of) the unit interval [0 , d i ’s naturally satisfy thecondition in Remark 4.1. Thus, the numerical computation v MGi of v i = ( A − d i I ) − f , i = 1 , . . . , k is a stable process.When β = 1, using the notation (15), we observe that the most time-consuming procedure is thederivation of v MG that corresponds to inverting A ( v = A − f ). In our numerical experiments, weuse algebraic multigrid (AMG) [13] as a preconditioner in a CG method. The same preconditionercan be used to operators like A− d i I . We already observed that, provided m = k , all the coefficients c i are positive and sum to p m /q m (see (24)). The ratio p m /q m increases with α (see Tables 2 - 4 for m = k = 5) but seems to always be O (1). Therefore, for this setting the accuracy of the numericalderivation of u MGr is proportional to the accuracy of inverting A . Hence, numerics are trustworthyin general and unsubstantial additional errors for u MGr − u are accumulated.4.3. Further analysis in the case β > . When β >
1, since k( A β ) = k( A ) β , we need to findan approximate solution of a system with much worst condition number than the original system.This could cause loss of stability (or loss of accuracy). Since we solve A β v = f iteratively via β consecutive applications of the AMG solver as a preconditioner for A , we need to analyze thestability of such computational strategy. Lemma 4.2.
Let µ, ν, ε > be given and v = A − n f , where n ≥ . Assume that z MG is anumerical solution for z = A − ( n − f , while v MG is a numerical solution for ¯ v = A − z MG . Then (26) (cid:107) z MG − z (cid:107) A − (cid:107) z (cid:107) A − ≤ µε, (cid:107) v MG − ¯ v (cid:107) A (cid:107) ¯ v (cid:107) A ≤ νε imply (cid:107) v MG − v (cid:107) A (cid:107) v (cid:107) A ≤ ( µ + ν + µνε ) ε. Proof.
Applying triangle inequality, we derive (cid:107) v MG − v (cid:107) A ≤ (cid:107) v MG − ¯ v (cid:107) A + (cid:107) ¯ v − v (cid:107) A ≤ νε (cid:107) ¯ v (cid:107) A + (cid:107) ¯ v − v (cid:107) A ≤ νε (cid:107) v (cid:107) A + (1 + νε ) (cid:107) ¯ v − v (cid:107) A ; (cid:107) ¯ v − v (cid:107) A = (cid:107) A − ( z MG − z ) (cid:107) A = (cid:107) z MG − z (cid:107) A − ≤ µε (cid:107) z (cid:107) A − = µε (cid:107)A v (cid:107) A − = µε (cid:107) v (cid:107) A . These imply (cid:107) v MG − v (cid:107) A ≤ ( µ + ν + µνε ) ε (cid:107) v (cid:107) A , which completes the proof. (cid:3) Note that z serves as a right-hand-side for the linear system A v = z . Thus, for the stability incomputing v the z -related quantities need to be measured in the (cid:107) · (cid:107) A − norm.Now we are ready to quantify the accuracy of the numerical derivation of v ,β = A − β f under theproposed above computational strategy. Iteratively, we define v MG , to be the output of the appliedoptimal solver (the notation MG doesn’t mean that only multigrid solver can be used) to the linearsystem v , = A − f and v MG ,j +1 – numerical solution of A v = v MG ,j , j = 1 , . . . , β − . Corollary 4.3.
Let β ≥ and v MG ,β be derived as described above. Assume that a numerical solverfor the system A v = z , with arbitrary v , z ∈ R N has computed v MG with guaranteed relative error ε , i.e. (cid:107) v MG − v (cid:107) A / (cid:107) v (cid:107) A ≤ ε. Then (cid:107) v MG ,β − v ,β (cid:107) A / (cid:107) v ,β (cid:107) A ≤ a β ε with a β +1 = 1 + (1 + ε )k( A ) a β , a = 1 . Consequently, a β = O (k( A ) β − ) . Proof.
The proof is by induction. For β = 1, we have that v MG , is the solver output for v , = A − f .Thus, a = 1 follows directly from the assumption on the solver accuracy. Now, let (cid:107) v MG ,β − v ,β (cid:107) A / (cid:107) v ,β (cid:107) A ≤ a β ε holds true for β . Denote by ¯ v ,β +1 := A − v MG ,β . Again, due to the solver accuracy, we have (cid:107) v MG ,β +1 − ¯ v ,β +1 (cid:107) A / (cid:107) ¯ v ,β +1 (cid:107) A ≤ ε. Applying the obvious inequalities Λ (cid:104)A − v , v (cid:105) ≤ (cid:104)A v , v (cid:105) ≤ Λ N (cid:104)A − v , v (cid:105) we obtain (cid:107) v MG ,β − v ,β (cid:107) A − (cid:107) v ,β (cid:107) A − ≤ Λ − (cid:107) v MG ,β − v ,β (cid:107) A Λ − N (cid:107) v ,β (cid:107) A ≤ k( A ) (cid:107) v MG ,β − v ,β (cid:107) A (cid:107) v ,β (cid:107) A ≤ k( A ) a β ε. The result follows from Lemma 4.2 which we apply with µ = k( A ) a β and ν = 1. (cid:3) The coefficient c ,β = ± E α ( m, k ; β ), due to (23). Therefore the numerical accuracy for thecomputation of the term c ,β v ,β in u r (see (25)) depends on the product E α ( m, k ; β )k( A ) β − . For β > β grows: the first decreases(see, Table 1) while the second increases. As a result, we conclude that computing with β = 1 is areasonable practical choice.4.4. Multi-step BURA approximation.
From Table 1 we observe that when α increases, sodoes the error E α ( k, k ; β ). In particular, for β = 1, the quantities E . ( k, k ; 1), E . ( k, k ; 1) and E . ( k, k ; 1) are all of different order. This is due to the steeper slopes in a neighborhood of zero forthe function t − α , which results in higher and more frequent oscillations of the residual r α ( t ) − t − α there. Apart from such theoretical drawbacks, there are also additional numerical difficulties withthe convergence of Algorithm 3.1 as the set of extreme points { η i } k +21 for the residual clusteraround zero (see [24, Theorem 4]). Indeed, higher numerical precision is needed for the correctseparation of the extreme points, as well as more internal and external iterations are executed forsolving the ill-conditioned linear systems in Step 3(ii) and for reaching the stopping criterion of theAlgorithm, respectively. As an alternative approach, we study the possibility to replace the actionof r α by the joint action of several r α i rational functions, where each α i is smaller than the original α , thus r α i is cheaper to be generated and its approximation error E α i ( k, k ; 1) is smaller.Our idea is to apply a multi-step procedure, based on the identity A − α f = A − α n ◦ A − α n − ◦ · · · ◦ A − α f , n (cid:88) i =1 α i = α. First, we approximate A − α f by u := r α ( A ) A − f . Then we approximate A − α ◦ A − α f by u := r α ( A ) A − u and so on. Finally, we approximate u = A − α f by u n = r α n ( A ) A − u n − .Following (13) and setting γ = 1 we are interested in the theoretical and numerical behavior of theerror ratio (cid:107) u n − u (cid:107) A / (cid:107) f (cid:107) A − .The theoretical error analysis is based on Lemma 2.1. Denote by ε i ( t ) the residual of r α i ( t ) withrespect to t − α i . Then, for each i = 1 , . . . , n we have(27) r α i ( t ) = t − α i + ε i ( t ) , | ε i ( t ) | ≤ E α i ( k, k ; 1) ∀ t ∈ [0 , . The multi-step approximation u n can be rewritten in the form u n = n (cid:89) i =1 r α i ( A ) A − n f = (cid:32) A n − n (cid:89) i =1 r α i ( A ) (cid:33) A − f . Therefore the proof of Lemma 2.1 implies that we need to estimate the approximation error(28) E α ...α n ( k, k ; 1) := max t ∈{ Λ i } N (cid:12)(cid:12)(cid:12)(cid:12) r α ( t ) r α ( t ) . . . r α n ( t ) t n − − t − α (cid:12)(cid:12)(cid:12)(cid:12) . OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 15
Consider n = 2. Using (27) we obtain r α ( t ) r α ( t ) t = t − α + t − α ε ( t ) + t − α ε ( t ) + t − ε ( t ) ε ( t ) . Denote by E α i = E α i ( k, k ; 1), i = 1 ,
2. Since the spectrum of A is normalized and Λ N ≤
1, so that t − α i ≤ Λ − α i ≈ k ( A ) α i and from (28) we conclude(29) E α α ( k, k ; 1) ≤ E α k( A ) α + E α k( A ) α + E α E α k( A ) . Comparing with the error estimate E α ( k, k ; 1) we observe that unlike the direct approach, the two-step approximation error depends on the condition number of A , thus is not dimension-invariant ingeneral. Therefore, the benefits from the proposed multi-step procedure for computing A − α f arelimited in theory since the overall error is magnified by the condition number of A .On the other hand, the multi-step BURA method possesses several interesting properties thatare worth investigating further. First of all, it provides good approximation on the high part of thespectrum of A . From Lemma 2.5 we know that (cid:0) r α i ( t ) − t − α i (cid:1) (cid:12)(cid:12) t =1 = − E α i , i = 1 , , so for f = Ψ N we get (cid:107) u − A − α Ψ N (cid:107) A / (cid:107) Ψ N (cid:107) A − ≈ E α + E α − E α E α , since Λ N ≈
1. This is much better than E α ( k, k ; 1). Thus, as N → ∞ , the approach is beneficialwhen f ∈ span { Ψ N − (cid:96) , . . . , Ψ N } for some (cid:96) (cid:28) N . Second of all, especially if α (cid:54) = α , E α α mightremain significantly smaller than the right-hand-side of (29). Indeed, the extreme points of r α and r α are with high probability disjoint sets, therefore | ε ( t ) | and | ε ( t ) | cannot simultaneouslyattend their maximums, meaning that the factor E α E α in front of k( A ) is an overestimate.5. Numerical tests
A comparative analysis of the numerical accuracy of the proposed solvers and the related theo-retical estimates are presented in this section. The first group of tests concerns normalized matricesobtained from a three-point approximation of the Poisson equation in one space dimension. For thissetting we are able to directly compute the exact solution u = A − α f , as well as the approximatesolution u r = r βα ( A ) A − β f , thus no additional numerical errors are accumulated in the process. Thefirst experimental set is devoted to the numerical validation of Lemma 2.1. The second one studiespossible improvements in the accuracy of the approximation u r for larger α when a multi-stepapproximation process that involves smaller α ’s is applied. A third experiment deals with a 2Dfractional Laplacian operator and illustrates the rescaling effect in (13). Finally, we confirm the ac-curacy analysis in Section 4.2 by running 3D numerical experiments, where u is unknown, while u r is computed by a numerical solver that uses algebraic multigrid as a preconditioner in the conjugategradient method. In all the experiments we take m + 1 = k + β , thus we solve m + 1 linear systems α E α (5 ,
5; 1) E α (5 ,
4; 2) E α (5 ,
3; 3) E α (7 ,
7; 1) E α (7 ,
6; 2) E α (7 ,
5; 3)0.75 2.7348E-3 3.8415E-6 4.6657E-7 7.8650E-4 2.0108E-7 6.6194E-90.50 2.6896E-4 2.0349E-6 4.0421E-7 4.6037E-5 7.8577E-8 4.3899E-90.25 2.8676E-5 6.2333E-7 1.8958E-7 3.2566E-6 1.8043E-8 1.5792E-90.10 4.9432E-6 1.7490E-7 6.7114E-8 4.5139E-7 4.2824E-9 4.7675E-10
Table 6.
Errors E α ( m, m + 1 − β ; β ) of BURA P ∗ m ( t ) /Q ∗ m +1 − β ( t ) of t β − α on [0 , m = 5 , u r (see (16)). We consider m = { , } , β = { , } , and α = { . , . , . } .For each of the corresponding BURA functions all the zeros d i of the denominator are real and ofmultiplicity one, so only systems of the type ( A − d i I ) − f appear. The approximation errors are summarized in Table 6. In the discussion below the Euclidean norm of a vector in R N is denotedas (cid:96) -norm.5.1. Numerical validation of Lemma 2.1.
We consider the N × N stiffness matrix A , corre-sponding to a three-point finite difference approximation (or FE approximation with linear ele-ments) of the operator L u = − u (cid:48)(cid:48) with zero Dirichet boundary conditions on a uniform partitioningof (0 ,
1) with mesh-size h = 1 / ( N + 1). The tridiagonal matrix is normalized so that its spec-trum lies inside [0 ,
1] and has enrties 1 / − / A areΛ i = sin (cid:18) iπ N + 1) (cid:19) , Ψ i = (cid:26) sin ikπN + 1 (cid:27) Nk =1 , i = 1 , . . . , N. Note that all the eigenvectors Ψ i are of the same length, due to (cid:107) Ψ i (cid:107) = (cid:104) Ψ i , Ψ i (cid:105) = N (cid:88) k =1 sin ikπN + 1 = N − N (cid:88) k =1 cos 2 ikπN + 1 = N + 12so we do not normalize them.Numerical results for m = 7, β = 1 , (cid:107) u r − u (cid:107) A / (cid:107) f (cid:107) A − for β = 1 and the relative error (cid:107) u r − u (cid:107) A / (cid:107) f (cid:107) A − for β = 2. We use as input the coefficient vector of f with respect to the basis { Ψ i } Ni =1 , so thederivation of the exact solution u as well as the computation of the norms (cid:107) f (cid:107) A − , respectively (cid:107) f (cid:107) A − , is straightforward. In order to compute the approximated solution u r , we first generatethe coefficient vector of f with respect to the standard basis { δ ik } Ni,k =1 and then solve exactly thecorresponding m + β tridiagonal linear systems that originate from the fractional decomposition of t − β r βα ( t ). Randomness is with respect to the entries of the input coefficient vector.We study four different error quantities: the maximal error over the eigenvectors { Ψ i } Ni =1 , whichcoincides with the true estimate of the approximation error; the maximal error over a randomized setof 1000 f ’s, which is the numerical approach for estimating the former; the averaged error over theeigenvectors; and the averaged error over the random right-hand-side set. The last two quantitiesprovide information about the general behavior of the error and its expectation value. The mainobservation is that the errors, related to the eigenvectors set behave quite stably with respect tothe size of A , unlike the errors related to random vector input. Such “dimension-invariance” of theresults from the first class is due to the almost uniform distribution of the eigenvalues { Λ i } Ni =1 of A along the interval [0 ,
1] and that for every β -BURA function the endpoints 0 and 1 are extremepoints for the residual r βα ( t ) − t β − α (i.e., { , } ⊂ { η i } m + k +21 ). As N → ∞ , we have Λ → N → E α ( m, k ; β ). The β -BURA functions oscillate mainly close to 0 andare stable close to 1, making rapid convergence | r βα (Λ N ) − Λ β − αN | → E α ( m, k ; β ). Therefore, theplacement of the remaining spectrum { Λ i } N − i =1 of A with respect to { η i } is not significant for thisquantity. On the other hand, it is practically impossible to generate (a rescaled version of) Ψ N at random, so the randomized errors heavily depend on the placement of the whole spectrum of A with respect to the extreme points of the β -BURA function. As a result, both maximal andaveraged random errors can be anywhere in the interval between the minimal and maximal error ofthe eigenvectors. We generated various random sets of different size (e.g., 10 , ) and checked thatfor a fixed N the two errors behave stably with respect to the choice of randomness. This allowsus to conclude that the “dimension-instability” phenomenon is indeed fully due to the specifics ofthe spatial distribution of the spectrum of A .5.2. Multi-step 1-BURA approximation for α = { . , . } . The second series of numericalexperiments are devoted to the multi-step generalization of the method. The presented numericalexperiments for A = tridiag ( − . , . , − .
25) as in Section 5.1, k = { , } and α = { . , . } confirm the theoretical analysis in Section 4.4. The related results are summarized in Fig. 2. When OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 17 e rr o r h -1 α =0.25 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,7;1) 2.0E-09 4.0E-09 6.0E-09 8.0E-09 1.0E-08 1.2E-08 1.4E-08 1.6E-08 1.8E-082 e rr o r h -1 α =0.25 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,6;2) 2.5E-05 3.0E-05 3.5E-05 4.0E-05 4.5E-052 e rr o r h -1 α =0.50 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,7;1) 4.0E-08 4.5E-08 5.0E-08 5.5E-08 6.0E-08 6.5E-08 7.0E-08 7.5E-08 8.0E-082 e rr o r h -1 α =0.50 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,6;2) 3.5E-04 4.0E-04 4.5E-04 5.0E-04 5.5E-04 6.0E-04 6.5E-04 7.0E-04 7.5E-04 8.0E-042 e rr o r h -1 α =0.75 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,7;1) 5.0E-08 1.0E-07 1.5E-07 2.0E-072 e rr o r h -1 α =0.75 Maximal basis errorAveraged basis errorMaximal random errorAveraged random errorE (7,6;2)
Figure 1.
1D numerical validation of Lemma 2.1 for m = 7. Left: (7 ,
7; 1). Right:(7 ,
6; 2). Error is measured as indicated in (12). e rr o r h -1 α =0.5, m=k=5 Maximal basis errorAveraged basis errorMaximal random errorAveraged random error ε k(A)+2 ε k(A) E (5,5;1) e rr o r h -1 α =0.5, m=k=7 Maximal basis errorAveraged basis errorMaximal random errorAveraged random error ε k(A)+2 ε k(A) E (7,7;1) e rr o r h -1 α =0.75, m=k=5 (5,5;1) e rr o r h -1 α =0.75, m=k=7 (7,7;1) Figure 2.
1D numerical error analysis for the multi-step case. The relative errors (cid:107) u r − u (cid:107) A / (cid:107) f (cid:107) A − are plotted. OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 19 α = 0 .
5, we study the two step procedure based on α = α = 0 .
25. When α = 0 .
75, we investigateboth the two step procedure with ( α , α ) = (0 . , .
25) and the three step procedure based on α = α = α = 0 .
25. In the latter case, it is straightforward to derive the three-step analogousformula to (29), which in the particular setup implies(30) E .
25 0 .
25 0 . ( k, k ; 1) = E . k ( A ) + 3 E . k / ( A ) + 3 E . k / ( A ) . Again, as in (29), we use the short notation E . for E . ( k, k ; 1).We set h − = N +1 = 2 i , i ∈ { , , . . . , } and plot various errors. Namely, for α = 0 . { Ψ i } Ni =1 , the averaged two-step error over the eigenvectors,the maximal and averaged two-step errors over a thousand randomly generated vectors, and theone-step error E . ( k, k ; 1). For α = 0 .
75 we plot the theoretically estimated two and three steperrors (29)-(30), the true two and three step error estimates, the corresponding averaged errors ofrandom right-hand-side data, and the one-step error E . ( k, k ; 1).From the plots in Fig. 2 we observe that when the α i ’s coincide the true multi-step error estimatereaches the theoretical bound for particular sizes of A . This happens when Λ hits an extreme pointof r . , i.e., | r . (Λ ) − Λ / | ≈ E . , as for h = 2 − , k = 5 and h = { − , − } , k = 7. When α (cid:54) = α we confirm that the theoretical bound E .
25 0 . is an overestimation of the true maximalerror, since the sets of internal extreme points for r . and r . are disjoint and it is not possiblefor Λ to simultaneously hit both. Note that Λ tends to zero as N → ∞ but it never reaches zeroand the heavy oscillations of the residual in this area do not allow the maximal basis error to reach E .
25 0 . even for A of size 1023 × f . This is due to the specifics of the multi-step procedure andthe existence of pole at zero for the product rational approximation. Hence, whenever (cid:104) f , Ψ (cid:105) (cid:54) = 0this component dominates the overall error value. Because of that, the averaged basis error do notprovide reliable information about the error in the general (worst) case, since the eigenvectors of A are mutually orthogonal. This error remains substantially below the one-step error E . in allconducted experiments.As k increases, the two-step error remains better than the one-step error for a larger set of matrixsizes. For α = 0 . k = 5 the two-step error overpasses E . for h = 2 − , while for k = 7 thishappens for h = 2 − . For α = 0 .
75 the benefits of the two-step process are bigger, as the two-step error remains in vicinity of E . even for h = 2 − . However, as in the theoretical analysis,we clearly see the dimension-dependence of the multi-step errors. Nevertheless, with respect tocontrolling the ratio (cid:107) u r − u (cid:107) A / (cid:107) f (cid:107) A − in the cases when r α cannot be numerically computed, theproposed two-step procedure seems a better asymptotic choice than r α , since (cid:107) r α ( A ) A − f − A − α f (cid:107) A (cid:107) f (cid:107) A − ≤ (cid:107) r α ( A ) A − f − A − α f (cid:107) A (cid:107) f (cid:107) A − k( A ) ≤ E α ( k, k ; 2)k( A )and, comparing to (29), we experimentally observe that E α ( k, k ; 1) E α ( k, k ; 1) < E α ( k, k ; 2) if α + α = α .5.3. Comparison BURA and the method of Bonito and Pasciak, [3] . In this Subsectionwe experimentally compare the numerical efficiency of the BURA solver with the one, developedin [3] on a test example taken from their paper [3]. We consider the problem(31) ( − ∆) α u = f, u ∂ Ω = 0 , Ω = [0 , × [0 , h = 1 / ( N + 1).This leads to a 5-point stencil approximation A of − ∆ of the form A = h − tridiag ( − I N , A i,i , − I N ) , A i,i = tridiag ( − , , − , ∀ i = 1 , . . . , N. α = 0 . k = 9 α = 0 . k = 8 α = 0 . k = 7 Figure 3. A − α (cid:101) f for h = 2 − Then we have the algebraic problem (10) with A = h A /
8, where A is an N × N SPD matrixwith spectrum in the interval (0 ,
1] and f is the vector of the values of f ( x, y ) at the grid pointsscaled by h / f we use the checkerboard function on Ω \ ∂ Ω = (0 , × (0 , f ( x, y ) = (cid:26) , if ( x − . y − . > , − , otherwise . In [3, Remark 3.1] it is observed that in order to balance all the three exponential terms in theirerror estimate, the optimal quadrature approximate of A − α (cid:101) f is u Q = 2 k (cid:48) sin( πα ) π M (cid:88) (cid:96) = − m e α − (cid:96)k (cid:48) (cid:16) e − (cid:96)k (cid:48) I + A (cid:17) − f , m = (cid:24) π αk (cid:48) (cid:25) , M = (cid:24) π − α ) k (cid:48) (cid:25) , where k (cid:48) > u Q can be trivially estimated via M + m + 1 ≥ k Q + 1 , k Q := π α (1 − α ) k (cid:48) . For the ( k, k ) 1-BURA approximation u r we need to solve k + 1 linear systems. Note that inboth approaches all systems correspond to positive diagonal shifts of A , thus they possess similarcomputational complexity. Therefore, in order to perform comparison analysis on the numericalefficiency of the two solvers we need to take k Q ∼ k .As a reference solution u ref for (31) we consider the solution u Q for h = 2 − with k (cid:48) = 1 / − ) error, see [3, Table 3].In the numerical experiments we use the following parameters: h = 2 − ≈ − and k = { , , } for α = { . , . , . } , respectively. The corresponding 1-BURA-approximations of A − α (cid:101) f areillustrated on Fig. 3. Furthermore, we restrict our analysis to integer k Q . Note that both k Q andthe positive shifts e − (cid:96)k (cid:48) are continuous functions of k (cid:48) , meaning that there is a whole interval ofvalues k (cid:48) leading to the same number of systems to be solved for u Q and each k (cid:48) gives rise todifferent shift parameters, thus different quadrature rule, respectively approximation error. For usit is not clear which choice of k (cid:48) will lead to the smallest (cid:107) u − u Q (cid:107) / (cid:107) (cid:101) f (cid:107) .Straighforward computations for the considered three choices of α and u Q imply (cid:100) (1 − α ) k Q (cid:101) + (cid:100) αk Q (cid:101) + 1 = (cid:26) k Q + 1 + (cid:100) k Q (mod 4) (cid:101) , α = { . , . } ; k Q + 1 + (cid:100) k Q (mod 2) (cid:101) , α = 0 . . Therefore, for α = 0 . u Q can never consist of k + 1 = 10 summands, like u r ; for α = 0 . k Q = { , } lead to k + 1 = 9 linear systems for u Q ; for α = 0 .
75 only k Q = 6 leads to k + 1 = 8linear systems for u Q .Relative (cid:96) errors are documented in Table 7. To get the approximate solution u Q we consider k Q = { , , } , when α = { . , . , . } . We observe that in all three cases the relative error of OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 21
Table 7.
Relative (cid:96) errors for 2-D fractional diffusion. Top: ( k, k ) 1-BURA. Bot-tom: The [3] solver. α = 0 . α = 0 . α = 0 . (cid:107) u ref − u r (cid:107) / (cid:107) (cid:101) f (cid:107) . . . (cid:107) u ref − u Q (cid:107) / (cid:107) (cid:101) f (cid:107) . . . α = 0 . k = 9 α = 0 . k = 8 α = 0 . k = 7
21 31 400.180.511.5x 10 −3 Q errorBURA error −3 Q errorBURA error −3 Q errorBURA error
Figure 4.
Relative (cid:96) errors for u Q as functions on the number of solved linear systems..the BURA approximate solution u r is smaller than the error of u Q . Furthermore, for each α wekeep on increasing k Q by one until the corresponding relative (cid:96) error of u Q becomes smaller thanthe approximate solution u r obtained for ( k, k ) 1-BURA method. For α = 0 .
25, we get k Q = 38to be the smallest such integer, meaning that we need to solve 4 times more linear systems (40compared to 10) in order to beat the numerical accuracy of BURA. For α = 0 .
5, we get k Q = 20,thus we need to solve 21 linear systems if we apply [3] instead of 9, when we apply the BURAsolver. Finally, for α = 0 .
75, we get k Q = 13 and 15 linear systems to be solved, compared to8 in the BURA case. The relative errors as functions of the number of linear systems in u Q arepresented on Figure 4.5.4.
1D and 3D numerical tests with approximate solving of A u = f by PCG method. In higher spatial dimensions, when the domain Ω in (2) is a subset of R d , d >
1, we cannot computethe exact solution u = A − α f . In general, we don’t have explicitly the eigenvalues and eigenvectorsof A . Thus, in the analysis of the numerical tests, we cannot apply error estimates of the form(12). In order to numerically validate our theoretical estimates we use the two-step procedure fromSection 4.4 through the following obvious identity(33) f = A (cid:16) A − (1 − α ) (cid:0) A − α f (cid:1)(cid:17) , that holds true for an arbitrary vector f ∈ R N . In (29) we argued that taking a product oftwo rational functions as an approximation of t β − α with β = 1 the corresponding approximationerror depends on the condition number of A . For the multivariate validation of Lemma 2.1 suchdimension-dependence is not acceptable, so we choose β = 2 here. In particular, we treat r − α r α as a 2-uniform rational approximation (this is not BURA, since the approximation error is notoptimal) for the function t − = t on the unit interval (0 , γ = 2 we deduce(34) (cid:107) r − α ( A ) r α ( A ) A − f − A − f (cid:107) A ≤ E − α,α ( k, k ; 2) (cid:107) f (cid:107) A − . To avoid additional numerical inaccuracies, it is more convenient from a computational point ofview to introduce the approximation vector f r of f in (33), i.e., f r = A (cid:16) r − α ( A ) r α ( A ) A − f (cid:17) := A u r . Then we can rewrite the estimate (34) in the form(35) (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) ≤ E − α,α ( k, k ; 2) . We estimate the two-step error E − α,α ( k, k ; 2) with the help of (27), as in (28). The function r − α r α has no poles in [0 , A : E − α,α ( k, k ; 2) = max t ∈ [0 , | r − α ( t ) r α ( t ) − t | = max t ∈ [0 , | t α ε α ( t ) + t − α ε − α ( t ) + ε α ( t ) ε − α ( t ) |≤ E − α ( k, k ; 1) + E α ( k, k ; 1) + E − α ( k, k ; 1) E α ( k, k ; 1) . In practice, however, we observe that the residuals ε α and ε − α are negative and monotonicallydecreasing in [0 . , ε α (1) = − E α ( k, k ; 1), ε − α (1) = − E − α ( k, k ; 1), we can improve the two-step error estimate and conclude that:(36) (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) ≤ E − α + E α − E − α E α . The last estimate has been numerically confirmed as sharp for k = { , } . Inspecting closelythe above proof, we realize that the two-step residual is of order E − α E α or lower around zero.Furthermore, unlike the equioscillation BURA setting for the one-step process, the two-step residualreaches its maximum in absolute value only at t = 1 and the amplitudes of its other oscillationsgradually decrease as t approaches zero. As a result, the numerically computed values for themaximal and averaged random errors w.r.t. (36) are closer to E − α E α than to E − α + E α , meaningthat the typical error is significantly smaller than the worst case scenario E − α,α .We investigate two particular choices for f , namely f = (1 , . . . ,
1) and f = (1 , , . . . , A remains tridiag ( − . , . , − .
25) and as before we exactly solve the corrspondinglinear systems. In 3D, we use the finite element method in space with linear conforming tetrahedralfinite elements and an algebraic multigrid (AMG) preconditioner in the PCG solutions of thecorresponding linear systems. To be more precise, the BoomerAMG implementation, e.g. [13], isutilized in the presented numerical tests. We consider Ω = [0 , and A to be the stiffness matrixfrom the FE discretization of the problem (2) with a = a ( x ) I , with I the identity matrix in R d and a ( x ) is a piece-wise constant function in Ω. In this case, the jump of the coefficient a ( x ) isintroduced via the scaling factor 0 < µ ≤ f ’s comes from the 1D case. Since for i = 1 , . . . , N (cid:104) Ψ i , f (cid:105) = (cid:26) , i is evencot ( iπh/ , i is odd (cid:104) Ψ i , f (cid:105) = sin( iπh ) = Ψ ,i , the decompositions of the two vectors with respect to the eigen-vectors { Ψ i } Ni =1 are f = (cid:88) i is odd h cot( iπh/ Ψ i = (cid:88) i is even iπ iπh/ iπh/ Ψ i ; f = N (cid:88) i =1 h sin( iπh ) Ψ i . (37)Therefore, from x/ tan( x ) < , π/ f -decomposition(37) rapidly decay as i increases, meaning that the Ψ component dominates and the two-stepresidual at Λ determines the behavior of the error ratio (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) . To summarize(38) (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) ≈ (cid:12)(cid:12) r − α (Λ ) r α (Λ ) − Λ (cid:12)(cid:12) −−−→ h → E − α E α . Such asymptotic behavior of the (cid:96) -norm error ratio of f is numerically confirmed by the con-ducted 1D and 3D numerical experiments with k = { , } and α = { . , . } , as illustrated onthe left of Fig. 5. In 3D we run simulations up to h = 2 − , which corresponds to N = 6( h − + 1) ,while in 1D we go up to h = 2 − , and corresponding number of degrees of freedom N = h − + 1.Clearly, the error tends to E − α E α as Λ →
0. We observe that in the 3D case with homogeneouscoefficient ( L is the Laplacian) the error for mesh-size h − ∈ [2 , ] mimics the 1D error for mesh-size h − ∈ [2 , ]. For the heterogeneous 3D case, when in a half the unit cube the diffusioncoefficients are scaled by µ = 10 − , k( A ) is increased (approximately) by a factor of µ − , implyingthat Λ is closer to zero than in the homogeneous case. As a result, the 3D error for the mesh-size h − ∈ [2 , ] mimics the 1D error for the mesh-size h − ∈ [2 , ]. OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 23 e rr o r h -1 A(A -0.25 (A -0.75 f)), m=k=5, f=(1,1,...,1)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (5,5;1)*E (5,5;1) e rr o r h -1 A(A -0.25 (A -0.75 f)), m=k=5, f=(1,0,...,0)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (5,5;1)*E (5,5;1) e rr o r h -1 A(A -0.25 (A -0.75 f)), m=k=7, f=(1,1,...,1)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (7,7;1)*E (7,7;1) e rr o r h -1 A(A -0.25 (A -0.75 f)), m=k=7, f=(1,0,...,0)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (7,7;1)*E (7,7;1) e rr o r h -1 A(A -0.5 (A -0.5 f)), m=k=5, f=(1,1,...,1)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (5,5;1)*E (5,5;1) e rr o r h -1 A(A -0.5 (A -0.5 f)), m=k=5, f=(1,0,...,0)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (5,5;1)*E (5,5;1) e rr o r h -1 A(A -0.5 (A -0.5 f)), m=k=7, f=(1,1,...,1)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (7,7;1)*E (7,7;1) e rr o r h -1 A(A -0.5 (A -0.5 f)), m=k=7, f=(1,0,...,0)
1D error, µ =13D error, µ =13D error, µ =10 -3 E (7,7;1)*E (7,7;1) Figure 5.
1D and 3D numerical error analysis. Left: f = (1 , , . . . , f = (1 , , . . . , (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) are plotted. The coefficients in the f -decomposition (37) have symmetry due to the relations 2 h sin( iπh ) =2 h sin(( N + 1 − i ) πh with h = 1 / ( N + 1). Unlike for f case, the contribution of Ψ here isnegligible and the behavior of the relative error (cid:107) f r − f (cid:107) / (cid:107)A − f (cid:107) is dominated by the error | r − α (0 . r α (0 . − . | . This effect weakens with h →
0, because the coefficients depend onthe mesh-size and their distribution spreads away (standard deviation increases) when the grid isrefined. The two-step residual is stable around t = 1 / h − .The numerical results perfectly agree with this argument and, similar to the 1D–3D correspondencefor f , we observe that the slope of the error decay is steeper in 3D and the error decreases in thecase of piece-wise constant coefficient a ( x ).In the presented numerical tests, as a stopping criteria for the BoomerAMG PCG solver we haveused a relative error less or equal to 10 − . However, we want to note that the numerical results arepractically not affected by using stopping criteria 10 − , instead. Furthermore, for precision 10 − the order of applying the 1-BURA functions r . and r . (i.e., taking α = 0 .
25 or α = 0 .
75 first)seems irrelevant and the corresponding relative errors have the same first five meaningful digits.This implies that the main numerical difficulties are related to the performance of Algorithm 3.1and the correctness of the subsequent representation of r βα as a sum of partial fractions.6. Concluding remarks
In this paper we propose algorithms of optimal complexity for solving the linear algebraic system A α u = f , 0 < α <
1, where A is a sparse SPD matrix. The target class of applied problems A obtained by a finite difference or finite element discretization of a second order elliptic problem. Ourmain assumption is that the system A u = f can be solved with optimal computational complexity,e.g. by multi-grid, multi-level or other efficient solution technique. The proposed in the papermethod is applicable also when the matrix is not given explicitly, but one has at hand an optimalsolution procedure for the linear system A u = f and a upper bound for the spectrum of A .The method is based on best uniform rational approximations (BURA) of t β − α for 0 ≤ t ≤ β . Bigger β means stronger regularity assumptions and this is the reason to concentrateour considerations mostly to the cases β ∈ { , } . Depending on α , β and the degree k of bestuniform rational approximation a relative accuracy of the method between O (10 − ) and O (10 − )can be obtained for k ∈ { , , } . Then solution of A α u = f reduces to solving k + β problems withsparse SPD matrices of the form A + c I , c ≥ α = 0 . α ∈ (0 , m = k , in the case of smaller α . For larger α , the multi-step algorithm has some promising features. Future theoretical and experimentalinvestigations are needed for better understanding the observed superior convergence of two-stepBURA when α (cid:54) = α .The method and the integral quadrature formula method from [3] have been experimentallycompared. The test setup has been taken from [3, Section 4.1] with h = 2 − ≈ − . The BURAmethod performs better in all numerical experiments and this effect increases as α decreases. Acknowledgement
This research has been partially supported by the Bulgarian National Science Fund under grantNo. BNSF-DN12/1. The work of R. Lazarov has been partially supported by the grant NSF-DMS
OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 25
References [1] O. G. Bakunin.
Turbulence and diffusion: scaling versus equations . Springer Science & Business Media, 2008.[2] P. W. Bates. On some nonlocal evolution equations arising in materials science.
Nonlinear dynamics and evolutionequations , 48:13–52, 2006.[3] A. Bonito and J. Pasciak. Numerical approximation of fractional powers of elliptic operators.
Mathematics ofComputation , 84(295):2083–2110, 2015.[4] A. Bonito and J. Pasciak. Numerical approximation of fractional powers of regularly accretive operators.
IMA JNumer Anal , 37(3):1245–1273, 2017.[5] L. Chen, R. Nochetto, O. Enrique, and A. J. Salgado. Multilevel methods for nonuniformly elliptic operatorsand fractional diffusion.
Mathematics of Computation , 85:2583–2607, 2016.[6] E. W. Cheney and M. J. D. Powell. The differential correction algorithm for generalized rational functions.
Constructive Approximation , 3(1):249–256, 1987.[7] V. Druskin and L. Knizhnerman. Extended Krylov subspaces: approximation of the matrix square root andrelated functions.
SIAM Journal on Matrix Analysis and Applications , 19(3):755–771, 1998.[8] C. B. Dunham. Difficulties in rational Chebyshev approximation. In
Conference on Constructive Theory ofFunctions, Varna, Bulgaria , pages 319–327, 1984.[9] I. Gavrilyuk, W. Hackbusch, and B. Khoromskij. Hierarchical tensor-product approximation to the inverse andrelated operators for high dimensional elliptic problems.
Computing , 74(2):131–157, 2005.[10] G. Gilboa and S. Osher. Nonlocal operators with applications to image processing.
Multiscale Modeling & Sim-ulation , 7(3):1005–1028, 2008.[11] S. Harizanov and S. Margenov. Positive approximations of the inverse of fractional powers of SPD M-matrices.submitted, posted as arXiv:1706.07620v1, June 2017.[12] S. Harizanov, S. Margenov, P. Marinov, and Y. Vutov. Volume constrained 2-phase segmentation method utilizinga linear system solver based on the best uniform polynomial approximation of x − / . Journal of Computationaland Applied Mathematics , 310:115–128, 2017.[13] V. E. Henson and U. M. Yang. Boomeramg: A parallel algebraic multigrid solver and preconditioner.
AppliedNumerical Mathematics , 41(1):155 – 177, 2002.[14] N. J. Higham. Stable iterations for the matrix square root.
Numerical Algorithms , 15(2):227–242, 1997.[15] M. Ili´c, I. W. Turner, and V. Anh. A numerical solution using an adaptively preconditioned Lanczos methodfor a class of linear systems related with the fractional Poisson equation.
International Journal of StochasticAnalysis , 2008, 2009.[16] C. Kenney and A. J. Laub. Rational iterative methods for the matrix sign function.
SIAM J. Matrix Anal. Appl. ,12(2):273–291, Mar. 1991.[17] A. Kilbas, H. Srivastava, and J. Trujillo.
Theory and Applications of Fractional Differential Equations . Elsevier,Amsterdam, 2006.[18] R. Lazarov and P. Vabishchevich. A numerical study of the homogeneous elliptic equation with fractional orderboundary conditions.
Fractional Calculus and Applied Analysis , 20(2):337–351, 2017.[19] P. G. Marinov and A. S. Andreev. A modified Remez algorithm for approximate determination of the rationalfunction of the best approximation in Hausdorff metric.
Comptes rendus de l’Academie bulgare des Scieces ,40(3):13–16, 1987.[20] M. Matsuki and T. Ushijima. A note on the fractional powers of operators approximating a positive definiteselfadjoint operator.
J. Fac. Sci. Univ. Tokyo Sect. IA Math. , 40(2):517–528, 1993.[21] B. McCay and M. Narasimhan. Theory of nonlocal electromagnetic fluids.
Archives of Mechanics , 33(3):365–384,1981.[22] R. Metzler, J.-H. Jeon, A. G. Cherstvy, and E. Barkai. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking.
Physical Chemistry ChemicalPhysics , 16(44):24128–24164, 2014.[23] S. Nepomnyaschikh. Mesh theorems on traces, normalizations of function traces and their inversion.
Sov. J.Numer. Anal. Math. Modelling , 6(3):223–242, 1991.[24] E. Saff and H. Stahl. Asymptotic distribution of poles and zeros of best rational approximants to x α on [0, 1]. Sta , 299(4):2, 1992.[25] S. A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces.
Journal of the Mechanicsand Physics of Solids , 48(1):175–209, 2000.[26] H. R. Stahl. Best uniform rational approximation of x α on [0 , Acta Mathematica , 190(2):241–306, 2003.[27] P. N. Vabishchevich. Numerical solving the boundary value problem for fractional powers of elliptic operators.
CoRR , abs/1402.1636, 2014.[28] P. N. Vabishchevich. Numerically solving an equation for fractional powers of elliptic operators.
Journal ofComputational Physics , 282:289–302, 2015.[29] R. S. Varga and A. J. Carpenter. Some numerical results on best uniform rational approximation ofx α on [0,1]. Numerical Algorithms , 2(2):171–185, 1992.[30] J. Xu and L. Zikatanov. The method of alternating projections and the method of subspace corrections in Hilbertspace.
Journal of the American Mathematical Society , 15(3):573–597, 2002. [31] G. M. Zaslavsky. Chaos, fractional kinetics, and anomalous transport.
Physics Reports , 371(6):461–580, 2002.[32] X. Zhao, X. Hu, W. Cai, and G. E. Karniadakis. Adaptive finite element method for fractional differentialequations using hierarchical matrices.
Computer Methods in Applied Mechanics and Engineering , 325(SupplementC):56 – 76, 2017. AppendixTable 8.
The coefficients in the representation (20) of 1-BURA P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , j α = 0 . α = 0 . α = 0 . c j d j c j d j c j d j Table 9.
The coefficients in the representation (20) of 2-BURA P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , j α = 0 . α = 0 . α = 0 . c j d j c j d j c j d j Table 10.
The coefficients in the representation (20) of 2-BURA P ∗ ( t ) /Q ∗ ( t ) of t − α on [0 , j α = 0 . α = 0 . α = 0 . c j d j c j d j c j d j Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, Acad.G. Bonchev, bl. 25A, 1113 Sofia, Bulgaria ([email protected])Deptartment of Mathematics, Texas A&M University, College Station, TX 77843-3368, USA ([email protected])and Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Acad. G. Bonchev,bl. 8, 1113 Sofia, BulgariaInstitute of Information and Communication Technologies, Bulgarian Academy of Sciences, Acad.G. Bonchev, bl. 25A, 1113 Sofia, Bulgaria ([email protected])Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, Acad.G. Bonchev, bl. 25A, 1113 Sofia, Bulgaria ([email protected])
OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 27
Table 11.
Numerical error ε in (12) for ( m, k ) = (7 ,
7) and β = 1. For each α the leftcolumn shows the results for f consisting of eigenvectors, while the right column is the errorfor f taken as 1000 random eigenvector combinations. In each box averaged error (top) andthe maximal error (bottom) are reported. h − α = 0 . α = 0 . α = 0 . E α (7 ,
7; 1) =
E-06 E α (7 ,
7; 1) =
E-05 E α (7 ,
7; 1) =
E-04 { Ψ i } Ni =1 rand1000 { Ψ i } Ni =1 rand1000 { Ψ i } Ni =1 rand10008 1.9565E-06 1.7915E-06 2.8431E-05 2.7493E-05 5.6995E-04 6.9670E-043.2061E-06 2.9876E-06 4.6024E-05 4.5356E-05 7.6989E-04 7.6437E-0416 2.0032E-06 2.4354E-06 2.6759E-05 2.3804E-05 5.2532E-04 6.5826E-043.1948E-06 2.9711E-06 4.5812E-05 3.7222E-05 7.8518E-04 7.0151E-0432 2.1362E-06 2.7011E-06 2.9166E-05 3.3229E-05 5.0559E-04 5.9927E-043.2522E-06 2.9153E-06 4.5720E-05 3.7850E-05 7.8286E-04 6.4667E-0464 2.0616E-06 1.9800E-06 2.9487E-05 3.9642E-05 4.9911E-04 4.8683E-043.2564E-06 2.8259E-06 4.6035E-05 4.3436E-05 7.8744E-04 6.3295E-04128 2.0567E-06 1.5541E-06 2.9505E-05 4.0691E-05 4.9672E-04 3.5982E-043.2565E-06 2.6089E-06 4.6033E-05 4.4048E-05 7.8922E-04 6.7316E-04256 2.0675E-06 1.8786E-06 2.9370E-05 3.8431E-05 4.9953E-04 4.4143E-043.2566E-06 2.6313E-06 4.6029E-05 4.3061E-05 7.8959E-04 7.0087E-04512 2.0734E-06 2.6468E-06 2.9312E-05 3.0132E-05 5.0116E-04 6.3817E-043.2566E-06 3.0750E-06 4.6036E-05 3.9441E-05 7.8965E-04 7.4051E-041024 2.0736E-06 2.3078E-06 2.9288E-05 2.3497E-05 5.0152E-04 6.3749E-043.2566E-06 2.9537E-06 4.6037E-05 3.7366E-05 7.8965E-04 7.1731E-04 Table 12.
Numerical error ε in (12) for ( m, k ) = (7 ,
6) and β = 2. For each α the leftcolumn shows the results for f consisting of eigenvectors, while the right column is the errorfor f taken as 1000 random eigenvector combinations. In each box averaged error (top) andthe maximal error (bottom) are reported. h − α = 0 . α = 0 . α = 0 . E α (7 ,
6; 2) =
E-08 E α (7 ,
6; 2) =
E-08 E α (7 ,
6; 2) =
E-07 { Ψ i } Ni =1 rand1000 { Ψ i } Ni =1 rand1000 { Ψ i } Ni =1 rand10008 1.0632E-08 1.0631E-08 5.9400E-08 6.0261E-08 1.0966E-07 1.0468E-071.5547E-08 1.2642E-08 7.8416E-08 6.9105E-08 2.01E00-07 1.5912E-0716 9.5732E-09 1.3119E-08 5.4368E-08 5.3487E-08 1.1609E-07 1.7237E-071.7828E-08 1.3955E-08 7.7768E-08 7.6782E-08 2.0032E-07 1.8421E-0732 1.1235E-08 3.6625E-09 5.1080E-08 7.5291E-08 1.2520E-07 6.2571E-081.8049E-08 1.5797E-08 7.8585E-08 7.8456E-08 2.0100E-07 1.9484E-0764 1.1198E-08 3.2757E-09 5.1082E-08 7.7033E-08 1.2412E-07 2.8021E-081.8040E-08 1.0433E-08 7.8644E-08 7.8491E-08 2.0090E-07 9.7132E-08128 1.1468E-08 1.6223E-08 5.0090E-08 6.7111E-08 1.2665E-07 5.8966E-081.8066E-08 1.7177E-08 7.8573E-08 7.7379E-08 2.0110E-07 1.4951E-07256 1.1446E-08 3.8759E-09 4.9977E-08 4.4867E-08 1.2795E-07 1.8929E-071.8062E-08 1.7213E-08 7.8638E-08 6.9319E-08 2.0111E-07 2.0039E-07512 1.1473E-08 1.1522E-08 4.9994E-08 3.1754E-08 1.2783E-07 6.5736E-081.8065E-08 1.2699E-08 7.8647E-08 7.0909E-08 2.0111E-07 1.9958E-07 Institute of Information and Communication Technologies, Bulgarian Academy of Sciences, Acad.G. Bonchev, bl. 25A, 1113 Sofia, Bulgaria ([email protected])
Table 13.
Numerical error ε in (12) for the multi-step case ( m, k ) = (5 ,
5) and β = 1.For each setup the left column shows the results for f consisting of eigenvectors, while theright column is the error for f taken as 1000 random eigenvector combinations. In each boxaveraged error (top) and the maximal error (bottom) are reported. h − α = 0 . α = 0 . α = 0 . u r = A − . (cid:0) A − . f (cid:1) u r = A − . (cid:0) A − . (cid:0) A − . f (cid:1) (cid:1) u r = A − . (cid:0) A − . f (cid:1) { Ψ i ) } Ni =1 rand1000 { Ψ i ) } Ni =1 rand1000 { Ψ i ) } Ni =1 rand100016 4.1065E-05 4.1827E-05 8.7666E-05 1.2954E-04 1.9972E-04 1.7173E-049.4745E-05 7.2813E-05 2.3891E-04 2.2187E-04 3.9098E-04 3.3170E-0432 4.7118E-05 5.1558E-05 1.0939E-04 2.0585E-04 2.1337E-04 2.3155E-041.2093E-04 1.0361E-04 4.8074E-04 3.9256E-04 5.0294E-04 4.1310E-0464 5.1447E-05 1.0267E-04 1.4007E-04 7.7816E-04 2.3385E-04 4.5311E-042.0432E-04 1.5404E-04 1.1396E-03 9.9375E-04 5.6349E-04 5.4832E-04128 5.7422E-05 3.2390E-04 2.0400E-04 4.1774E-03 2.3648E-04 4.0115E-044.5572E-04 4.3974E-04 6.2524E-03 6.0171E-03 5.9066E-04 5.4683E-04256 5.9104E-05 4.3670E-04 2.3325E-04 7.4754E-03 2.3848E-04 7.7769E-045.3831E-04 5.2302E-04 1.0097E-02 9.7472E-03 1.0680E-03 1.0328E-03512 5.9310E-05 5.9445E-04 2.4906E-04 1.3916E-02 2.4811E-04 3.8679E-037.7978E-04 7.5465E-04 1.9729E-02 1.9016E-02 5.8996E-03 5.6791E-031024 5.9021E-05 4.1906E-04 2.4148E-04 9.6104E-03 2.5064E-04 5.3558E-038.0305E-04 7.2124E-04 1.9733E-02 1.7263E-02 6.8130E-03 6.5918E-03 Table 14.
Numerical error ε in (12) for the multi-step case ( m, k ) = (7 ,
7) and β = 1.For each setup the left column shows the results for f consisting of eigenvectors, while theright column is the error for f taken as 1000 random eigenvector combinations. In each boxaveraged error (top) and the maximal error (bottom) are reported. h − α = 0 . α = 0 . α = 0 . u r = A − . (cid:0) A − . f (cid:1) u r = A − . (cid:0) A − . (cid:0) A − . f (cid:1) (cid:1) u r = A − . (cid:0) A − . f (cid:1) { Ψ i ) } Ni =1 rand1000 { Ψ i ) } Ni =1 rand1000 { Ψ i ) } Ni =1 rand10008 5.0192E-06 5.1628E-06 1.0010E-05 1.1877E-05 3.6246E-05 4.0301E-059.8128E-06 9.6002E-06 2.5169E-05 2.4547E-05 7.0309E-05 6.8761E-0516 5.8605E-06 1.4026E-05 1.5501E-05 6.6467E-05 3.4149E-05 3.7829E-051.9834E-05 1.9260E-05 9.7937E-05 9.4949E-05 6.6285E-05 5.2376E-0532 6.7145E-06 2.0145E-05 2.0351E-05 1.2619E-04 3.9973E-05 8.4991E-052.5714E-05 2.5005E-05 1.7685E-04 1.7140E-04 1.1430E-04 1.1152E-0464 6.3275E-06 1.8390E-05 1.9405E-05 1.4991E-04 4.3354E-05 1.6584E-042.6754E-05 2.3775E-05 1.8260E-04 1.7894E-04 2.2211E-04 2.1671E-04128 6.3162E-06 1.7655E-05 1.9982E-05 1.9293E-04 4.5322E-05 2.5209E-042.7195E-05 2.1459E-05 2.3671E-04 2.3119E-04 3.2535E-04 3.1700E-04256 6.5595E-06 3.5085E-05 2.4035E-05 5.9986E-04 4.5112E-05 2.3879E-044.5202E-05 4.4046E-05 8.6567E-04 8.3841E-04 3.1927E-04 2.9358E-04512 6.7797E-06 8.1866E-05 2.9931E-05 2.1139E-03 4.6091E-05 4.4815E-041.1361E-04 1.0915E-04 3.1081E-03 2.9789E-03 5.7114E-04 5.5602E-041024 6.8055E-06 9.0098E-05 3.1472E-05 3.0128E-03 4.7830E-05 1.1992E-031.1366E-04 1.0352E-04 3.9165E-03 3.8048E-03 1.6865E-03 1.6425E-03 OLUTION OF SYSTEMS WITH FRACTIONAL POWER OF SPD MATRICES 29
Table 15.
Numerical (cid:96) -error ε for h − = 2 n based on (36) for f reconstruction with( m, k ) = (5 ,
5) and β = 1. For each setup the left column shows the results for f consistingof eigenvectors, while in the middle and right columns are the errors for f taken as 1000random eigenvector combinations. In each box averaged error (top) and the maximal error(bottom) are reported. h − A ( A − . ( A − . f )) A ( A − . ( A − . f )) E . + E . − E . E . = E-03 2 E . − E . = E-04 { Ψ i ) } Ni =1 rand1000 rand1000 { Ψ i ) } Ni =1 rand1000 rand100016 4.7354E-04 7.5219E-05 7.5147E-05 1.0598E-04 8.7507E-06 8.6706E-062.4855E-03 8.1696E-05 8.2122E-05 4.8149E-04 1.3993E-05 1.5610E-0532 4.7350E-04 2.4231E-05 2.4240E-05 1.0689E-04 4.0082E-06 3.9951E-062.6751E-03 2.5452E-05 2.5517E-05 5.1127E-04 6.8881E-06 6.7313E-0664 4.7775E-04 3.4729E-06 3.4744E-06 1.0850E-04 6.9353E-06 6.9349E-062.7268E-03 4.0236E-06 3.9443E-06 5.3095E-04 7.3211E-06 7.2493E-06128 4.7955E-04 6.0366E-07 6.0448E-07 1.0878E-04 6.2346E-06 6.2351E-062.7402E-03 8.7980E-07 9.4323E-07 5.3609E-04 6.2777E-06 6.2911E-06256 4.8034E-04 5.7357E-07 5.7343E-07 1.0891E-04 2.5050E-06 2.5049E-062.7437E-03 5.8741E-07 5.8982E-07 5.3740E-04 2.5149E-06 2.5144E-06512 4.8074E-04 1.2493E-06 1.2493E-06 1.0895E-04 6.7062E-07 6.7059E-072.7445E-03 1.2496E-06 1.2496E-06 5.3773E-04 6.7463E-07 6.7458E-071024 4.8094E-04 2.0673E-07 2.0677E-07 1.0898E-04 6.9352E-07 6.9357E-072.7447E-03 2.0885E-07 2.0924E-07 5.3781E-04 6.9625E-07 6.9637E-072048 4.8104E-04 4.1620E-07 4.1620E-07 1.0899E-04 2.8010E-07 2.8010E-072.7448E-03 4.1646E-07 4.1649E-07 5.3783E-04 2.8034E-07 2.8027E-074096 4.8109E-04 5.3275E-07 5.3275E-07 1.0900E-04 3.6167E-08 3.6185E-082.7448E-03 5.3281E-07 5.3282E-07 5.3784E-04 3.7534E-08 3.7333E-08 Table 16.
Numerical (cid:96) -error ε for h − = 2 n based on (36) for f reconstruction with( m, k ) = (7 ,
7) and β = 1. For each setup the left column shows the results for f consistingof eigenvectors, while in the middle and right columns are the errors for f taken as 1000random eigenvector combinations. In each box averaged error (top) and the maximal error(bottom) are reported. h − A ( A − . ( A − . f )) A ( A − . ( A − . f )) E . + E . − E . E . = E-04 2 E . − E . = E-05 { Ψ i ) } Ni =1 rand1000 rand1000 { Ψ i ) } Ni =1=1