Issues in EigenvalueDecomposition.php - New Issues - Inspection of "Linear Discrimant Analysis (LDA) (#82)" - php-ai/php-ml - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( 12b8b1...5b373f )

by Arkadiusz

created 2017-04-25 06:58 UTC

Math/LinearAlgebra/EigenvalueDecomposition.php (1 issue)

Labels

Duplication 1

Severity

Informational 1

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		*
6		* Class to obtain eigenvalues and eigenvectors of a real matrix.
7		*
8		* If A is symmetric, then A = VDV' where the eigenvalue matrix D
9		* is diagonal and the eigenvector matrix V is orthogonal (i.e.
10		* A = V.times(D.times(V.transpose())) and V.times(V.transpose())
11		* equals the identity matrix).
12		*
13		* If A is not symmetric, then the eigenvalue matrix D is block diagonal
14		* with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
15		* lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The
16		* columns of V represent the eigenvectors in the sense that AV = VD,
17		* i.e. A.times(V) equals V.times(D). The matrix V may be badly
18		* conditioned, or even singular, so the validity of the equation
19		* A = VDinverse(V) depends upon V.cond().
20		*
21		* @author Paul Meagher
22		* @license PHP v3.0
23		* @version 1.1
24		*
25		* Slightly changed to adapt the original code to PHP-ML library
26		* @date 2017/04/11
27		* @author Mustafa Karabulut
28		*/
29
30		namespace Phpml\Math\LinearAlgebra;
31
32		use Phpml\Math\Matrix;
33
34		class EigenvalueDecomposition
35		{
36
37		/**
38		* Row and column dimension (square matrix).
39		* @var int
40		*/
41		private $n;
42
43		/**
44		* Internal symmetry flag.
45		* @var int
46		*/
47		private $issymmetric;
48
49		/**
50		* Arrays for internal storage of eigenvalues.
51		* @var array
52		*/
53		private $d = [];
54		private $e = [];
55
56		/**
57		* Array for internal storage of eigenvectors.
58		* @var array
59		*/
60		private $V = [];
61
62		/**
63		* Array for internal storage of nonsymmetric Hessenberg form.
64		* @var array
65		*/
66		private $H = [];
67
68		/**
69		* Working storage for nonsymmetric algorithm.
70		* @var array
71		*/
72		private $ort;
73
74		/**
75		* Used for complex scalar division.
76		* @var float
77		*/
78		private $cdivr;
79		private $cdivi;
80
81
82		/**
83		* Symmetric Householder reduction to tridiagonal form.
84		*/
85		private function tred2()
86		{
87		// This is derived from the Algol procedures tred2 by
88		// Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
89		// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
90		// Fortran subroutine in EISPACK.
91		$this->d = $this->V[$this->n-1];
92		// Householder reduction to tridiagonal form.
93		for ($i = $this->n-1; $i > 0; --$i) {
94		$i_ = $i -1;
95		// Scale to avoid under/overflow.
96		$h = $scale = 0.0;
97		$scale += array_sum(array_map('abs', $this->d));
98		if ($scale == 0.0) {
99		$this->e[$i] = $this->d[$i_];
100		$this->d = array_slice($this->V[$i_], 0, $i_);
101		for ($j = 0; $j < $i; ++$j) {
102		$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
103		}
104		} else {
105		// Generate Householder vector.
106		for ($k = 0; $k < $i; ++$k) {
107		$this->d[$k] /= $scale;
108		$h += pow($this->d[$k], 2);
109		}
110		$f = $this->d[$i_];
111		$g = sqrt($h);
112		if ($f > 0) {
113		$g = -$g;
114		}
115		$this->e[$i] = $scale * $g;
116		$h = $h - $f * $g;
117		$this->d[$i_] = $f - $g;
118		for ($j = 0; $j < $i; ++$j) {
119		$this->e[$j] = 0.0;
120		}
121		// Apply similarity transformation to remaining columns.
122		for ($j = 0; $j < $i; ++$j) {
123		$f = $this->d[$j];
124		$this->V[$j][$i] = $f;
125		$g = $this->e[$j] + $this->V[$j][$j] * $f;
126		for ($k = $j+1; $k <= $i_; ++$k) {
127		$g += $this->V[$k][$j] * $this->d[$k];
128		$this->e[$k] += $this->V[$k][$j] * $f;
129		}
130		$this->e[$j] = $g;
131		}
132		$f = 0.0;
133		if ($h === 0 \|\| $h < 1e-32) {
134		$h = 1e-32;
135		}
136		for ($j = 0; $j < $i; ++$j) {
137		$this->e[$j] /= $h;
138		$f += $this->e[$j] * $this->d[$j];
139		}
140		$hh = $f / (2 * $h);
141		for ($j=0; $j < $i; ++$j) {
142		$this->e[$j] -= $hh * $this->d[$j];
143		}
144		for ($j = 0; $j < $i; ++$j) {
145		$f = $this->d[$j];
146		$g = $this->e[$j];
147		for ($k = $j; $k <= $i_; ++$k) {
148		$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
149		}
150		$this->d[$j] = $this->V[$i-1][$j];
151		$this->V[$i][$j] = 0.0;
152		}
153		}
154		$this->d[$i] = $h;
155		}
156
157		// Accumulate transformations.
158		for ($i = 0; $i < $this->n-1; ++$i) {
159		$this->V[$this->n-1][$i] = $this->V[$i][$i];
160		$this->V[$i][$i] = 1.0;
161		$h = $this->d[$i+1];
162		if ($h != 0.0) {
163	View Code Duplication	for ($k = 0; $k <= $i; ++$k) {
164		$this->d[$k] = $this->V[$k][$i+1] / $h;
165		}
166		for ($j = 0; $j <= $i; ++$j) {
167		$g = 0.0;
168	View Code Duplication	for ($k = 0; $k <= $i; ++$k) {
169		$g += $this->V[$k][$i+1] * $this->V[$k][$j];
170		}
171	View Code Duplication	for ($k = 0; $k <= $i; ++$k) {
172		$this->V[$k][$j] -= $g * $this->d[$k];
173		}
174		}
175		}
176		for ($k = 0; $k <= $i; ++$k) {
177		$this->V[$k][$i+1] = 0.0;
178		}
179		}
180
181		$this->d = $this->V[$this->n-1];
182		$this->V[$this->n-1] = array_fill(0, $j, 0.0);
183		$this->V[$this->n-1][$this->n-1] = 1.0;
184		$this->e[0] = 0.0;
185		}
186
187
188		/**
189		* Symmetric tridiagonal QL algorithm.
190		*
191		* This is derived from the Algol procedures tql2, by
192		* Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
193		* Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
194		* Fortran subroutine in EISPACK.
195		*/
196		private function tql2()
197		{
198		for ($i = 1; $i < $this->n; ++$i) {
199		$this->e[$i-1] = $this->e[$i];
200		}
201		$this->e[$this->n-1] = 0.0;
202		$f = 0.0;
203		$tst1 = 0.0;
204		$eps = pow(2.0, -52.0);
205
206		for ($l = 0; $l < $this->n; ++$l) {
207		// Find small subdiagonal element
208		$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
209		$m = $l;
210		while ($m < $this->n) {
211		if (abs($this->e[$m]) <= $eps * $tst1) {
212		break;
213		}
214		++$m;
215		}
216		// If m == l, $this->d[l] is an eigenvalue,
217		// otherwise, iterate.
218		if ($m > $l) {
219		$iter = 0;
220		do {
221		// Could check iteration count here.
222		$iter += 1;
223		// Compute implicit shift
224		$g = $this->d[$l];
225		$p = ($this->d[$l+1] - $g) / (2.0 * $this->e[$l]);
226		$r = hypot($p, 1.0);
227		if ($p < 0) {
228		$r *= -1;
229		}
230		$this->d[$l] = $this->e[$l] / ($p + $r);
231		$this->d[$l+1] = $this->e[$l] * ($p + $r);
232		$dl1 = $this->d[$l+1];
233		$h = $g - $this->d[$l];
234		for ($i = $l + 2; $i < $this->n; ++$i) {
235		$this->d[$i] -= $h;
236		}
237		$f += $h;
238		// Implicit QL transformation.
239		$p = $this->d[$m];
240		$c = 1.0;
241		$c2 = $c3 = $c;
242		$el1 = $this->e[$l + 1];
243		$s = $s2 = 0.0;
244		for ($i = $m-1; $i >= $l; --$i) {
245		$c3 = $c2;
246		$c2 = $c;
247		$s2 = $s;
248		$g = $c * $this->e[$i];
249		$h = $c * $p;
250		$r = hypot($p, $this->e[$i]);
251		$this->e[$i+1] = $s * $r;
252		$s = $this->e[$i] / $r;
253		$c = $p / $r;
254		$p = $c * $this->d[$i] - $s * $g;
255		$this->d[$i+1] = $h + $s * ($c * $g + $s * $this->d[$i]);
256		// Accumulate transformation.
257		for ($k = 0; $k < $this->n; ++$k) {
258		$h = $this->V[$k][$i+1];
259		$this->V[$k][$i+1] = $s * $this->V[$k][$i] + $c * $h;
260		$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
261		}
262		}
263		$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
264		$this->e[$l] = $s * $p;
265		$this->d[$l] = $c * $p;
266		// Check for convergence.
267		} while (abs($this->e[$l]) > $eps * $tst1);
268		}
269		$this->d[$l] = $this->d[$l] + $f;
270		$this->e[$l] = 0.0;
271		}
272
273		// Sort eigenvalues and corresponding vectors.
274		for ($i = 0; $i < $this->n - 1; ++$i) {
275		$k = $i;
276		$p = $this->d[$i];
277		for ($j = $i+1; $j < $this->n; ++$j) {
278		if ($this->d[$j] < $p) {
279		$k = $j;
280		$p = $this->d[$j];
281		}
282		}
283		if ($k != $i) {
284		$this->d[$k] = $this->d[$i];
285		$this->d[$i] = $p;
286	View Code Duplication	for ($j = 0; $j < $this->n; ++$j) {
		0 ignored issues – show Duplication introduced 2017-04-24 14:04 UTC by Report Bug Copy Issue Report Show Similar Issues like this This code seems to be duplicated across your project. Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository. Loading history...
287		$p = $this->V[$j][$i];
288		$this->V[$j][$i] = $this->V[$j][$k];
289		$this->V[$j][$k] = $p;
290		}
291		}
292		}
293		}
294
295
296		/**
297		* Nonsymmetric reduction to Hessenberg form.
298		*
299		* This is derived from the Algol procedures orthes and ortran,
300		* by Martin and Wilkinson, Handbook for Auto. Comp.,
301		* Vol.ii-Linear Algebra, and the corresponding
302		* Fortran subroutines in EISPACK.
303		*/
304		private function orthes()
305		{
306		$low = 0;
307		$high = $this->n-1;
308
309		for ($m = $low+1; $m <= $high-1; ++$m) {
310		// Scale column.
311		$scale = 0.0;
312	View Code Duplication	for ($i = $m; $i <= $high; ++$i) {
313		$scale = $scale + abs($this->H[$i][$m-1]);
314		}
315		if ($scale != 0.0) {
316		// Compute Householder transformation.
317		$h = 0.0;
318		for ($i = $high; $i >= $m; --$i) {
319		$this->ort[$i] = $this->H[$i][$m-1] / $scale;
320		$h += $this->ort[$i] * $this->ort[$i];
321		}
322		$g = sqrt($h);
323		if ($this->ort[$m] > 0) {
324		$g *= -1;
325		}
326		$h -= $this->ort[$m] * $g;
327		$this->ort[$m] -= $g;
328		// Apply Householder similarity transformation
329		// H = (I -u * u' / h) * H * (I -u * u') / h)
330	View Code Duplication	for ($j = $m; $j < $this->n; ++$j) {
331		$f = 0.0;
332		for ($i = $high; $i >= $m; --$i) {
333		$f += $this->ort[$i] * $this->H[$i][$j];
334		}
335		$f /= $h;
336		for ($i = $m; $i <= $high; ++$i) {
337		$this->H[$i][$j] -= $f * $this->ort[$i];
338		}
339		}
340	View Code Duplication	for ($i = 0; $i <= $high; ++$i) {
341		$f = 0.0;
342		for ($j = $high; $j >= $m; --$j) {
343		$f += $this->ort[$j] * $this->H[$i][$j];
344		}
345		$f = $f / $h;
346		for ($j = $m; $j <= $high; ++$j) {
347		$this->H[$i][$j] -= $f * $this->ort[$j];
348		}
349		}
350		$this->ort[$m] = $scale * $this->ort[$m];
351		$this->H[$m][$m-1] = $scale * $g;
352		}
353		}
354
355		// Accumulate transformations (Algol's ortran).
356		for ($i = 0; $i < $this->n; ++$i) {
357		for ($j = 0; $j < $this->n; ++$j) {
358		$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
359		}
360		}
361		for ($m = $high-1; $m >= $low+1; --$m) {
362		if ($this->H[$m][$m-1] != 0.0) {
363	View Code Duplication	for ($i = $m+1; $i <= $high; ++$i) {
364		$this->ort[$i] = $this->H[$i][$m-1];
365		}
366		for ($j = $m; $j <= $high; ++$j) {
367		$g = 0.0;
368		for ($i = $m; $i <= $high; ++$i) {
369		$g += $this->ort[$i] * $this->V[$i][$j];
370		}
371		// Double division avoids possible underflow
372		$g = ($g / $this->ort[$m]) / $this->H[$m][$m-1];
373		for ($i = $m; $i <= $high; ++$i) {
374		$this->V[$i][$j] += $g * $this->ort[$i];
375		}
376		}
377		}
378		}
379		}
380
381
382		/**
383		* Performs complex division.
384		*/
385		private function cdiv($xr, $xi, $yr, $yi)
386		{
387		if (abs($yr) > abs($yi)) {
388		$r = $yi / $yr;
389		$d = $yr + $r * $yi;
390		$this->cdivr = ($xr + $r * $xi) / $d;
391		$this->cdivi = ($xi - $r * $xr) / $d;
392		} else {
393		$r = $yr / $yi;
394		$d = $yi + $r * $yr;
395		$this->cdivr = ($r * $xr + $xi) / $d;
396		$this->cdivi = ($r * $xi - $xr) / $d;
397		}
398		}
399
400
401		/**
402		* Nonsymmetric reduction from Hessenberg to real Schur form.
403		*
404		* Code is derived from the Algol procedure hqr2,
405		* by Martin and Wilkinson, Handbook for Auto. Comp.,
406		* Vol.ii-Linear Algebra, and the corresponding
407		* Fortran subroutine in EISPACK.
408		*/
409		private function hqr2()
410		{
411		// Initialize
412		$nn = $this->n;
413		$n = $nn - 1;
414		$low = 0;
415		$high = $nn - 1;
416		$eps = pow(2.0, -52.0);
417		$exshift = 0.0;
418		$p = $q = $r = $s = $z = 0;
419		// Store roots isolated by balanc and compute matrix norm
420		$norm = 0.0;
421
422		for ($i = 0; $i < $nn; ++$i) {
423		if (($i < $low) or ($i > $high)) {
424		$this->d[$i] = $this->H[$i][$i];
425		$this->e[$i] = 0.0;
426		}
427		for ($j = max($i-1, 0); $j < $nn; ++$j) {
428		$norm = $norm + abs($this->H[$i][$j]);
429		}
430		}
431
432		// Outer loop over eigenvalue index
433		$iter = 0;
434		while ($n >= $low) {
435		// Look for single small sub-diagonal element
436		$l = $n;
437		while ($l > $low) {
438		$s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]);
439		if ($s == 0.0) {
440		$s = $norm;
441		}
442		if (abs($this->H[$l][$l-1]) < $eps * $s) {
443		break;
444		}
445		--$l;
446		}
447		// Check for convergence
448		// One root found
449		if ($l == $n) {
450		$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
451		$this->d[$n] = $this->H[$n][$n];
452		$this->e[$n] = 0.0;
453		--$n;
454		$iter = 0;
455		// Two roots found
456		} elseif ($l == $n-1) {
457		$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
458		$p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0;
459		$q = $p * $p + $w;
460		$z = sqrt(abs($q));
461		$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
462		$this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift;
463		$x = $this->H[$n][$n];
464		// Real pair
465		if ($q >= 0) {
466		if ($p >= 0) {
467		$z = $p + $z;
468		} else {
469		$z = $p - $z;
470		}
471		$this->d[$n-1] = $x + $z;
472		$this->d[$n] = $this->d[$n-1];
473		if ($z != 0.0) {
474		$this->d[$n] = $x - $w / $z;
475		}
476		$this->e[$n-1] = 0.0;
477		$this->e[$n] = 0.0;
478		$x = $this->H[$n][$n-1];
479		$s = abs($x) + abs($z);
480		$p = $x / $s;
481		$q = $z / $s;
482		$r = sqrt($p * $p + $q * $q);
483		$p = $p / $r;
484		$q = $q / $r;
485		// Row modification
486	View Code Duplication	for ($j = $n-1; $j < $nn; ++$j) {
487		$z = $this->H[$n-1][$j];
488		$this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j];
489		$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
490		}
491		// Column modification
492	View Code Duplication	for ($i = 0; $i <= $n; ++$i) {
493		$z = $this->H[$i][$n-1];
494		$this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n];
495		$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
496		}
497		// Accumulate transformations
498	View Code Duplication	for ($i = $low; $i <= $high; ++$i) {
499		$z = $this->V[$i][$n-1];
500		$this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n];
501		$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
502		}
503		// Complex pair
504		} else {
505		$this->d[$n-1] = $x + $p;
506		$this->d[$n] = $x + $p;
507		$this->e[$n-1] = $z;
508		$this->e[$n] = -$z;
509		}
510		$n = $n - 2;
511		$iter = 0;
512		// No convergence yet
513		} else {
514		// Form shift
515		$x = $this->H[$n][$n];
516		$y = 0.0;
517		$w = 0.0;
518		if ($l < $n) {
519		$y = $this->H[$n-1][$n-1];
520		$w = $this->H[$n][$n-1] * $this->H[$n-1][$n];
521		}
522		// Wilkinson's original ad hoc shift
523		if ($iter == 10) {
524		$exshift += $x;
525		for ($i = $low; $i <= $n; ++$i) {
526		$this->H[$i][$i] -= $x;
527		}
528		$s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]);
529		$x = $y = 0.75 * $s;
530		$w = -0.4375 * $s * $s;
531		}
532		// MATLAB's new ad hoc shift
533		if ($iter == 30) {
534		$s = ($y - $x) / 2.0;
535		$s = $s * $s + $w;
536		if ($s > 0) {
537		$s = sqrt($s);
538		if ($y < $x) {
539		$s = -$s;
540		}
541		$s = $x - $w / (($y - $x) / 2.0 + $s);
542		for ($i = $low; $i <= $n; ++$i) {
543		$this->H[$i][$i] -= $s;
544		}
545		$exshift += $s;
546		$x = $y = $w = 0.964;
547		}
548		}
549		// Could check iteration count here.
550		$iter = $iter + 1;
551		// Look for two consecutive small sub-diagonal elements
552		$m = $n - 2;
553		while ($m >= $l) {
554		$z = $this->H[$m][$m];
555		$r = $x - $z;
556		$s = $y - $z;
557		$p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1];
558		$q = $this->H[$m+1][$m+1] - $z - $r - $s;
559		$r = $this->H[$m+2][$m+1];
560		$s = abs($p) + abs($q) + abs($r);
561		$p = $p / $s;
562		$q = $q / $s;
563		$r = $r / $s;
564		if ($m == $l) {
565		break;
566		}
567		if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) <
568		$eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) {
569		break;
570		}
571		--$m;
572		}
573		for ($i = $m + 2; $i <= $n; ++$i) {
574		$this->H[$i][$i-2] = 0.0;
575		if ($i > $m+2) {
576		$this->H[$i][$i-3] = 0.0;
577		}
578		}
579		// Double QR step involving rows l:n and columns m:n
580		for ($k = $m; $k <= $n-1; ++$k) {
581		$notlast = ($k != $n-1);
582		if ($k != $m) {
583		$p = $this->H[$k][$k-1];
584		$q = $this->H[$k+1][$k-1];
585		$r = ($notlast ? $this->H[$k+2][$k-1] : 0.0);
586		$x = abs($p) + abs($q) + abs($r);
587		if ($x != 0.0) {
588		$p = $p / $x;
589		$q = $q / $x;
590		$r = $r / $x;
591		}
592		}
593		if ($x == 0.0) {
594		break;
595		}
596		$s = sqrt($p * $p + $q * $q + $r * $r);
597		if ($p < 0) {
598		$s = -$s;
599		}
600		if ($s != 0) {
601		if ($k != $m) {
602		$this->H[$k][$k-1] = -$s * $x;
603		} elseif ($l != $m) {
604		$this->H[$k][$k-1] = -$this->H[$k][$k-1];
605		}
606		$p = $p + $s;
607		$x = $p / $s;
608		$y = $q / $s;
609		$z = $r / $s;
610		$q = $q / $p;
611		$r = $r / $p;
612		// Row modification
613	View Code Duplication	for ($j = $k; $j < $nn; ++$j) {
614		$p = $this->H[$k][$j] + $q * $this->H[$k+1][$j];
615		if ($notlast) {
616		$p = $p + $r * $this->H[$k+2][$j];
617		$this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z;
618		}
619		$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
620		$this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y;
621		}
622		// Column modification
623	View Code Duplication	for ($i = 0; $i <= min($n, $k+3); ++$i) {
624		$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1];
625		if ($notlast) {
626		$p = $p + $z * $this->H[$i][$k+2];
627		$this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r;
628		}
629		$this->H[$i][$k] = $this->H[$i][$k] - $p;
630		$this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q;
631		}
632		// Accumulate transformations
633	View Code Duplication	for ($i = $low; $i <= $high; ++$i) {
634		$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1];
635		if ($notlast) {
636		$p = $p + $z * $this->V[$i][$k+2];
637		$this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r;
638		}
639		$this->V[$i][$k] = $this->V[$i][$k] - $p;
640		$this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q;
641		}
642		} // ($s != 0)
643		} // k loop
644		} // check convergence
645		} // while ($n >= $low)
646
647		// Backsubstitute to find vectors of upper triangular form
648		if ($norm == 0.0) {
649		return;
650		}
651
652		for ($n = $nn-1; $n >= 0; --$n) {
653		$p = $this->d[$n];
654		$q = $this->e[$n];
655		// Real vector
656		if ($q == 0) {
657		$l = $n;
658		$this->H[$n][$n] = 1.0;
659		for ($i = $n-1; $i >= 0; --$i) {
660		$w = $this->H[$i][$i] - $p;
661		$r = 0.0;
662		for ($j = $l; $j <= $n; ++$j) {
663		$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
664		}
665		if ($this->e[$i] < 0.0) {
666		$z = $w;
667		$s = $r;
668		} else {
669		$l = $i;
670		if ($this->e[$i] == 0.0) {
671		if ($w != 0.0) {
672		$this->H[$i][$n] = -$r / $w;
673		} else {
674		$this->H[$i][$n] = -$r / ($eps * $norm);
675		}
676		// Solve real equations
677		} else {
678		$x = $this->H[$i][$i+1];
679		$y = $this->H[$i+1][$i];
680		$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
681		$t = ($x * $s - $z * $r) / $q;
682		$this->H[$i][$n] = $t;
683		if (abs($x) > abs($z)) {
684		$this->H[$i+1][$n] = (-$r - $w * $t) / $x;
685		} else {
686		$this->H[$i+1][$n] = (-$s - $y * $t) / $z;
687		}
688		}
689		// Overflow control
690		$t = abs($this->H[$i][$n]);
691		if (($eps * $t) * $t > 1) {
692	View Code Duplication	for ($j = $i; $j <= $n; ++$j) {
693		$this->H[$j][$n] = $this->H[$j][$n] / $t;
694		}
695		}
696		}
697		}
698		// Complex vector
699		} elseif ($q < 0) {
700		$l = $n-1;
701		// Last vector component imaginary so matrix is triangular
702		if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) {
703		$this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1];
704		$this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1];
705		} else {
706		$this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q);
707		$this->H[$n-1][$n-1] = $this->cdivr;
708		$this->H[$n-1][$n] = $this->cdivi;
709		}
710		$this->H[$n][$n-1] = 0.0;
711		$this->H[$n][$n] = 1.0;
712		for ($i = $n-2; $i >= 0; --$i) {
713		// double ra,sa,vr,vi;
714		$ra = 0.0;
715		$sa = 0.0;
716		for ($j = $l; $j <= $n; ++$j) {
717		$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1];
718		$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
719		}
720		$w = $this->H[$i][$i] - $p;
721		if ($this->e[$i] < 0.0) {
722		$z = $w;
723		$r = $ra;
724		$s = $sa;
725		} else {
726		$l = $i;
727		if ($this->e[$i] == 0) {
728		$this->cdiv(-$ra, -$sa, $w, $q);
729		$this->H[$i][$n-1] = $this->cdivr;
730		$this->H[$i][$n] = $this->cdivi;
731		} else {
732		// Solve complex equations
733		$x = $this->H[$i][$i+1];
734		$y = $this->H[$i+1][$i];
735		$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
736		$vi = ($this->d[$i] - $p) * 2.0 * $q;
737		if ($vr == 0.0 & $vi == 0.0) {
738		$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
739		}
740		$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
741		$this->H[$i][$n-1] = $this->cdivr;
742		$this->H[$i][$n] = $this->cdivi;
743		if (abs($x) > (abs($z) + abs($q))) {
744		$this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x;
745		$this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x;
746		} else {
747		$this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q);
748		$this->H[$i+1][$n-1] = $this->cdivr;
749		$this->H[$i+1][$n] = $this->cdivi;
750		}
751		}
752		// Overflow control
753		$t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n]));
754		if (($eps * $t) * $t > 1) {
755		for ($j = $i; $j <= $n; ++$j) {
756		$this->H[$j][$n-1] = $this->H[$j][$n-1] / $t;
757		$this->H[$j][$n] = $this->H[$j][$n] / $t;
758		}
759		}
760		} // end else
761		} // end for
762		} // end else for complex case
763		} // end for
764
765		// Vectors of isolated roots
766		for ($i = 0; $i < $nn; ++$i) {
767		if ($i < $low \| $i > $high) {
768		for ($j = $i; $j < $nn; ++$j) {
769		$this->V[$i][$j] = $this->H[$i][$j];
770		}
771		}
772		}
773
774		// Back transformation to get eigenvectors of original matrix
775		for ($j = $nn-1; $j >= $low; --$j) {
776		for ($i = $low; $i <= $high; ++$i) {
777		$z = 0.0;
778		for ($k = $low; $k <= min($j, $high); ++$k) {
779		$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
780		}
781		$this->V[$i][$j] = $z;
782		}
783		}
784		} // end hqr2
785
786
787		/**
788		* Constructor: Check for symmetry, then construct the eigenvalue decomposition
789		*
790		* @param array $Arg
791		*/
792		public function __construct(array $Arg)
793		{
794		$this->A = $Arg;
795		$this->n = count($Arg[0]);
796
797		$issymmetric = true;
798		for ($j = 0; ($j < $this->n) & $issymmetric; ++$j) {
799	View Code Duplication	for ($i = 0; ($i < $this->n) & $issymmetric; ++$i) {
800		$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
801		}
802		}
803
804		if ($issymmetric) {
805		$this->V = $this->A;
806		// Tridiagonalize.
807		$this->tred2();
808		// Diagonalize.
809		$this->tql2();
810		} else {
811		$this->H = $this->A;
812		$this->ort = [];
813		// Reduce to Hessenberg form.
814		$this->orthes();
815		// Reduce Hessenberg to real Schur form.
816		$this->hqr2();
817		}
818		}
819
820		/**
821		* Return the eigenvector matrix
822		*
823		* @access public
824		* @return array
825		*/
826		public function getEigenvectors()
827		{
828		$vectors = $this->V;
829
830		// Always return the eigenvectors of length 1.0
831		$vectors = new Matrix($vectors);
832		$vectors = array_map(function ($vect) {
833		$sum = 0;
834	View Code Duplication	for ($i=0; $i<count($vect); $i++) {
835		$sum += $vect[$i] ** 2;
836		}
837		$sum = sqrt($sum);
838	View Code Duplication	for ($i=0; $i<count($vect); $i++) {
839		$vect[$i] /= $sum;
840		}
841		return $vect;
842		}, $vectors->transpose()->toArray());
843
844		return $vectors;
845		}
846
847
848		/**
849		* Return the real parts of the eigenvalues<br>
850		* d = real(diag(D));
851		*
852		* @return array
853		*/
854		public function getRealEigenvalues()
855		{
856		return $this->d;
857		}
858
859
860		/**
861		* Return the imaginary parts of the eigenvalues <br>
862		* d = imag(diag(D))
863		*
864		* @return array
865		*/
866		public function getImagEigenvalues()
867		{
868		return $this->e;
869		}
870
871
872		/**
873		* Return the block diagonal eigenvalue matrix
874		*
875		* @return array
876		*/
877		public function getDiagonalEigenvalues()
878		{
879		for ($i = 0; $i < $this->n; ++$i) {
880		$D[$i] = array_fill(0, $this->n, 0.0);
881		$D[$i][$i] = $this->d[$i];
882		if ($this->e[$i] == 0) {
883		continue;
884		}
885		$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
886		$D[$i][$o] = $this->e[$i];
887		}
888		return $D;
889		}
890		} // class EigenvalueDecomposition
891

php-ai / php-ml

Push — master ( 12b8b1...5b373f )

Math/LinearAlgebra/EigenvalueDecomposition.php (1 issue)

Labels

Severity

Introduced By

Upgrade to new PHP Analysis Engine

Duplication Side-by-Side

Filter issues like