Gradient method

Tsung-Ming Huang


Matrix Computation, 2016, NTNU

1

Plan

Gradient method

Conjugate gradient method Preconditioner

2

Gradient method

3

Theorem

4

Ax = b, A : s.p.d

Inner product

< x, y > = x

y for any x, y ∈!

g(x) = < x, Ax > −2 < x,b > = x

Ax − 2x

b

Define

Theorem

A : s.p.d

x

is the sol. of Ax = b g(x

) = min

g(x)

A^⊤ = A

x^⊤Ax > 0, ∀ x ≠ 0

A : symmetric positive definite if

Definition

Proof

5

g(x) = < x, Ax > −2 < x,b >

Assume

x

is the sol. of Ax = b

= < x − x

, A(x − x

) > + < x, Ax

> + < x

, Ax >

− < x

, Ax

> −2 < x,b >

= < x − x

, A(x − x

) > − < x

, Ax

>

+ 2 < x, Ax

> −2 < x,b >

= < x − x

, A(x − x

) > − < x

, Ax

> +2 < x, Ax

− b >

Ax

= b

= < x − x

, A(x − x

) > − < x

, Ax

>

g(x

) = min

g(x)

< x − x

, A(x − x

) > ≥ 0

Proof

6

Assume

g(x

) = min

g(x)

Fixed vectors and , for any

α

f ( α ⁾ ≡ g(x + α ^v) = < x + α ^{v, Ax} + α ^Av > −2 < x + α ^v,b >

= < x, Ax > + α < v, Ax > + α < x, Av > + α

< v, Av >

−2 < x,b > −2 α < v,b >

= < x, Ax > −2 < x,b > +2 α < v, Ax >

−2 α < v,b > + α

< v, Av >

= g(x) + 2 α < v, Ax − b > + α

< v, Av >

g(x + ˆ α ^v) = f ( ˆ α ⁾ = g(x) − 2 < v,b − Ax >

< v, Av > < v,b − Ax >

+ < v,b − Ax >

< v, Av >

⎛ ⎝⎜ ⎞

⎠⎟

2

< v, Av >

Proof

7

f ( α ⁾ = g(x) + 2 α < v, Ax − b > + α

< v, Av >

is a quadratic function of

f

α

A : s.p.d

has a minimal value when

α

f ′( ˆ α ⁾ = 2 < v, Ax − b > +2 ˆ α < v, Av > = 0 α ˆ = − < v, Ax − b >

< v, Av > = < v,b − Ax >

< v, Av >

= g(x) − < v,b − Ax >

< v, Av >

Proof

8

g(x + ˆ α ^v) = g(x) − < v,b − Ax >

< v, Av >

g(x + ˆ α ^v) = g(x) if < v,b − Ax > = 0 g(x + ˆ α ^v) < g(x) if < v,b − Ax > ≠ 0

∀ v ≠ 0

Suppose that

g(x

) = min

g(x)

g(x

+ ˆ α ^v) ≥ g(x

) for any v

< v,b − Ax

> = 0, ∀ v Ax

= b

9

α = < v,b − Ax >

< v, Av > = < v,r >

< v, Av > , r ≡ b − Ax

If r ≠ 0 and < v,r > ≠ 0

g(x +

α

< v, Av > < g(x)

is closer to than is

x + α ^v x^{∗ x}

Given and

k = 1,2,3,!

For

α _k = < v^{(k )},b − Ax^(k−1) >

< v^{(k )}, Av^{(k )} > , x^{(k )} = x^(k−1) + α_k^{v(k )}

Choose a new search direction

Steepest descent

10

Question

How to choose s.t. rapidly? { }

{ }

Let be a differentiable function on

x

Φ(x + ε ^p) − Φ(x)

ε ^{= ∇Φ(x)}

⊤

p +O( ε ⁾

p = − ∇Φ (x)

‖ ∇Φ(x)‖ (i.e., the largest descent) The right hand side takes minimum at

for all with (neglect ) p ‖p‖ = 1 ^O( ε ⁾

Steepest descent direction of g

11

Denote x = [x

, x

, !, x

]

g(x) = < x, Ax > −2 < x,b > = a

j=1

∑

∑

^x

^x

^{− 2 x}

∑

^b

∂g

∂x

(x) = 2 a

i=1

∑

^x

^{− 2b}

= 2 A(k,:)x − b (

)

∇g(x) = ∂ g

∂x

(x), ∂g

∂x

, !, ∂g

∂x

(x)

⎡

⎣ ⎢ ⎤

⎦ ⎥

⊤

= 2(Ax − b) = −2r

Steepest descent method (gradient method)

12

Given

r_k−1 = b − Ax^(k−1)

α

< r_k−1, Ar_k−1 >

x^{(k )} = x^(k−1) +

α

Else

If , then

Stop;

End

k = 1,2,3,!

For

End

Convergence Theorem

λ₁ ≥ λ₂ ≥! ≥ λ_n > 0 : eigenvalues

x^{(k )}, x^(k−1) : approx. sol.

x^∗ : exact sol.

where ‖x‖_A = x^⊤Ax

‖x^{(k )} − x^*‖_A

≤ λ₁ − λ_n λ₁ + λ_n

⎛

⎝⎜

⎞

⎠⎟‖x^(k−1) − x^*‖_A

Conjugate gradient method

13

A-orthogonal

14

If κ ^(A) = λ₁

λ_n ^{is large}

λ

− λ

λ

+ λ

^{≈ 1}

Convergence is very slow

Improvement

Choose A-orthogonal search directions

Definition

p,q ∈!ⁿ

are called A-orthogonal (A-conjugate) if

p^⊤Aq = 0

NOT recommend it

Lemma

15

v

, …,v

≠ 0 : pairwisely A-conjugate

Proof

0 = c

j=1

∑

^v

0 = (v

)

A c

j=1

∑

^v

⎛

⎝⎜

⎞

⎠⎟ = c

j=1

∑

^(v

⁾

^Av

^{= c}

^(v

⁾

c

= 0, k = 1,…,n

v

, …,v

: linearly independent

v

, …,v

: linearly independent

Theorem

16

A : symmetric positive definite

v

, …,v

≠ 0 ∈!

: pairwisely A-conjugate x

: given

For , let k = 1,…,n

α

< v_k, Av_k >

x_k = x_k−1 + α _k^v_k

Then

Ax_n = b

< b − Ax_k ,v_j > = 0, for j = 1,2,…,k

Proof

17

x_k = x_k−1 + α_k^v_k Ax_n = Ax_n−1 + α_n^Av_n

= ! = Ax₀ + α₁^Av₁ + α₂^Av₂ +!+ α_n^Av_n

= (Ax_n−2 + α_n−1^Av_n−1⁾ + α_n^Av_n

< Ax_n − b,v_k >

= < Ax₀ − b,v_k > +α₁ < Av₁,v_k > +!+ α_n < Av_n,v_k >

= < Ax₀ − b,v_k > +α₁ < v₁, Av_k > +!+ α_n < v_n, Av_k >

= < Ax₀ − b,v_k > +α _k < v_k , Av_k >

= < Ax₀ − b,v_k > + < v_k,b − Ax_k−1 >

< v_k , Av_k > < v_k, Av_k >

= < Ax₀ − b,v_k > + < v_k ,b − Ax_k−1 >

Proof

18

< Ax_n − b,v_k > = < Ax₀ − b,v_k > + < v_k ,b − Ax_k−1 >

= < Ax₀ − b,v_k >

+ < v_k ,b − Ax₀ + Ax₀ − Ax₁ +!− Ax_k−2 + Ax_k−2 − Ax_k−1 >

= < Ax₀ − b,v_k > + < v_k ,b − Ax₀ >

+ < v_k , Ax₀ − Ax₁ > +!+ < v_k, Ax_k−2 − Ax_k−1 >

= < v_k, Ax₀ − Ax₁ > +!+ < v_k , Ax_k−2 − Ax_k−1 >

x_i = x_i−1 + α_i^v_i^, ∀i Ax_i = Ax_i−1 + α_i^Av_i Ax_i−1 − Ax_i = −α_i^Av_i

< Ax_n − b,v_k > = −α₁ < v_k , Av₁ > −!− α_k−1 < v_k , Av_k−1 > = 0 Ax_n = b

Proof

19

< b − Ax_k ,v_j > = 0, for j = 1,2,…,k

< r_k−1,v_j > = 0, for j = 1,2,…,k −1

Assume

r_k = b − Ax_k = b − A(x_k−1 + α_k^v_k ⁾ = r_k−1 − α _k ^Av_k

< r_k ,v_k > = < r_k−1,v_k > −α_k < Av_k ,v_k >

= < r_k−1,v_k > − < v_k ,b − Ax_k−1 >

< v_k , Av_k > < Av_k,v_k > = 0

For

< r_k ,v_j > = < r_k−1,v_j > −α _k < Av_k ,v_j > = 0

Assumption A-conjugate

which is completed the proof by the mathematic induction.

Method of conjugate directions

20

r_k = r_k−1 −

α

α

< v_k, Av_k > , x^{(k )} = x^(k−1) +

α

k = 1,…,n

For

End

Given : pairwisely A-orthogonal

r₀ = b − Ax⁽⁰⁾

Question

How to find A-orthogonal search directions?

A-orthogonalization

21

v₁ v₂

α^v₁

!v₂

!v

= v

− α ^v

⊥ v

0 = v

!v

= v

v

− α ^v

^v

α = v

v

v

v

A-orthogonal

!v

= v

− α ^v

⊥

v

0 = v

A !v

= v

Av

− α ^v

^Av

α = v

Av

v

Av

A-orthogonalization

22

!v

= v

− v

Av

v

Av

v

⊥

v

v

, v

{ } { ^v

^{, !v}

} : A-orthogonal v

, v

, v

{ } { ^v

^{, !v}

^{, !v}

} : A-orthogonal

!v

= v

− α

^v

− α

^!v

⊥

{ v

, !v

}

0 = v

A !v

= v

Av

− α

^v

^Av

α

= v

Av

/ v

Av

0 = !v

A !v

= !v

Av

− α

^!v

^{A !v}

α

= !v

Av

/ !v

A !v

Practical Implementation

23

Given

α

< v₁, Av₁ > , x⁽¹⁾ = x⁽⁰⁾ +

α

r₁ = r₀ − α₁^Av₁ steepest descent direction

v₁, r₁

{ }

v₂ = r₁ + β₁^v₁^, β₁ ^{= − < v1, Ar1 >}

< v₁, Av₁ >

Construct A-orthogonal vector

α

< v₂, Av₂ > , x⁽²⁾ = x⁽¹⁾ +

α

Construct A-orthogonal vector

24

v₁, v₂,r₂

{ }

r₁ = r₀ −

α

v₂^⊤r₂ = v₂^⊤r₁ −

α

v₂^⊤Av₂ v₂^⊤Av₂ = 0

v₁^⊤Ar₂ = r₂^⊤Av₁ =

α

(

)

0 = v₂^⊤r₂ = r

(

β

)

^β

v₃ = r₂ +

β

β

β

v₁^⊤Av₁ ,

β

= r₁^⊤r₂ + β₁^v₁^⊤

(

)

= r₁^⊤r₂ +

β

(

α

)

^β

1, Av₁ > v₁^⊤Av₁

⎛

⎝⎜

⎞

= r₁^⊤r₂ ⎠⎟

25

r₁ = r₀ − α₁^Av₁^, α₁ = < v₁,r₀ >

< v₁, Av₁ >

< v₁,r₁ > = < v₁,r₀ > −α₁ < v₁, Av₁ > = 0

< r₂,r₀ > = < r₂,v₁ > = < r₁,v₁ > −α₂ < Av₂,v₁ > = 0 v₁^⊤Ar₂ =

α

(

)

β₂₁ = − v₁^⊤Ar₂

v₁^⊤Av₁ = 0

v₃ = r₂ + β₂^v₂^, β₂ = − v₂^⊤Ar₂ v₂^⊤Av₂

In general case

26

v_k = r_k−1 + β_k−1^v_k−1 ^{if r}_k−1 ≠ 0

(i). r

{

}

0 = < v_k−1, Av_k > = < v_k−1, Ar_k−1 +

β

= < v_k−1, Ar_k−1 > +β_k−1 < v_k−1, Av_k−1 >

β_k−1 = − < v_k−1, Ar_k−1 >

< v_k−1, Av_k−1 >

Theorem

(ii). v

{

}

Reformula

27

α_k = < v_k,r_k−1 >

< v_k , Av_k > = < r_k−1 + β_k−1^v_k−1^,r_k−1 >

< v_k , Av_k >

= < r_k−1,r_k−1 >

< v_k , Av_k > +

β

< v_k, Av_k > = < r_k−1,r_k−1 >

< v_k, Av_k >

< r_k−1,r_k−1 > =

α

α

^, β

v_k = r_k−1 + β_k−1^v_k−1

r_k = r_k−1 −

α

< r_k ,r_k > = < r_k−1,r_k > −

α

α

= < r_k ,r_k >

< r_k−1,r_k−1 >

β_k = − < v_k , Ar_k >

< v_k , Av_k > = − < r_k , Av_k >

< v_k , Av_k >

Algorithm (Conjugate Gradient Method)

28

r_k+1 = r_k − α _k ^Av_k α _k = < r_k ,r_k >

< v_k, Av_k > , x^(k+1) = x^{(k )} + α_k^v_k

k = 0,1,…

For

End

Given compute

If , then

Stop;

End

β_k = < r_k+1,r_k+1 >

< r_k,r_k > , v_k+1 = r_k+1 + β_k^v_k

Else

well-conditioned

‖r _n‖< tol Ax_n = b

Theorem

ill-conditioned

‖r_k‖< tol k > n

Conjugate Gradient Method

29

Convergence Theorem

λ₁ ≥ λ₂ ≥! ≥ λ_n > 0 : eigenvalues

x^{(k )}

{ }

‖x_G^{(k )} − x^*‖_A ≤

λ

λ

λ

λ

⎛

⎝⎜

⎞

⎠⎟

k

‖x_G⁽⁰⁾ − x^*‖_A =

κ

κ

⎛⎝⎜ ⎞

⎠⎟

k

‖x_G⁽⁰⁾ − x^*‖_A

‖x^{(k )} − x^*‖_A ≤ 2 κ −1 κ +1

⎛

⎝⎜

⎞

⎠⎟

k

‖x₀ − x^*‖_A, κ = λ₁ λ_n x_G^{(k )}

{ }

CG is much better than Gradient method

κ −1

κ +1 > κ −1 κ +1

Preconditioner

30

C⁻¹A x = C⁻¹b

31

Ax = b C^−⊤C^⊤

Choose such thatC κ ^{(C−1AC−⊤ )} < κ ^(A)

Goal

!A !x !b

!A!x = !b

Apply CG method to ^Get !x ^Solve x = C^−⊤ !x

Nothing NEW

Apply CG method to !A!x = !b ^Get x

Question

Algorithm (Conjugate Gradient Method)

32

!r_k+1 = !b − !A!x^(k+1)

α

< !v_k , !A !v_k >

!x^(k+1) = !x^{(k )} + !

α

k = 0,1,…

For

End

Given compute

If , then Stop

β

< !r_k , !r_k >

!v_k+1 = !r_k+1 + !

β

!r_k+1 = C⁻¹b − C

(

)

= C⁻¹(b − Ax_k+1) = C⁻¹r_k+1

= C⁻¹r_k+1

Let

!v_k = C^⊤v_k, w_k = C⁻¹r_k

= < w_k+1,w_k+1 >

< w_k ,w_k >

= < w_k+1,w_k+1 >

< w_k ,w_k >

β!_k = < C⁻¹r_k+1,C⁻¹r_k+1 >

< C⁻¹r_k ,C⁻¹r_k >

Algorithm (Conjugate Gradient Method)

33

α

< !v_k , !A !v_k >

!x^(k+1) = !x^{(k )} + !

α

k = 0,1,…

For

End

Given compute

If , then Stop

< C^⊤v_k ,C⁻¹Av_k >

= v_k^⊤CC⁻¹Av_k = v_k^⊤Av_k

!r_k+1 = C⁻¹r_k+1

β

< w_k ,w_k >

!v_k+1 = !r_k+1 + !

β

α! _k = < C⁻¹r_k ,C⁻¹r_k >

< C^⊤v_k ,C⁻¹AC^−⊤C^⊤v_k >

= < w_k,w_k >

< C^⊤v_k,C⁻¹Av_k >

α! _k ^{= < wk,wk >}

< v_k , Av_k >

= < w_k ,w_k >

< v_k , Av_k >

Algorithm (Conjugate Gradient Method)

34

α

< v_k , Av_k >

!x^(k+1) = !x^{(k )} + !

α

k = 0,1,…

For

End

Given compute

If , then Stop

C⁻¹r_k+1 = C⁻¹r_k − !α _k^{C−1AC−⊤C⊤v}_k

β

< w_k ,w_k >

!v_k+1 = !r_k+1 + !

β

C^⊤x^(k+1) = C^⊤x^{(k )} + !α _k^C⊤v_k x^(k+1) = x^{(k )} + !α _k^v_k

!r_k+1 = !r_k − !

α

C^⊤v_k+1 = C ⁻¹r_k+1 + !β_k^C⊤v_k

!v_k = C^⊤v_k, w_k = C⁻¹r_k

= C^−⊤w_k+1 + !β_k^v_k

v_k+1 = C^−⊤C⁻¹r_k+1 + !β_k^v_k

Algorithm (Conjugate Gradient Method)

35

α

< v_k , Av_k >

x^(k+1) = x^{(k )} + !

α

k = 0,1,…

For

End

Given compute

If , then Stop

β

< w_k ,w_k >

v_k+1 = C^−⊤w_k+1 + !

β

r_k+1 = r_k − !

α

−⊤w_k+1 + !β_k^v_k

need w₀

w₀ = C⁻¹r₀ = C⁻¹(b − Ax⁽⁰⁾) w_k = C⁻¹r_k

need v₀

v₀ = C^−⊤w₀

Solve C w_k+1 = r_k+1

Algorithm (CG Method with preconditioner C)

36

α

x^(k+1) = x^{(k )} +

α

k = 0,1,…

For

End

If , then Stop

β

v_k+1 = z_k+1 +

β

α

Given and compute

Solve and

Solve and

r_k+1 = CC^⊤z_k+1 ≡ Mz_k+1 β_k = < C⁻¹r_k+1,C⁻¹r_k+1 >

< C⁻¹r_k ,C⁻¹r_k >

= < z_k+1,r_k+1 >

< z_k ,r_k >

α _k = < C⁻¹r_k ,C⁻¹r_k >

< C^⊤v_k ,C⁻¹Av_k >

= < z_k ,r_k >

< v_k , Av_k >

Algorithm (CG Method with preconditioner M)

37

k = 0,1,…

For

End

If , then Stop

Given and compute

Solve and set

Solve

α _k = < z_k ,r_k > / < v_k , Av_k >

Compute Compute

β_k = < z_k+1,r_k+1 > / < z_k ,r_k >

Compute

Compute

v_k+1 = z_k+1 +

β

Compute

Choices of M (Criterion)

cond is nearly by 1, i.e.,

38

(M ^−1/2AM ^−1/2 )

M ^−1/2AM ^−1/2 ≈ I, A ≈ M

The linear system must be easily solved. e.g.

is symmetric positive definite

Preconditioner M

Jacobi method

39

A = D + (L +U), M = D

x_k+1 = −D⁻¹(L +U)x_k + D⁻¹b

= −D⁻¹(A − D)x_k + D⁻¹b

= x_k + D⁻¹r_k

Gauss-Seidel

= (D + L)⁻¹(D + L − A)x_k + (D + L)⁻¹b

= x_k + (D + L)⁻¹r_k