Passed
Push — master ( e83f7b...d953ef )
by Arkadiusz
03:28
created

src/Phpml/Helper/Optimizer/ConjugateGradient.php (1 issue)

Upgrade to new PHP Analysis Engine

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(array_sum($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
        [, $gradient] = parent::gradient($theta);
0 ignored issues
show
The variable $gradient does not exist. Did you forget to declare it?

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.

Loading history...
67
68
        return $gradient;
69
    }
70
71
    /**
72
     * Returns the value of f(x) for given solution
73
     */
74
    protected function cost(array $theta): float
75
    {
76
        [$cost] = parent::gradient($theta);
77
78
        return array_sum($cost) / $this->sampleCount;
79
    }
80
81
    /**
82
     * Calculates alpha that minimizes the function f(θ + α.d)
83
     * by performing a line search that does not rely upon the derivation.
84
     *
85
     * There are several alternatives for this function. For now, we
86
     * prefer a method inspired from the bisection method for its simplicity.
87
     * This algorithm attempts to find an optimum alpha value between 0.0001 and 0.01
88
     *
89
     * Algorithm as follows:
90
     *  a) Probe a small alpha  (0.0001) and calculate cost function
91
     *  b) Probe a larger alpha (0.01) and calculate cost function
92
     *		b-1) If cost function decreases, continue enlarging alpha
93
     *		b-2) If cost function increases, take the midpoint and try again
94
     */
95
    protected function getAlpha(float $d): float
96
    {
97
        $small = 0.0001 * $d;
98
        $large = 0.01 * $d;
99
100
        // Obtain θ + α.d for two initial values, x0 and x1
101
        $x0 = MP::adds($this->theta, $small);
102
        $x1 = MP::adds($this->theta, $large);
103
104
        $epsilon = 0.0001;
105
        $iteration = 0;
106
        do {
107
            $fx1 = $this->cost($x1);
108
            $fx0 = $this->cost($x0);
109
110
            // If the difference between two values is small enough
111
            // then break the loop
112
            if (abs($fx1 - $fx0) <= $epsilon) {
113
                break;
114
            }
115
116
            if ($fx1 < $fx0) {
117
                $x0 = $x1;
118
                $x1 = MP::adds($x1, 0.01); // Enlarge second
119
            } else {
120
                $x1 = MP::divs(MP::add($x1, $x0), 2.0);
121
            } // Get to the midpoint
122
123
            $error = $fx1 / $this->dimensions;
124
        } while ($error <= $epsilon || $iteration++ < 10);
125
126
        //  Return α = θ / d
127
        if ($d == 0) {
128
            return $x1[0] - $this->theta[0];
129
        }
130
131
        return ($x1[0] - $this->theta[0]) / $d;
132
    }
133
134
    /**
135
     * Calculates new set of solutions with given alpha (for each θ(k)) and
136
     * gradient direction.
137
     *
138
     * θ(k+1) = θ(k) + α.d
139
     */
140
    protected function getNewTheta(float $alpha, array $d): array
141
    {
142
        $theta = $this->theta;
143
144
        for ($i = 0; $i < $this->dimensions + 1; ++$i) {
145
            if ($i === 0) {
146
                $theta[$i] += $alpha * array_sum($d);
147
            } else {
148
                $sum = 0.0;
149
                foreach ($this->samples as $si => $sample) {
150
                    $sum += $sample[$i - 1] * $d[$si] * $alpha;
151
                }
152
153
                $theta[$i] += $sum;
154
            }
155
        }
156
157
        return $theta;
158
    }
159
160
    /**
161
     * Calculates new beta (β) for given set of solutions by using
162
     * Fletcher–Reeves method.
163
     *
164
     * β = ||f(x(k+1))||²  ∕  ||f(x(k))||²
165
     *
166
     * See:
167
     *  R. Fletcher and C. M. Reeves, "Function minimization by conjugate gradients", Comput. J. 7 (1964), 149–154.
168
     */
169
    protected function getBeta(array $newTheta): float
170
    {
171
        $dNew = array_sum($this->gradient($newTheta));
172
        $dOld = array_sum($this->gradient($this->theta)) + 1e-100;
173
174
        return $dNew ** 2 / $dOld ** 2;
175
    }
176
177
    /**
178
     * Calculates the new conjugate direction
179
     *
180
     * d(k+1) =–∇f(x(k+1)) + β(k).d(k)
181
     */
182
    protected function getNewDirection(array $theta, float $beta, array $d): array
183
    {
184
        $grad = $this->gradient($theta);
185
186
        return MP::add(MP::muls($grad, -1), MP::muls($d, $beta));
187
    }
188
}
189
190
/**
191
 * Handles element-wise vector operations between vector-vector
192
 * and vector-scalar variables
193
 */
194
class MP
195
{
196
    /**
197
     * Element-wise <b>multiplication</b> of two vectors of the same size
198
     */
199
    public static function mul(array $m1, array $m2): array
200
    {
201
        $res = [];
202
        foreach ($m1 as $i => $val) {
203
            $res[] = $val * $m2[$i];
204
        }
205
206
        return $res;
207
    }
208
209
    /**
210
     * Element-wise <b>division</b> of two vectors of the same size
211
     */
212
    public static function div(array $m1, array $m2): array
213
    {
214
        $res = [];
215
        foreach ($m1 as $i => $val) {
216
            $res[] = $val / $m2[$i];
217
        }
218
219
        return $res;
220
    }
221
222
    /**
223
     * Element-wise <b>addition</b> of two vectors of the same size
224
     */
225
    public static function add(array $m1, array $m2, int $mag = 1): array
226
    {
227
        $res = [];
228
        foreach ($m1 as $i => $val) {
229
            $res[] = $val + $mag * $m2[$i];
230
        }
231
232
        return $res;
233
    }
234
235
    /**
236
     * Element-wise <b>subtraction</b> of two vectors of the same size
237
     */
238
    public static function sub(array $m1, array $m2): array
239
    {
240
        return self::add($m1, $m2, -1);
241
    }
242
243
    /**
244
     * Element-wise <b>multiplication</b> of a vector with a scalar
245
     */
246
    public static function muls(array $m1, float $m2): array
247
    {
248
        $res = [];
249
        foreach ($m1 as $val) {
250
            $res[] = $val * $m2;
251
        }
252
253
        return $res;
254
    }
255
256
    /**
257
     * Element-wise <b>division</b> of a vector with a scalar
258
     */
259
    public static function divs(array $m1, float $m2): array
260
    {
261
        $res = [];
262
        foreach ($m1 as $val) {
263
            $res[] = $val / ($m2 + 1e-32);
264
        }
265
266
        return $res;
267
    }
268
269
    /**
270
     * Element-wise <b>addition</b> of a vector with a scalar
271
     */
272
    public static function adds(array $m1, float $m2, int $mag = 1): array
273
    {
274
        $res = [];
275
        foreach ($m1 as $val) {
276
            $res[] = $val + $mag * $m2;
277
        }
278
279
        return $res;
280
    }
281
282
    /**
283
     * Element-wise <b>subtraction</b> of a vector with a scalar
284
     */
285
    public static function subs(array $m1, float $m2): array
286
    {
287
        return self::adds($m1, $m2, -1);
288
    }
289
}
290