metodes iteratius, Apuntes de Álgebra

¡Descarga metodes iteratius y más Apuntes en PDF de Álgebra solo en Docsity!

INTRODUCCI O´

Introducci´o

Els m`etodes de Jacobi, Gauss-Seidel i SOR

• Obtenci´o m`etodes iteratius tipus lineal:
  1. descomposici´o A = M − N , M aproximaci´o de A,
  2. equaci´o x = M −^1 (N x + b), equivalent a Ax = b
  3. a partir x(0), aproximaci´o inicial, iteraci´o x(k+1)^ = M −^1 (N x(k)^ + b).
• M −^1 (N x(k)^ + b): sistema amb matriu precondicionadora M i terme independent N x(k)^ + b.
• sistemes amb M m´es senzills que amb A per efectivitat.

Introducci´o

SOR sim`etric, SSOR

• Important A SPD ⇒ precond. SPD.
• Jacobi OK, per`o no GS o SOR.
• Simetritzem SOR
• SOR invers: precond. M = (^) ω^1 D + U
• Suposem A SPD i fem un pas de SOR seguit d’un pas de SOR invers ( U = LT^ ): x ˜(k+1)^ = (^1 ω D + L)−^1 (−((1 − 1 ω )D + LT^ )x(k)^ + b)

x(k+1)^ = ( ω^1 D + LT^ )−^1 (−((1 − (^) ω^1 )D + L)˜x(k+1)^ + b) = · · · + (^1 ω D + LT^ )−^1 2 −^ ω ω D(^1 ω D + L)−^1 ︸ ︷︷ ︸ M−^1

b

• Dedu¨ım

M = (

ω

D + L)

ω 2 − ω

D−^1 (

ω

D + LT^ )

• M SPD si A SPD.

Introducci´o

Convergencia dels metodes iteratius lineals

• M`etode iteratiu convergeix:

x(k)^ →soluci´o ∀x(0).

• Val. pr. Q = M −^1 N clau converg`encia:
  1. M`etode convergent ⇔ ρ(Q) < 1.
  2. A est. diag. dom. ⇒ J, G-S converg.
  3. A SPD, ω ∈ (0, 2) ⇒ SOR convergeix (G-S tamb´e).
• En molts casos: velocitat SOR > vel. G-S > vel. Jacobi

Introducci´o

Avantatges i inconvenients m`etodes iteratius

• Podem aprofitar dispersitat de la matriu: multiplicar per N i resoldre sistemes amb matriu M costa O(no zeros de A), sense problemes de fill-in i sense emmag. matrius.
• Sovint, mem`oria requerida = O(n)
• Memoria O(n) molt gran si n = O(10^6 ) ⇒ no hi ha memoria per emmagatzemar la matriu, ni amb estructures de dades sofisticades ⇒ m`etodes directes no aplicables.
• Metode iteratiu rapid per matriu dispersa ⇒ reducci´o important cost computacional resp. m`etode directe o alternativa ´unica.
• Inconvenient: Cal garantir convergencia a priori ⇒ limitaci´o de capacitat d’´us als casos anomenats. (pero aquests s´on suficientment importants)

GRADIENT CONJUGAT

Gradient conjugat

Implementaci´o i cost

• Cost:
  • axpy ≡ y = y + ax ⇒ 2 n
  • dot ≡ (x, y) = xT^ y ⇒ 2 n
• Cost aproximat iteraci´o

cost(A × v) + 3 cost(axpy) + 2 cost(dot) ≈ cost(A × v) + 10n

• memoria necessaria implementaci´o eficient: 3 vectors adicionals ri, qi = Api, pi.

Gradient conjugat

Propietats fonamentals

• A SPD ⇒ (x, y)A = (x, Ay) prod. escalar.
• p(x) = pnxn^ + · · · + p 1 x + p 0 ⇒ p(A) = pnAn^ + · · · + p 1 A + p 0 I.
• r 0 ∈ Rn, Ki = {p(A)r 0 : grau(p) ≤ i}.
• Proposici´o.
  • ∃ fi, gi, hi de grau i
    1. xi+1 = x 0 + fi(A)r 0 ∈ x 0 + Ki
    2. pi = gi(A)r 0 ∈ Ki
    3. ri = hi(A)r 0 ∈ Ki, hi(0) = 1
  • (ri, rj) = 0, (pi, pj)A = 0, ∀i 6 = j.
  • Si ri 6 = 0
    1. {r 0 ,... , ri} base ortogonal de Ki respecte (·, ·).
    2. {p 0 ,... , pi} base ortogonal de Ki respecte (·, ·)A.

Gradient conjugat

Converg`encia

• 0 < λ 1 ≤ · · · ≤ λn, valors propis de A
• κ = λ λn 1 condici´o en norma 2.
• ei = x∗ − xi, x∗ soluci´o del sistema Ax = b
• ‖x‖Z =

(x, Zx) (Z = A, A−^1 )

• Proposici´o
  1. ‖ei‖A = ‖ri‖A− 1
  2. ri = hi(A)r 0 ⇒ ‖ri‖A− 1 = ‖hi(A)r 0 ‖A− 1 = min grau(hi)≤i,hi(0)=

‖h(A)r 0 ‖A− 1.

3. ‖ei‖A ≤

√κ−^1 κ+

)i ‖e 0 ‖A

4. ‖ei‖ ≤

(√κ− 1 √ κ+

)i κ‖e 0 ‖

Gradient conjugat

Converg`encia

R =

κ − 1 √ κ + 1

• limi ‖ei‖ = 0 ⇒o xi → x∗ i fita sobre decreixement de ‖ei‖.
• Com

lim κ→ 1

κ − 1 √ κ + 1

= 0 (^) κlim→∞

κ − 1 √ κ + 1

⇒ convergencia m´es rapida quan κ ↘ 1.

• En general

‖ei‖ ‖e 0 ‖

≤ ε, ∀i > log(ε/κ)/ log(R)

• ε = 10−^6 κ = 100 ⇒ R ≈ 0. 8 ⇒ i = log(10−^6 / 102 )/ log(0.8) ≈ 83 iteracions.
• ε = 10−^6 , κ = 10^8 ⇒ R ≈ 0. 9998 ⇒ i = log(10−^6 / 108 )/ log(0.9998) ≈ 1. 6 e5 its!

GRADIENT CONJUGAT

PRECONDICIONAT

Gradient conjugat precondicionat

Precondicionament

• ‖ei‖ ≤

(√κ− 1 √ κ+

)i κ‖e 0 ‖ ⇒ converg`encia controlada per κ.

• κ petit ⇒ m`etode competitiu.
• Malhauradament κ sol ser gran!
• Idea: aconseguir sistema “equivalent” amb millor condici´o.
• M matriu SPD, M ≈ A ⇒ sistema equivalent M −^1 Ax = M −^1 b.
• M −^1 A no ´es SPD ⇒ no podem aplicar CG.

Gradient conjugat precondicionat

Converg`encia

• Valors propis de M −^1 A: λ 1 ≤ · · · ≤ λn
• Suposem M bona aproximaci´o a A (λi ≈ 1, ∀i)
• A ´es semblant a M −^1 A:

C−^1 AC = C−^1 C−T^ AC−^1 C = M −^1 A,

⇒ A i M −^1 A mateixos valors propis ⇒ κ(A) = λ λn 1 ≈ 1 (ja que λ 1 i λn ≈ 1) ⇒ convergencia molt rapida.

Gradient conjugat precondicionat

Gradient conjugat precondicionat (PCG)

• algoritme dels gradients conjugats per a Ax = b, aproximaci´o inicial x 0 = Cx 0.
• Notaci´o:
  1. ri: residual per a xi
  2. pi: vectors auxiliars de CG.

inici: A, b, x 0 r 0 = b − Ax 0 p 0 = r 0 for i = 0,... αi = (ri, ri)/(pi, Api) xi+1 = xi + αipi ri+1 = ri − αiApi βi = (ri+1, ri+1)/(ri, ri) pi+1 = ri+1 + βipi end xt = C−^1 xt

• t ≡ ´ındex parada bucle