Jingga/phpOMS
1
0
Fork 0
mirror of https://github.com/Karaka-Management/phpOMS.git synced 2026-02-09 13:38:41 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Fix matrix bugs

Browse Source
This commit is contained in:
Dennis Eichhorn 2019-03-03 22:03:59 +01:00
parent 169d361522
commit d33f24aa6b
4 changed files with 93 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Math/Matrix/EigenvalueDecomposition.php
View File

@ -432,7 +432,7 @@ final class EigenvalueDecomposition
} }
} }
for ($m = $high - 1; $m >= $low + 1; --$m) { for ($m = $high - 1; $m > $low; --$m) {
if ($this->H[$m][$m - 1] != 0) { if ($this->H[$m][$m - 1] != 0) {
for ($i = $m + 1; $i <= $high; ++$i) { for ($i = $m + 1; $i <= $high; ++$i) {
$this->ort[$i] = $this->H[$i][$m - 1]; $this->ort[$i] = $this->H[$i][$m - 1];
@ -463,7 +463,7 @@ final class EigenvalueDecomposition
$d = $yr + $r * $yi; $d = $yr + $r * $yi;
$this->cdivr = ($xr + $r * $xi) / $d; $this->cdivr = ($xr + $r * $xi) / $d;
$this->cdivi = ($xi + $r * $xr) / $d; $this->cdivi = ($xi - $r * $xr) / $d;
} else { } else {
$r = $yr / $yi; $r = $yr / $yi;
$d = $yi + $r * $yr; $d = $yi + $r * $yr;
@ -479,7 +479,7 @@ final class EigenvalueDecomposition
$n = $nn - 1; $n = $nn - 1;
$low = 0; $low = 0;
$high = $nn - 1; $high = $nn - 1;
$eps = 0.0001; $eps = 0.00001;
$exshift = 0.0; $exshift = 0.0;
$p = 0; $p = 0;
$q = 0; $q = 0;
@ -654,15 +654,13 @@ final class EigenvalueDecomposition
$r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
$x = \abs($p) + \abs($q) + \abs($r); $x = \abs($p) + \abs($q) + \abs($r);
if ($x != 0) { if ($x == 0) {
$p /= $x; continue;
$q /= $x;
$r /= $x;
} }
}
if ($x == 0) { $p /= $x;
break; $q /= $x;
$r /= $x;
} }
$s = $p < 0 ? -\sqrt($p * $p + $q * $q + $r * $r) : \sqrt($p * $p + $q * $q + $r * $r); $s = $p < 0 ? -\sqrt($p * $p + $q * $q + $r * $r) : \sqrt($p * $p + $q * $q + $r * $r);
@ -709,7 +707,7 @@ final class EigenvalueDecomposition
$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1]; $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
if ($notlast) { if ($notlast) {
$p = $p + $z * $this->V[$i][$k + 2]; $p += $z * $this->V[$i][$k + 2];
$this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r; $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
} }
$this->V[$i][$k] = $this->V[$i][$k] - $p; $this->V[$i][$k] = $this->V[$i][$k] - $p;
@ -853,7 +851,7 @@ final class EigenvalueDecomposition
$min = \min($j, $high); $min = \min($j, $high);
for ($k = $low; $k <= $min; ++$k) { for ($k = $low; $k <= $min; ++$k) {
$z = $z + $this->V[$i][$k] * $this->H[$k][$j]; $z += $this->V[$i][$k] * $this->H[$k][$j];
} }
$this->V[$i][$j] = $z; $this->V[$i][$j] = $z;

95
Math/Matrix/SingularValueDecomposition.php
View File

@ -66,6 +66,13 @@ final class SingularValueDecomposition
*/ */
private $n = 0; private $n = 0;
/**
* Constructor.
*
* @param Matrix $M Matrix
*
* @since 1.0.0
*/
public function __construct(Matrix $M) public function __construct(Matrix $M)
{ {
$A = $M->toArray(); $A = $M->toArray();
@ -78,7 +85,9 @@ final class SingularValueDecomposition
$nrt = \max(0, \min($this->n - 2, $this->m)); $nrt = \max(0, \min($this->n - 2, $this->m));
$eps = 0.00001; $eps = 0.00001;
for ($k = 0; $k < \max($nct, $nrt); ++$k) { $maxNctNrt = \max($nct, $nrt);
for ($k = 0; $k < $maxNctNrt; ++$k) {
if ($k < $nct) { if ($k < $nct) {
$this->S[$k] = 0; $this->S[$k] = 0;
for ($i = $k; $i < $this->m; ++$i) { for ($i = $k; $i < $this->m; ++$i) {
@ -207,7 +216,7 @@ final class SingularValueDecomposition
$this->U[$i][$k] = -$this->U[$i][$k]; $this->U[$i][$k] = -$this->U[$i][$k];
} }
$this->U[$k][$k] = 1.0 + $this->U[$k][$k]; $this->U[$k][$k] += 1.0;
for ($i = 0; $i < $k - 1; ++$i) { for ($i = 0; $i < $k - 1; ++$i) {
$this->U[$i][$k] = 0.0; $this->U[$i][$k] = 0.0;
} }
@ -245,7 +254,7 @@ final class SingularValueDecomposition
$pp = $p - 1; $pp = $p - 1;
$iter = 0; $iter = 0;
while ($p > 0) { while (true) {
for ($k = $p - 2; $k >= -1; --$k) { for ($k = $p - 2; $k >= -1; --$k) {
if ($k === -1) { if ($k === -1) {
break; break;
@ -257,8 +266,9 @@ final class SingularValueDecomposition
} }
} }
$case = 0;
if ($k === $p - 2) { if ($k === $p - 2) {
$kase = 4; $case = 4;
} else { } else {
for ($ks = $p - 1; $ks >= $k; --$ks) { for ($ks = $p - 1; $ks >= $k; --$ks) {
if ($ks === $k) { if ($ks === $k) {
@ -266,6 +276,7 @@ final class SingularValueDecomposition
} }
$t = ($ks !== $p ? \abs($e[$ks]) : 0) + ($ks !== $k + 1 ? \abs($e[$ks - 1]) : 0); $t = ($ks !== $p ? \abs($e[$ks]) : 0) + ($ks !== $k + 1 ? \abs($e[$ks - 1]) : 0);
if (\abs($this->S[$ks]) <= $eps * $t) { if (\abs($this->S[$ks]) <= $eps * $t) {
$this->S[$ks] = 0.0; $this->S[$ks] = 0.0;
break; break;
@ -273,23 +284,23 @@ final class SingularValueDecomposition
} }
if ($ks === $k) { if ($ks === $k) {
$kase = 3; $case = 3;
} elseif ($ks === $p - 1) { } elseif ($ks === $p - 1) {
$kase = 1; $case = 1;
} else { } else {
$kase = 2; $case = 2;
$k = $ks; $k = $ks;
} }
} }
++$k; ++$k;
switch ($kase) { switch ($case) {
case 1: case 1:
$f = $e[$p - 2]; $f = $e[$p - 2];
$e[$p - 2] = 0.0; $e[$p - 2] = 0.0;
for ($j = $p - 2; $j >= $k; --$j) { for ($j = $p - 2; $j >= $k; --$j) {
$t = Triangle::getHypot($this->S[$j],$f); $t = Triangle::getHypot($this->S[$j], $f);
$cs = $this->S[$j] / $t; $cs = $this->S[$j] / $t;
$sn = $f / $t; $sn = $f / $t;
$this->S[$j] = $t; $this->S[$j] = $t;
@ -362,7 +373,7 @@ final class SingularValueDecomposition
$g = $sk * $ek; $g = $sk * $ek;
for ($j = $k; $j < $p - 1; ++$j) { for ($j = $k; $j < $p - 1; ++$j) {
$t = Triangle::getHypot($f,$g); $t = Triangle::getHypot($f, $g);
$cs = $f / $t; $cs = $f / $t;
$sn = $g / $t; $sn = $g / $t;
@ -370,10 +381,10 @@ final class SingularValueDecomposition
$e[$j - 1] = $t; $e[$j - 1] = $t;
} }
$f = $cs * $this->S[$j] + $sn * $e[$j]; $f = $cs * $this->S[$j] + $sn * $e[$j];
$e[$j] = $cs * $e[$j] - $sn * $this->S[$j]; $e[$j] = $cs * $e[$j] - $sn * $this->S[$j];
$g = $sn * $this->S[$j + 1]; $g = $sn * $this->S[$j + 1];
$this->S[$j + 1] = $cs * $this->S[$j + 1]; $this->S[$j + 1] *= $cs;
for ($i = 0; $i < $this->n; ++$i) { for ($i = 0; $i < $this->n; ++$i) {
$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1]; $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1];
@ -381,7 +392,7 @@ final class SingularValueDecomposition
$this->V[$i][$j] = $t; $this->V[$i][$j] = $t;
} }
$t = Triangle::getHypot($f,$g); $t = Triangle::getHypot($f, $g);
$cs = $f / $t; $cs = $f / $t;
$sn = $g / $t; $sn = $g / $t;
$this->S[$j] = $t; $this->S[$j] = $t;
@ -400,7 +411,7 @@ final class SingularValueDecomposition
} }
$e[$p - 2] = $f; $e[$p - 2] = $f;
$iter = $iter + 1; ++$iter;
break; break;
case 4: case 4:
if ($this->S[$k] <= 0.0) { if ($this->S[$k] <= 0.0) {
@ -445,6 +456,13 @@ final class SingularValueDecomposition
} }
} }
/**
* Get U matrix
*
* @return Matrix
*
* @since 1.0.0
*/
public function getU() : Matrix public function getU() : Matrix
{ {
$matrix = new Matrix(); $matrix = new Matrix();
@ -453,6 +471,13 @@ final class SingularValueDecomposition
return $matrix; return $matrix;
} }
/**
* Get V matrix
*
* @return Matrix
*
* @since 1.0.0
*/
public function getV() : Matrix public function getV() : Matrix
{ {
$matrix = new Matrix(); $matrix = new Matrix();
@ -461,6 +486,13 @@ final class SingularValueDecomposition
return $matrix; return $matrix;
} }
/**
* Get S matrix
*
* @return Matrix
*
* @since 1.0.0
*/
public function getS() : Matrix public function getS() : Matrix
{ {
$S = [[]]; $S = [[]];
@ -468,15 +500,23 @@ final class SingularValueDecomposition
for ($j = 0; $j < $this->n; ++$j) { for ($j = 0; $j < $this->n; ++$j) {
$S[$i][$j] = 0.0; $S[$i][$j] = 0.0;
} }
$S[$i][$i] = $this->S[$i]; $S[$i][$i] = $this->S[$i];
} }
$matrix = new Matrix(); $matrix = new Matrix();
$matrix->setMatrix($this->V); $matrix->setMatrix($S);
return $matrix; return $matrix;
} }
/**
* Get singular Values
*
* @return Vector
*
* @since 1.0.0
*/
public function getSingularValues() : Vector public function getSingularValues() : Vector
{ {
$vector = new Vector(); $vector = new Vector();
@ -485,16 +525,37 @@ final class SingularValueDecomposition
return $vector; return $vector;
} }
/**
* Get norm
*
* @return float
*
* @since 1.0.0
*/
public function norm2() : float public function norm2() : float
{ {
return $this->S[0]; return $this->S[0];
} }
/**
* Get condition
*
* @return float
*
* @since 1.0.0
*/
public function cond() : float public function cond() : float
{ {
return $this->S[0] / $this->S[\min($this->m, $this->n) - 1]; return $this->S[0] / $this->S[\min($this->m, $this->n) - 1];
} }
/**
* Get rank
*
* @return int
*
* @since 1.0.0
*/
public function rank() : int public function rank() : int
{ {
$eps = 0.00001; $eps = 0.00001;

8
tests/Math/Matrix/EigenvalueDecompositionTest.php
View File

@ -60,9 +60,9 @@ class EigenvalueDecompositionTest extends \PHPUnit\Framework\TestCase
self::assertEquals([-5, 3, 6], $eig->getRealEigenvalues()->toArray(), '', 0.2); self::assertEquals([-5, 3, 6], $eig->getRealEigenvalues()->toArray(), '', 0.2);
self::assertEquals([ self::assertEquals([
[sqrt(2/3), sqrt(2/7), 1/sqrt(293)], [-sqrt(2/3), sqrt(2/7), -1/sqrt(293)],
[-1/sqrt(6), -3/sqrt(14), 6/sqrt(293)], [-1/sqrt(6), -3/sqrt(14), -6/sqrt(293)],
[1/sqrt(6), -1/sqrt(14), 16/sqrt(293)], [1/sqrt(6), -1/sqrt(14), -16/sqrt(293)],
], $eig->getV()->toArray(), '', 0.2); ], $eig->getV()->toArray(), '', 0.2);
self::assertEquals([ self::assertEquals([
@ -107,7 +107,7 @@ class EigenvalueDecompositionTest extends \PHPUnit\Framework\TestCase
$A->toArray(), $A->toArray(),
$eig->getV() $eig->getV()
->mult($eig->getD()) ->mult($eig->getD())
->mult($eig->getV()->transpose()) ->mult($eig->getV()->inverse())
->toArray(), ->toArray(),
'', 0.2 '', 0.2
); );

2
tests/Math/Matrix/SingularValueDecompositionTest.php
View File

@ -32,7 +32,7 @@ class SingularValueDecompositionTest extends \PHPUnit\Framework\TestCase
self::assertEquals([ self::assertEquals([
[-0.3092, -0.9510], [-0.3092, -0.9510],
[-0.9510, -0.3092], [-0.9510, 0.3092],
], $svd->getU()->toArray(), '', 0.2); ], $svd->getU()->toArray(), '', 0.2);
self::assertEquals([ self::assertEquals([