forked from TheAlgorithms/Python
- Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathconjugate_gradient.py
More file actions
Latest commit
180 lines (143 loc) · 5.03 KB
/
Copy pathconjugate_gradient.py
File metadata and controls
180 lines (143 loc) · 5.03 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
"""
Resources:
- https://en.wikipedia.org/wiki/Conjugate_gradient_method
- https://en.wikipedia.org/wiki/Definite_symmetric_matrix
"""
fromtypingimportAny
importnumpyasnp
def_is_matrix_spd(matrix: np.ndarray) ->bool:
"""
Returns True if input matrix is symmetric positive definite.
Returns False otherwise.
For a matrix to be SPD, all eigenvalues must be positive.
>>> import numpy as np
>>> matrix = np.array([
... [4.12401784, -5.01453636, -0.63865857],
... [-5.01453636, 12.33347422, -3.40493586],
... [-0.63865857, -3.40493586, 5.78591885]])
>>> _is_matrix_spd(matrix)
True
>>> matrix = np.array([
... [0.34634879, 1.96165514, 2.18277744],
... [0.74074469, -1.19648894, -1.34223498],
... [-0.7687067 , 0.06018373, -1.16315631]])
>>> _is_matrix_spd(matrix)
False
"""
# Ensure matrix is square.
assertnp.shape(matrix)[0] ==np.shape(matrix)[1]
# If matrix not symmetric, exit right away.
ifnp.allclose(matrix, matrix.T) isFalse:
returnFalse
# Get eigenvalues and eignevectors for a symmetric matrix.
eigen_values, _=np.linalg.eigh(matrix)
# Check sign of all eigenvalues.
# np.all returns a value of type np.bool_
returnbool(np.all(eigen_values>0))
def_create_spd_matrix(dimension: int) ->Any:
"""
Returns a symmetric positive definite matrix given a dimension.
Input:
dimension gives the square matrix dimension.
Output:
spd_matrix is an diminesion x dimensions symmetric positive definite (SPD) matrix.
>>> import numpy as np
>>> dimension = 3
>>> spd_matrix = _create_spd_matrix(dimension)
>>> _is_matrix_spd(spd_matrix)
True
"""
rng=np.random.default_rng()
random_matrix=rng.normal(size=(dimension, dimension))
spd_matrix=np.dot(random_matrix, random_matrix.T)
assert_is_matrix_spd(spd_matrix)
returnspd_matrix
defconjugate_gradient(
spd_matrix: np.ndarray,
load_vector: np.ndarray,
max_iterations: int=1000,
tol: float=1e-8,
) ->Any:
"""
Returns solution to the linear system np.dot(spd_matrix, x) = b.
Input:
spd_matrix is an NxN Symmetric Positive Definite (SPD) matrix.
load_vector is an Nx1 vector.
Output:
x is an Nx1 vector that is the solution vector.
>>> import numpy as np
>>> spd_matrix = np.array([
... [8.73256573, -5.02034289, -2.68709226],
... [-5.02034289, 3.78188322, 0.91980451],
... [-2.68709226, 0.91980451, 1.94746467]])
>>> b = np.array([
... [-5.80872761],
... [ 3.23807431],
... [ 1.95381422]])
>>> conjugate_gradient(spd_matrix, b)
array([[-0.63114139],
[-0.01561498],
[ 0.13979294]])
"""
# Ensure proper dimensionality.
assertnp.shape(spd_matrix)[0] ==np.shape(spd_matrix)[1]
assertnp.shape(load_vector)[0] ==np.shape(spd_matrix)[0]
assert_is_matrix_spd(spd_matrix)
# Initialize solution guess, residual, search direction.
x0=np.zeros((np.shape(load_vector)[0], 1))
r0=np.copy(load_vector)
p0=np.copy(r0)
# Set initial errors in solution guess and residual.
error_residual=1e9
error_x_solution=1e9
error=1e9
# Set iteration counter to threshold number of iterations.
iterations=0
whileerror>tol:
# Save this value so we only calculate the matrix-vector product once.
w=np.dot(spd_matrix, p0)
# The main algorithm.
# Update search direction magnitude.
alpha=np.dot(r0.T, r0) /np.dot(p0.T, w)
# Update solution guess.
x=x0+alpha*p0
# Calculate new residual.
r=r0-alpha*w
# Calculate new Krylov subspace scale.
beta=np.dot(r.T, r) /np.dot(r0.T, r0)
# Calculate new A conjuage search direction.
p=r+beta*p0
# Calculate errors.
error_residual=np.linalg.norm(r-r0)
error_x_solution=np.linalg.norm(x-x0)
error=np.maximum(error_residual, error_x_solution)
# Update variables.
x0=np.copy(x)
r0=np.copy(r)
p0=np.copy(p)
# Update number of iterations.
iterations+=1
ifiterations>max_iterations:
break
returnx
deftest_conjugate_gradient() ->None:
"""
>>> test_conjugate_gradient() # self running tests
"""
# Create linear system with SPD matrix and known solution x_true.
dimension=3
spd_matrix=_create_spd_matrix(dimension)
rng=np.random.default_rng()
x_true=rng.normal(size=(dimension, 1))
b=np.dot(spd_matrix, x_true)
# Numpy solution.
x_numpy=np.linalg.solve(spd_matrix, b)
# Our implementation.
x_conjugate_gradient=conjugate_gradient(spd_matrix, b)
# Ensure both solutions are close to x_true (and therefore one another).
assertnp.linalg.norm(x_numpy-x_true) <=1e-6
assertnp.linalg.norm(x_conjugate_gradient-x_true) <=1e-6
if__name__=="__main__":
importdoctest
doctest.testmod()
test_conjugate_gradient()