Arndt - Algorithms for Programmers (523138), страница 50
Текст из файла (страница 50)
The square root hasto be chosen such that all its eigenvalues all have positive real parts. The sign decomposition is undefinedif A has eigenvalues on the imaginary axis. S is its own inverse (its eigenvalues are ±1). UseN0= 1Nk+1=Sk+1=S0 = Aµ¶1 − Sk2Nk 1 +2Sk Nk+1 → S(13.95a)→N(13.95b)(13.95c)While we are at it: Define a matrix A+ asA+¡¢ ¡¢:= (AAT )−1 AT = V Ω−2 V TV ΩU T = V Ω−1 U TThis looks suspiciously like the inverse of A.
In fact, this is the pseudo-inverse of A:¡¢ ¡¢A+ A = V Ω−1 U TU Ω V T = 1 but wait(13.96)(13.97)A+ has the nice property to exist even if A−1 does not. If A−1 exists, it is identical to A+ . If not,A+ A 6= 1 but A+ will give the best possible (in a least-square sense) solution x+ = A+ b of the equationA x = b (see [21], p.770ff). To find (AAT )−1 use the iteration for the inverse:¡¢Φ(x) = x 1 + (1 − dx) + (1 − dx)2 + . . .(13.98)2with d = A AT and the start value x0 = 2 − n (A AT )/ ||A AT || where n is the dimension of A.CHAPTER 13.
ARITHMETICAL ALGORITHMS13.7292n-th root by Goldschmidt’s algorithmTBD: show derivation (as root of 1) TBD: give numerical exampleTBD: parallel featureThe so-called Goldschmidt algorithm to approximate the a-th root of d can be stated as follows:setE0 := da−1x0 := d(13.99a)then iterate:1 − Ek→1a:= xk · rk:= Ek · rka→1rk(13.99b):= 1 +xk+1Ek+1(13.99c)(13.99d)until x close enough tox∞The invariant quantity is(xk ·r)a(Ek ·r a ) .xa0E0=dada−1da .(13.100)Clearlyxak+1Ek+1With1==(xk · r)axak=(Ek · ra )Ek(13.101)= d and E∞ = 1, therefore xa∞ = d. Convergence is quadratic.A variant for inverse roots is as follows:setx0:= 1E0 := d(13.102)then iterate as in formulas 13.99b..13.99dFor a = 1 we get:1d=∞Y(2 − Ek )(13.103)k=0(13.104)where Ek+1 := Ek (2 − Ek ).For a = 2 we get a iteration for the inverse square root:1√d=∞Y3 − Ek2(13.105)k=0(13.106)k 2where Ek+1 := Ek ( 3−E2 ) .
Cf. [55].¡Higher order iterations are obtained by appending higher terms to the expression 1 +definitions of rk+1 as suggested by equation 13.52 (and the identification y = 1 − E):µ1 − Ek1++a(1 + a) (1 − Ek )2[third order:] +2 a2(1 + a) (2 + a) (1 − Ek )3[fourth order:] +6 a3+ ... +¶(1 + a) (1 + 2w) . . .
(1 + n a) (1 − Ek )n[(n + 1)-th order:] +n! an1−Eka¢in the(13.107)CHAPTER 13. ARITHMETICAL ALGORITHMS293For those fond of products: for d > 0, d 6= 1√¶∞ µY11+qkd =where q0 =k=0d+1,d−1qk+1 = 2qk2 − 1(13.108)hk+1 = (hk + 2)2 (hk − 1) + 1(13.109)(convergence is quadratic) and√d =∞ µYk=021+hk¶where h0 =d+3,d−1(convergence is cubic). These are given in [56], the first is ascribed to Friedrich Engel.
The paper givesQk4d2hk+1 = d−1i=0 hi − 3. Note that for 13.108qk1qk==T2k (q0 )(d − 1)NPN ¡2 N ¢i=02i(13.110)diwhere N = 2k(13.111)where Tn is the n-th Chebychev polynomial of the first kind. One findsqk=T2k (1/c) whered=1−c,1+cc<1(13.112)andrwhich can be expressed in d =1−c1+c√1−c1+c≈1 − c U2k −1 (1/c)cT2k (1/c)(13.113)asd ≈1+d2 d U2k −1 ( 1−d )1+d1 − d T2k ( 1−d)where d > 1(13.114)Qk−1where Un is the n-th Chebychev polynomial of the second kind.
Note that U2k −1 (x) = 2k i=0 T2i (x).A computation might successively compute T2i = 2 T22i−1 − 1, accumulate the product U2i −1 =12 U2i−1 −1 T2i−1 until U2k −1 and T2k are obtained. Alternatively one might use Uk (x) = k+1∂z Tk+1 (x)and the recursion for the coefficients of T (cf. section 13.14).√√Tk (i a)/ikk (i a)Similarly, approximations for a2 + 1 are given by Rk (a) = Uk−1= −i UTk−1−1.(i a) where i =(i a)/ik−1The composition law Rm n (x) = Rm (Rn (x)) holds in analogy to Tm n (x) = Tm (Tn (x))k12345Rk (1)11327517124129Rk (a)a11 + 2 a22a3 a + 4 a31 + 4 a21 + 8 a2 + 8 a 44 a + 8 a35 a + 20 a3 + 16 a51 + 12 a2 + 16 a4(13.115)CHAPTER 13.
ARITHMETICAL ALGORITHMS13.8294Iterations for the inversion of a functionIn this section we will look at general forms of iterations for zeros9 x = r of a function f (x). Iterationsare themselves functions Φ(x) that, when ‘used’ asxk+1= Φ(xk )(13.116)will make x converge towards x∞ = r if x0 was chosen not too far away from r.The functions Φ(x) must be constructed so that they have an attracting fixed point where f (x) has azero: Φ(r) = r (fixed point) and |Φ0 (r)| < 1 (attracting).The order of convergence (or simply order ) of a given iteration can be defined as follows: let x = r · (1 + e)with |e| ¿ 1 and Φ(x) = r ·(1+αen +O(en+1 ), then the iteration Φ is called linear (or first order) if n = 1(and |α| < 1) and super-linear if n > 1.
Iterations of second order (n = 2) are often called quadratically-,those of third order cubically convergent. A linear iteration improves the result by (roughly) adding aconstant amount of correct digits with every step, a super-linear iteration of order n will multiply thenumber of correct digits by n.For n ≥ 2 the function Φ has a super-attracting fixed point at r: Φ0 (r) = 0.
Moreover, an iteration oforder n ≥ 2 hasΦ0 (r) = 0,Φ00 (r) = 0,...,Φ(n−1) (r) = 0(13.117)There seems to be no standard term for this in terms of fixed points, attracting of order n might beappropriate.To any iteration of order n for a function f one can add a term f (xk )n+1 · ϕ(x) (where ϕ is an arbitraryfunction that is analytic in a neighborhood of the root) without changing the order of convergence. It isassumed to be zero in what follows.nAny two iterations of (the same) order n differ in a term (x − r) ν(x) where ν(x) is a function that isfinite at r (cf. [10], p.
174, ex.3).Two general expressions, Householder’s formula and Schröder’s formula, can be found in the literature.Both allow the construction of iterations for a given function f (x) that converge at arbitrary order. Asimple construction that contains both of them as special cases is given.TBD: p-adic iterations13.8.1Householder’s formulaLet n ≥ 2, then³Φn (xk ) :=xk + (n − 1) ³g(xk )f (xk )g(xk )f (xk )´(n−2)n+1ϕ(x)´(n−1) + f (xk )(13.118)gives a n−th order iteration for a (simple) root r of f . g(x) must be a function that is analytic near theroot and is set to 1 in what follows (cf. [10] p.169).For n = 2 we get Newton’s formula:Φ2 (x)ff0(13.119)2f f 02f 02 − f f 00(13.120)= x−For n = 3 we get Halley’s formula:Φ3 (x) =9 orroots of the function: r so thatf (r) = 0x−CHAPTER 13. ARITHMETICAL ALGORITHMS295n = 4 and n = 5 result in:Φ4 (x) =Φ5 (x) =3f (f f 00 − 2f 02 )6f f 0 f 00 − 6f 03 − f f 000³´34f 6f 0 − 6f f 0 f 00 + f 2 f 000¢x + ¡ 3 0000f f − 24f 0 4 + 36f f 0 2 f 00 − 8f 2 f 0 f 000 − 6f 2 f 00 2x−Second order 13.118 with f (x) :=that require long divisions.1xa(13.121)(13.122)− d gives formula 13.56, but for higher orders one gets iterationsKalantari and Gerlach [57] give the iterationBm (x) = x − f (x)where m ≥ 2 andDm (x)=f 0 (x) f (x)det 0 . ..f 00 (x)2!f 0 (x)f (x)...00Dm−2 (x)Dm−1 (x)...............(13.123)f (m−1) (x)(m−1)!........f (m) (x)m!f (m−1) (x)(m−1)!.f (x)...00f (x)2!(13.124)f 0 (x)(and D0 = 1).
The iteration turns out to be identical to the one of Householder. A recursive definitionfor Dm (x) is given byDm (x) =mXi−1(−1)f (x)i−1i=1f (i) (x)Dm−i (x)i!(13.125)√Similar, the well-known derivation of Halley’s formula by applying Newton’s formula to f / f 0 can begeneralized to produce m-order iterations as follows: Let F1 (x) = f (x) and for m ≥ 2 letFm−1 (x)0Fm−1 (x)1/m(13.126)Fm−1 (x)0Fm−1(x)(13.127)Fm (x)=Gm (x)= x−Then Gm (x) = Dm (x) as shown in [57].13.8.2Schröder’s formulaLet n ≥ 2, and ϕ be an arbitrary (analytic near the root) function that is set to zero in what follows,then the expressionµ¶t−1nX11f (xk )t∂+ f (xk )n+1 ϕ(x)(13.128)Φn (xk ) :=(−1)t00t!f(x)f(x)kkt=0gives a n−th order iteration for a (simple) root r of f (cf.
[9] p.13). This is, explicitly,Φn=¢f2f 3 ¡ 002f00−·f−· 3f − f 0 f 000(13.129)003051! f2! f3! f¢f4 ¡−· 15f 003 − 10f 0 f 00 f 000 + f 02 f 00004! f 07¢f5 ¡−· 105f 004 − 105f 0 f 002 f 000 + 10f 02 f 0002 + 15f 02 f 00 f 0000 − f 03 f 00000 − . . .5! f 09x−CHAPTER 13. ARITHMETICAL ALGORITHMS296The second order iteration is the same as the corresponding iteration from 13.118 while all higher orderiterations are different. The third order iteration obtained upon truncation after the third term on theright hand side, written asµ¶ff f 00Φ3 = x − 0 1 − 02(13.130)f2fis sometimes referred to as ‘Householder’s method’.Cite from [9], (p.16, translation has a typo in the first formula):If we denote the general term by−f a χaa! f 02a−1(13.131)the numbers χa can be easily computed by the recurrenceχa+1=(2a − 1)f 00 χa − f 0 ∂χa(13.132).Formula 13.128 with f (x) := 1/xa − d gives the ‘division-free’ iteration 13.53 for arbitrary order.For f (x) := log(x) − d one gets the iteration 13.9.3.For f (x) := x2 − d one getsÃΦ(x)= x−!¡ 2¢2¡ 2¢3¡¢4x −dx −d5 x2 − dx2 − d++++ ...2x8x316x5128x7(13.133)µ¶125x2 − d34= x− y+· y2 +·y+·y+...where y :=(13.134)232x(2x)(2x)2x¡¢x2 − d= x − 2x · Y + Y 2 + 2 Y 3 + 5 Y 4 + 14 Y 5 + 42 Y 6 + .
. .where Y :=(13.135)(2x)2The connection between Householder’s and Schröder’s iterations is that the Taylor series of the k-th orderHouseholder iteration around f = 0 up to order k − 1 gives the k-th order Schröder iteration.13.8.3Dealing with multiple rootsThe iterations given so far will not converge at the stated order if f has a multiple root at r. As an exampleconsider the (for simple roots second order) iteration Φ(x) = x − f /f 0 for f (x) = (x2 − d)p , p ∈ N, p ≥ 2:√√2−d2Φ? (x) = x − xp2x. Its convergence is only linear: Φ( d(1 + e)) = d(1 + p−1p e + O(e ))Householder ([10] p.161 ex.6) givesΦ2 (x) = x − p ·ff0(13.136)as a second order iteration for functions f known a priori to have roots of multiplicity p.A general approach is to use the general10 expressions with F := f /f 0 instead of f . Both F and f havethe same set of roots, but the multiple roots of f are simple roots of F .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
