diff --git a/Math/Matrix/QRDecomposition.php b/Math/Matrix/QRDecomposition.php index 25459b0ea..5c6dbbb66 100644 --- a/Math/Matrix/QRDecomposition.php +++ b/Math/Matrix/QRDecomposition.php @@ -83,7 +83,7 @@ final class QRDecomposition $nrm = Triangle::getHypot($nrm, $this->QR[$i][$k]); } - if ($nrm != 0.0) { + if ($nrm != 0) { // Form k-th Householder vector. if ($this->QR[$k][$k] < 0) { $nrm = -$nrm; diff --git a/Math/Matrix/SingularValueDecomposition.php b/Math/Matrix/SingularValueDecomposition.php index b96a13157..14e113912 100644 --- a/Math/Matrix/SingularValueDecomposition.php +++ b/Math/Matrix/SingularValueDecomposition.php @@ -24,4 +24,109 @@ namespace phpOMS\Math\Matrix; */ final class SingularValueDecomposition { + /** + * U matrix. + * + * @var array[] + * @since 1.0.0 + */ + private $U = []; + + /** + * V matrix. + * + * @var array[] + * @since 1.0.0 + */ + private $V = []; + + /** + * Singular values. + * + * @var array + * @since 1.0.0 + */ + private $S = []; + + /** + * Dimension m + * + * @var int + * @since 1.0.0 + */ + private $m = 0; + + /** + * Dimension n + * + * @var int + * @since 1.0.0 + */ + private $n = 0; + + public function getU() : Matrix + { + $matrix = new Matrix(); + $matrix->setMatrix($this->U); + + return $matrix; + } + + public function getV() : Matrix + { + $matrix = new Matrix(); + $matrix->setMatrix($this->V); + + return $matrix; + } + + public function getS() : Matrix + { + $S = [[]]; + for ($i = 0; $i < $this->n; ++$i) { + 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); + + return $matrix; + } + + public function getSingularValues() : Vector + { + $vector = new Vector(); + $vector->setMatrix($this->S); + + return $vector; + } + + public function norm2() : float + { + return $this->S[0]; + } + + public function cond() : float + { + return $this->S[0] / $this->S[\min($this->m, $this->n) - 1]; + } + + public function rank() : int + { + $eps = 0.00001; + $tol = \max($this->m, $this->n) * $this->S[0] * $eps; + $r = 0; + + $length = \count($this->S); + for ($i = 0; $i < $length; ++$i) { + if ($this->S[$i] > $tol) { + ++$r; + } + } + + return $r; + } } diff --git a/tests/Math/Matrix/QRDecompositionTest.php b/tests/Math/Matrix/QRDecompositionTest.php index f6f7ddc45..46af3ef79 100644 --- a/tests/Math/Matrix/QRDecompositionTest.php +++ b/tests/Math/Matrix/QRDecompositionTest.php @@ -19,7 +19,7 @@ use phpOMS\Math\Matrix\QRDecomposition; class QRDecompositionTest extends \PHPUnit\Framework\TestCase { - /**public function testDecomposition() + public function testDecomposition() { $A = new Matrix(); $A->setMatrix([ @@ -32,24 +32,26 @@ class QRDecompositionTest extends \PHPUnit\Framework\TestCase self::assertTrue($QR->isFullRank()); + self::assertEquals([ + [-6 / 7, 69 / 175, 58 / 175], + [-3 / 7, -158 / 175, -6 / 175], + [2 / 7, -6 / 35, 33 / 35], + ], $QR->getQ()->toArray(), '', 0.2); + + self::assertEquals([ + [-14, -21, 14], + [0, -175, 70], + [0, 0, -35], + ], $QR->getR()->toArray(), '', 0.2); + + self::assertEquals($A->toArray(), $QR->getQ()->mult($QR->getR()), '', 0.2); + self::assertEquals([ [12, -69, -58 / 5], [6, 158, 6 / 5], [-4, 30, -33], ], $QR->getH()->toArray(), '', 0.2); - - self::assertEquals([ - [6 / 7, -69 / 175, -58 / 175], - [3 / 7, 158 / 175, 6 / 175], - [-2 / 7, 6 / 35, -33 / 35], - ], $QR->getQ()->toArray(), '', 0.2); - - self::assertEquals([ - [14, 21, -14], - [0, 175, -70], - [0, 0, 35], - ], $QR->getR()->toArray(), '', 0.2); - }*/ + } public function testSolve() { diff --git a/tests/Math/Matrix/SingularValueDecompositionTest.php b/tests/Math/Matrix/SingularValueDecompositionTest.php new file mode 100644 index 000000000..b394f11c9 --- /dev/null +++ b/tests/Math/Matrix/SingularValueDecompositionTest.php @@ -0,0 +1,26 @@ +