These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more
1 | <?php |
||
2 | |||
3 | declare(strict_types=1); |
||
4 | |||
5 | namespace Phpml\Helper\Optimizer; |
||
6 | |||
7 | use Closure; |
||
8 | |||
9 | /** |
||
10 | * Conjugate Gradient method to solve a non-linear f(x) with respect to unknown x |
||
11 | * See https://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method) |
||
12 | * |
||
13 | * The method applied below is explained in the below document in a practical manner |
||
14 | * - http://web.cs.iastate.edu/~cs577/handouts/conjugate-gradient.pdf |
||
15 | * |
||
16 | * However it is compliant with the general Conjugate Gradient method with |
||
17 | * Fletcher-Reeves update method. Note that, the f(x) is assumed to be one-dimensional |
||
18 | * and one gradient is utilized for all dimensions in the given data. |
||
19 | */ |
||
20 | class ConjugateGradient extends GD |
||
21 | { |
||
22 | public function runOptimization(array $samples, array $targets, Closure $gradientCb): array |
||
23 | { |
||
24 | $this->samples = $samples; |
||
25 | $this->targets = $targets; |
||
26 | $this->gradientCb = $gradientCb; |
||
27 | $this->sampleCount = count($samples); |
||
28 | $this->costValues = []; |
||
29 | |||
30 | $d = MP::muls($this->gradient($this->theta), -1); |
||
31 | |||
32 | for ($i = 0; $i < $this->maxIterations; ++$i) { |
||
33 | // Obtain α that minimizes f(θ + α.d) |
||
34 | $alpha = $this->getAlpha($d); |
||
35 | |||
36 | // θ(k+1) = θ(k) + α.d |
||
37 | $thetaNew = $this->getNewTheta($alpha, $d); |
||
38 | |||
39 | // β = ||∇f(x(k+1))||² ∕ ||∇f(x(k))||² |
||
40 | $beta = $this->getBeta($thetaNew); |
||
41 | |||
42 | // d(k+1) =–∇f(x(k+1)) + β(k).d(k) |
||
43 | $d = $this->getNewDirection($thetaNew, $beta, $d); |
||
44 | |||
45 | // Save values for the next iteration |
||
46 | $oldTheta = $this->theta; |
||
47 | $this->costValues[] = $this->cost($thetaNew); |
||
48 | |||
49 | $this->theta = $thetaNew; |
||
50 | if ($this->enableEarlyStop && $this->earlyStop($oldTheta)) { |
||
51 | break; |
||
52 | } |
||
53 | } |
||
54 | |||
55 | $this->clear(); |
||
56 | |||
57 | return $this->theta; |
||
58 | } |
||
59 | |||
60 | /** |
||
61 | * Executes the callback function for the problem and returns |
||
62 | * sum of the gradient for all samples & targets. |
||
63 | */ |
||
64 | protected function gradient(array $theta): array |
||
65 | { |
||
66 | [, $updates, $penalty] = parent::gradient($theta); |
||
0 ignored issues
–
show
|
|||
67 | |||
68 | // Calculate gradient for each dimension |
||
69 | $gradient = []; |
||
70 | for ($i = 0; $i <= $this->dimensions; ++$i) { |
||
71 | if ($i === 0) { |
||
72 | $gradient[$i] = array_sum($updates); |
||
73 | } else { |
||
74 | $col = array_column($this->samples, $i - 1); |
||
75 | $error = 0; |
||
76 | foreach ($col as $index => $val) { |
||
77 | $error += $val * $updates[$index]; |
||
78 | } |
||
79 | |||
80 | $gradient[$i] = $error + $penalty * $theta[$i]; |
||
81 | } |
||
82 | } |
||
83 | |||
84 | return $gradient; |
||
85 | } |
||
86 | |||
87 | /** |
||
88 | * Returns the value of f(x) for given solution |
||
89 | */ |
||
90 | protected function cost(array $theta): float |
||
91 | { |
||
92 | [$cost] = parent::gradient($theta); |
||
93 | |||
94 | return array_sum($cost) / $this->sampleCount; |
||
95 | } |
||
96 | |||
97 | /** |
||
98 | * Calculates alpha that minimizes the function f(θ + α.d) |
||
99 | * by performing a line search that does not rely upon the derivation. |
||
100 | * |
||
101 | * There are several alternatives for this function. For now, we |
||
102 | * prefer a method inspired from the bisection method for its simplicity. |
||
103 | * This algorithm attempts to find an optimum alpha value between 0.0001 and 0.01 |
||
104 | * |
||
105 | * Algorithm as follows: |
||
106 | * a) Probe a small alpha (0.0001) and calculate cost function |
||
107 | * b) Probe a larger alpha (0.01) and calculate cost function |
||
108 | * b-1) If cost function decreases, continue enlarging alpha |
||
109 | * b-2) If cost function increases, take the midpoint and try again |
||
110 | */ |
||
111 | protected function getAlpha(array $d): float |
||
112 | { |
||
113 | $small = MP::muls($d, 0.0001); |
||
114 | $large = MP::muls($d, 0.01); |
||
115 | |||
116 | // Obtain θ + α.d for two initial values, x0 and x1 |
||
117 | $x0 = MP::add($this->theta, $small); |
||
118 | $x1 = MP::add($this->theta, $large); |
||
119 | |||
120 | $epsilon = 0.0001; |
||
121 | $iteration = 0; |
||
122 | do { |
||
123 | $fx1 = $this->cost($x1); |
||
124 | $fx0 = $this->cost($x0); |
||
125 | |||
126 | // If the difference between two values is small enough |
||
127 | // then break the loop |
||
128 | if (abs($fx1 - $fx0) <= $epsilon) { |
||
129 | break; |
||
130 | } |
||
131 | |||
132 | if ($fx1 < $fx0) { |
||
133 | $x0 = $x1; |
||
134 | $x1 = MP::adds($x1, 0.01); // Enlarge second |
||
135 | } else { |
||
136 | $x1 = MP::divs(MP::add($x1, $x0), 2.0); |
||
137 | } // Get to the midpoint |
||
138 | |||
139 | $error = $fx1 / $this->dimensions; |
||
140 | } while ($error <= $epsilon || $iteration++ < 10); |
||
141 | |||
142 | // Return α = θ / d |
||
143 | // For accuracy, choose a dimension which maximize |d[i]| |
||
144 | $imax = 0; |
||
145 | for ($i = 1; $i <= $this->dimensions; ++$i) { |
||
146 | if (abs($d[$i]) > abs($d[$imax])) { |
||
147 | $imax = $i; |
||
148 | } |
||
149 | } |
||
150 | |||
151 | if ($d[$imax] == 0) { |
||
152 | return $x1[$imax] - $this->theta[$imax]; |
||
153 | } |
||
154 | |||
155 | return ($x1[$imax] - $this->theta[$imax]) / $d[$imax]; |
||
156 | } |
||
157 | |||
158 | /** |
||
159 | * Calculates new set of solutions with given alpha (for each θ(k)) and |
||
160 | * gradient direction. |
||
161 | * |
||
162 | * θ(k+1) = θ(k) + α.d |
||
163 | */ |
||
164 | protected function getNewTheta(float $alpha, array $d): array |
||
165 | { |
||
166 | return MP::add($this->theta, MP::muls($d, $alpha)); |
||
167 | } |
||
168 | |||
169 | /** |
||
170 | * Calculates new beta (β) for given set of solutions by using |
||
171 | * Fletcher–Reeves method. |
||
172 | * |
||
173 | * β = ||f(x(k+1))||² ∕ ||f(x(k))||² |
||
174 | * |
||
175 | * See: |
||
176 | * R. Fletcher and C. M. Reeves, "Function minimization by conjugate gradients", Comput. J. 7 (1964), 149–154. |
||
177 | */ |
||
178 | protected function getBeta(array $newTheta): float |
||
179 | { |
||
180 | $gNew = $this->gradient($newTheta); |
||
181 | $gOld = $this->gradient($this->theta); |
||
182 | $dNew = 0; |
||
183 | $dOld = 1e-100; |
||
184 | for ($i = 0; $i <= $this->dimensions; ++$i) { |
||
185 | $dNew += $gNew[$i] ** 2; |
||
186 | $dOld += $gOld[$i] ** 2; |
||
187 | } |
||
188 | |||
189 | return $dNew / $dOld; |
||
190 | } |
||
191 | |||
192 | /** |
||
193 | * Calculates the new conjugate direction |
||
194 | * |
||
195 | * d(k+1) =–∇f(x(k+1)) + β(k).d(k) |
||
196 | */ |
||
197 | protected function getNewDirection(array $theta, float $beta, array $d): array |
||
198 | { |
||
199 | $grad = $this->gradient($theta); |
||
200 | |||
201 | return MP::add(MP::muls($grad, -1), MP::muls($d, $beta)); |
||
202 | } |
||
203 | } |
||
204 | |||
205 | /** |
||
206 | * Handles element-wise vector operations between vector-vector |
||
207 | * and vector-scalar variables |
||
208 | */ |
||
209 | class MP |
||
210 | { |
||
211 | /** |
||
212 | * Element-wise <b>multiplication</b> of two vectors of the same size |
||
213 | */ |
||
214 | public static function mul(array $m1, array $m2): array |
||
215 | { |
||
216 | $res = []; |
||
217 | foreach ($m1 as $i => $val) { |
||
218 | $res[] = $val * $m2[$i]; |
||
219 | } |
||
220 | |||
221 | return $res; |
||
222 | } |
||
223 | |||
224 | /** |
||
225 | * Element-wise <b>division</b> of two vectors of the same size |
||
226 | */ |
||
227 | public static function div(array $m1, array $m2): array |
||
228 | { |
||
229 | $res = []; |
||
230 | foreach ($m1 as $i => $val) { |
||
231 | $res[] = $val / $m2[$i]; |
||
232 | } |
||
233 | |||
234 | return $res; |
||
235 | } |
||
236 | |||
237 | /** |
||
238 | * Element-wise <b>addition</b> of two vectors of the same size |
||
239 | */ |
||
240 | public static function add(array $m1, array $m2, int $mag = 1): array |
||
241 | { |
||
242 | $res = []; |
||
243 | foreach ($m1 as $i => $val) { |
||
244 | $res[] = $val + $mag * $m2[$i]; |
||
245 | } |
||
246 | |||
247 | return $res; |
||
248 | } |
||
249 | |||
250 | /** |
||
251 | * Element-wise <b>subtraction</b> of two vectors of the same size |
||
252 | */ |
||
253 | public static function sub(array $m1, array $m2): array |
||
254 | { |
||
255 | return self::add($m1, $m2, -1); |
||
256 | } |
||
257 | |||
258 | /** |
||
259 | * Element-wise <b>multiplication</b> of a vector with a scalar |
||
260 | */ |
||
261 | public static function muls(array $m1, float $m2): array |
||
262 | { |
||
263 | $res = []; |
||
264 | foreach ($m1 as $val) { |
||
265 | $res[] = $val * $m2; |
||
266 | } |
||
267 | |||
268 | return $res; |
||
269 | } |
||
270 | |||
271 | /** |
||
272 | * Element-wise <b>division</b> of a vector with a scalar |
||
273 | */ |
||
274 | public static function divs(array $m1, float $m2): array |
||
275 | { |
||
276 | $res = []; |
||
277 | foreach ($m1 as $val) { |
||
278 | $res[] = $val / ($m2 + 1e-32); |
||
279 | } |
||
280 | |||
281 | return $res; |
||
282 | } |
||
283 | |||
284 | /** |
||
285 | * Element-wise <b>addition</b> of a vector with a scalar |
||
286 | */ |
||
287 | public static function adds(array $m1, float $m2, int $mag = 1): array |
||
288 | { |
||
289 | $res = []; |
||
290 | foreach ($m1 as $val) { |
||
291 | $res[] = $val + $mag * $m2; |
||
292 | } |
||
293 | |||
294 | return $res; |
||
295 | } |
||
296 | |||
297 | /** |
||
298 | * Element-wise <b>subtraction</b> of a vector with a scalar |
||
299 | */ |
||
300 | public static function subs(array $m1, float $m2): array |
||
301 | { |
||
302 | return self::adds($m1, $m2, -1); |
||
303 | } |
||
304 | } |
||
305 |
This check marks access to variables or properties that have not been declared yet. While PHP has no explicit notion of declaring a variable, accessing it before a value is assigned to it is most likely a bug.