WebSVN – Moodle – Autoría – /lib/mlbackend/php/phpml/src/Phpml/Math/LinearAlgebra/EigenvalueDecomposition.php

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

Proyectos de Subversion Moodle

(root)/lib/mlbackend/php/phpml/src/Phpml/Math/LinearAlgebra/EigenvalueDecomposition.php – Rev 1