From d33f24aa6b4675a6eea2e59c75ebdec4605c205a Mon Sep 17 00:00:00 2001 From: Dennis Eichhorn Date: Sun, 3 Mar 2019 22:03:59 +0100 Subject: [PATCH] Fix matrix bugs --- Math/Matrix/EigenvalueDecomposition.php | 22 ++--- Math/Matrix/SingularValueDecomposition.php | 95 +++++++++++++++---- .../Matrix/EigenvalueDecompositionTest.php | 8 +- .../Matrix/SingularValueDecompositionTest.php | 2 +- 4 files changed, 93 insertions(+), 34 deletions(-) diff --git a/Math/Matrix/EigenvalueDecomposition.php b/Math/Matrix/EigenvalueDecomposition.php index 6e5ca0925..5a621d87b 100644 --- a/Math/Matrix/EigenvalueDecomposition.php +++ b/Math/Matrix/EigenvalueDecomposition.php @@ -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) { for ($i = $m + 1; $i <= $high; ++$i) { $this->ort[$i] = $this->H[$i][$m - 1]; @@ -463,7 +463,7 @@ final class EigenvalueDecomposition $d = $yr + $r * $yi; $this->cdivr = ($xr + $r * $xi) / $d; - $this->cdivi = ($xi + $r * $xr) / $d; + $this->cdivi = ($xi - $r * $xr) / $d; } else { $r = $yr / $yi; $d = $yi + $r * $yr; @@ -479,7 +479,7 @@ final class EigenvalueDecomposition $n = $nn - 1; $low = 0; $high = $nn - 1; - $eps = 0.0001; + $eps = 0.00001; $exshift = 0.0; $p = 0; $q = 0; @@ -654,15 +654,13 @@ final class EigenvalueDecomposition $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0); $x = \abs($p) + \abs($q) + \abs($r); - if ($x != 0) { - $p /= $x; - $q /= $x; - $r /= $x; + if ($x == 0) { + continue; } - } - if ($x == 0) { - break; + $p /= $x; + $q /= $x; + $r /= $x; } $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]; 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] = $this->V[$i][$k] - $p; @@ -853,7 +851,7 @@ final class EigenvalueDecomposition $min = \min($j, $high); 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; diff --git a/Math/Matrix/SingularValueDecomposition.php b/Math/Matrix/SingularValueDecomposition.php index 03f794606..8474efce9 100644 --- a/Math/Matrix/SingularValueDecomposition.php +++ b/Math/Matrix/SingularValueDecomposition.php @@ -66,6 +66,13 @@ final class SingularValueDecomposition */ private $n = 0; + /** + * Constructor. + * + * @param Matrix $M Matrix + * + * @since 1.0.0 + */ public function __construct(Matrix $M) { $A = $M->toArray(); @@ -78,7 +85,9 @@ final class SingularValueDecomposition $nrt = \max(0, \min($this->n - 2, $this->m)); $eps = 0.00001; - for ($k = 0; $k < \max($nct, $nrt); ++$k) { + $maxNctNrt = \max($nct, $nrt); + + for ($k = 0; $k < $maxNctNrt; ++$k) { if ($k < $nct) { $this->S[$k] = 0; for ($i = $k; $i < $this->m; ++$i) { @@ -207,7 +216,7 @@ final class SingularValueDecomposition $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) { $this->U[$i][$k] = 0.0; } @@ -245,7 +254,7 @@ final class SingularValueDecomposition $pp = $p - 1; $iter = 0; - while ($p > 0) { + while (true) { for ($k = $p - 2; $k >= -1; --$k) { if ($k === -1) { break; @@ -257,8 +266,9 @@ final class SingularValueDecomposition } } + $case = 0; if ($k === $p - 2) { - $kase = 4; + $case = 4; } else { for ($ks = $p - 1; $ks >= $k; --$ks) { 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); + if (\abs($this->S[$ks]) <= $eps * $t) { $this->S[$ks] = 0.0; break; @@ -273,23 +284,23 @@ final class SingularValueDecomposition } if ($ks === $k) { - $kase = 3; + $case = 3; } elseif ($ks === $p - 1) { - $kase = 1; + $case = 1; } else { - $kase = 2; + $case = 2; $k = $ks; } } ++$k; - switch ($kase) { + switch ($case) { case 1: $f = $e[$p - 2]; $e[$p - 2] = 0.0; 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; $sn = $f / $t; $this->S[$j] = $t; @@ -362,7 +373,7 @@ final class SingularValueDecomposition $g = $sk * $ek; for ($j = $k; $j < $p - 1; ++$j) { - $t = Triangle::getHypot($f,$g); + $t = Triangle::getHypot($f, $g); $cs = $f / $t; $sn = $g / $t; @@ -370,10 +381,10 @@ final class SingularValueDecomposition $e[$j - 1] = $t; } - $f = $cs * $this->S[$j] + $sn * $e[$j]; - $e[$j] = $cs * $e[$j] - $sn * $this->S[$j]; - $g = $sn * $this->S[$j + 1]; - $this->S[$j + 1] = $cs * $this->S[$j + 1]; + $f = $cs * $this->S[$j] + $sn * $e[$j]; + $e[$j] = $cs * $e[$j] - $sn * $this->S[$j]; + $g = $sn * $this->S[$j + 1]; + $this->S[$j + 1] *= $cs; for ($i = 0; $i < $this->n; ++$i) { $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j + 1]; @@ -381,7 +392,7 @@ final class SingularValueDecomposition $this->V[$i][$j] = $t; } - $t = Triangle::getHypot($f,$g); + $t = Triangle::getHypot($f, $g); $cs = $f / $t; $sn = $g / $t; $this->S[$j] = $t; @@ -400,7 +411,7 @@ final class SingularValueDecomposition } $e[$p - 2] = $f; - $iter = $iter + 1; + ++$iter; break; case 4: 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 { $matrix = new Matrix(); @@ -453,6 +471,13 @@ final class SingularValueDecomposition return $matrix; } + /** + * Get V matrix + * + * @return Matrix + * + * @since 1.0.0 + */ public function getV() : Matrix { $matrix = new Matrix(); @@ -461,6 +486,13 @@ final class SingularValueDecomposition return $matrix; } + /** + * Get S matrix + * + * @return Matrix + * + * @since 1.0.0 + */ public function getS() : Matrix { $S = [[]]; @@ -468,15 +500,23 @@ final class SingularValueDecomposition for ($j = 0; $j < $this->n; ++$j) { $S[$i][$j] = 0.0; } + $S[$i][$i] = $this->S[$i]; } $matrix = new Matrix(); - $matrix->setMatrix($this->V); + $matrix->setMatrix($S); return $matrix; } + /** + * Get singular Values + * + * @return Vector + * + * @since 1.0.0 + */ public function getSingularValues() : Vector { $vector = new Vector(); @@ -485,16 +525,37 @@ final class SingularValueDecomposition return $vector; } + /** + * Get norm + * + * @return float + * + * @since 1.0.0 + */ public function norm2() : float { return $this->S[0]; } + /** + * Get condition + * + * @return float + * + * @since 1.0.0 + */ public function cond() : float { return $this->S[0] / $this->S[\min($this->m, $this->n) - 1]; } + /** + * Get rank + * + * @return int + * + * @since 1.0.0 + */ public function rank() : int { $eps = 0.00001; diff --git a/tests/Math/Matrix/EigenvalueDecompositionTest.php b/tests/Math/Matrix/EigenvalueDecompositionTest.php index aec553e9e..a39e63723 100644 --- a/tests/Math/Matrix/EigenvalueDecompositionTest.php +++ b/tests/Math/Matrix/EigenvalueDecompositionTest.php @@ -60,9 +60,9 @@ class EigenvalueDecompositionTest extends \PHPUnit\Framework\TestCase self::assertEquals([-5, 3, 6], $eig->getRealEigenvalues()->toArray(), '', 0.2); self::assertEquals([ - [sqrt(2/3), sqrt(2/7), 1/sqrt(293)], - [-1/sqrt(6), -3/sqrt(14), 6/sqrt(293)], - [1/sqrt(6), -1/sqrt(14), 16/sqrt(293)], + [-sqrt(2/3), sqrt(2/7), -1/sqrt(293)], + [-1/sqrt(6), -3/sqrt(14), -6/sqrt(293)], + [1/sqrt(6), -1/sqrt(14), -16/sqrt(293)], ], $eig->getV()->toArray(), '', 0.2); self::assertEquals([ @@ -107,7 +107,7 @@ class EigenvalueDecompositionTest extends \PHPUnit\Framework\TestCase $A->toArray(), $eig->getV() ->mult($eig->getD()) - ->mult($eig->getV()->transpose()) + ->mult($eig->getV()->inverse()) ->toArray(), '', 0.2 ); diff --git a/tests/Math/Matrix/SingularValueDecompositionTest.php b/tests/Math/Matrix/SingularValueDecompositionTest.php index 4b72ceb0c..3b417c72b 100644 --- a/tests/Math/Matrix/SingularValueDecompositionTest.php +++ b/tests/Math/Matrix/SingularValueDecompositionTest.php @@ -32,7 +32,7 @@ class SingularValueDecompositionTest extends \PHPUnit\Framework\TestCase self::assertEquals([ [-0.3092, -0.9510], - [-0.9510, -0.3092], + [-0.9510, 0.3092], ], $svd->getU()->toArray(), '', 0.2); self::assertEquals([