From 39bb8c83d40c14baa5106996ee9ae77a033ad5c8 Mon Sep 17 00:00:00 2001 From: Dennis Eichhorn Date: Wed, 15 Nov 2017 13:04:43 +0100 Subject: [PATCH] Improve matrix implementations --- Math/Matrix/CholeskyDecomposition.php | 66 +++++++++++++++------------ Math/Matrix/LUDecomposition.php | 36 +++++++++------ Math/Matrix/Vector.php | 4 +- 3 files changed, 62 insertions(+), 44 deletions(-) diff --git a/Math/Matrix/CholeskyDecomposition.php b/Math/Matrix/CholeskyDecomposition.php index 050612678..cae3d41bd 100644 --- a/Math/Matrix/CholeskyDecomposition.php +++ b/Math/Matrix/CholeskyDecomposition.php @@ -21,8 +21,12 @@ class CholeskyDecomposition private $m = 0; + /** + * Is symmetric positive definite + */ private $isSpd = true; + // see http://www.aip.de/groups/soe/local/numres/bookfpdf/f2-9.pdf public function __construct(Matrix $M) { $this->L = $M->toArray(); @@ -33,14 +37,15 @@ class CholeskyDecomposition for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) { $sum -= $this->L[$i][$k] * $this->L[$j][$k]; } - if ($i == $j) { + + if ($i === $j) { if ($sum >= 0) { $this->L[$i][$i] = sqrt($sum); } else { $this->isSpd = false; } } else { - if ($this->L[$i][$i] != 0) { + if ($this->L[$i][$i] !== 0) { $this->L[$j][$i] = $sum / $this->L[$i][$i]; } } @@ -52,12 +57,12 @@ class CholeskyDecomposition } } - public function isSpd() + public function isSpd() : bool { return $this->isSpd; } - public function getL() + public function getL() : Matrix { $matrix = new Matrix(); $matrix->setMatrix($this->L); @@ -65,41 +70,44 @@ class CholeskyDecomposition return $matrix; } - public function solve(Matrix $B) + public function solve(Matrix $B) : Matrix { if ($B->getM() !== $this->m) { - // invalid dimension + throw new \Exception(); } if (!$this->isSpd) { - // is not positive definite + throw new \Exception(); } - $X = $B->toArray(); - $nx = $B->getN(); + $X = $B->toArray(); + $n = $B->getN(); - for ($k = 0; $k < $this->m; ++$k) { - for ($i = $k + 1; $i < $this->m; ++$i) { - for ($j = 0; $j < $nx; ++$j) { - $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k]; + // Solve L*Y = B; + for ($k = 0; $k < $this->m; $k++) { + for ($j = 0; $j < $n; $j++) { + for ($i = 0; $i < $k ; $i++) { + $X[$k][$j] -= $X[$i][$j] * $this->L[$k][$i]; } - } - for ($j = 0; $j < $nx; ++$j) { - $X[$k][$j] /= $this->L[$k][$k]; - } - } - - for ($k = $this->m - 1; $k >= 0; --$k) { - for ($j = 0; $j < $nx; ++$j) { - $X[$k][$j] /= $this->L[$k][$k]; - } - for ($i = 0; $i < $k; ++$i) { - for ($j = 0; $j < $nx; ++$j) { - $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i]; + + $X[$k][$j] /= $this->L[$k][$k]; + } + } + + // Solve L'*X = Y; + for ($k = $this->m - 1; $k >= 0; $k--) { + for ($j = 0; $j < $n; $j++) { + for ($i = $k + 1; $i < $this->m ; $i++) { + $X[$k][$j] -= $X[$i][$j] * $this->L[$i][$k]; } - } - } + + $X[$k][$j] /= $this->L[$k][$k]; + } + } - return new Matrix($X, $this->m, $nx); + $solution = new Matrix(); + $solution->setMatrix($X); + + return $solution; } } \ No newline at end of file diff --git a/Math/Matrix/LUDecomposition.php b/Math/Matrix/LUDecomposition.php index 124681658..b3e329518 100644 --- a/Math/Matrix/LUDecomposition.php +++ b/Math/Matrix/LUDecomposition.php @@ -40,40 +40,41 @@ class LUDecomposition $LUrowi = $LUcolj = []; for ($j = 0; $j < $this->n; ++$j) { - // Make a copy of the j-th column to localize references. for ($i = 0; $i < $this->m; ++$i) { $LUcolj[$i] = &$this->LU[$i][$j]; } - // Apply previous transformations. + for ($i = 0; $i < $this->m; ++$i) { $LUrowi = $this->LU[$i]; - // Most of the time is spent in the following dot product. $kmax = min($i, $j); $s = 0.0; + for ($k = 0; $k < $kmax; ++$k) { $s += $LUrowi[$k] * $LUcolj[$k]; } $LUrowi[$j] = $LUcolj[$i] -= $s; } - // Find pivot and exchange if necessary. + $p = $j; for ($i = $j + 1; $i < $this->m; ++$i) { if (abs($LUcolj[$i]) > abs($LUcolj[$p])) { $p = $i; } } + if ($p != $j) { for ($k = 0; $k < $this->n; ++$k) { $t = $this->LU[$p][$k]; $this->LU[$p][$k] = $this->LU[$j][$k]; $this->LU[$j][$k] = $t; } + $k = $this->piv[$p]; $this->piv[$p] = $this->piv[$j]; $this->piv[$j] = $k; $this->pivSign = $this->pivSign * -1; } - // Compute multipliers. + if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { for ($i = $j + 1; $i < $this->m; ++$i) { $this->LU[$i][$j] /= $this->LU[$j][$j]; @@ -82,7 +83,7 @@ class LUDecomposition } } - public function getL() + public function getL() : Matrix { $L = [[]]; @@ -104,7 +105,7 @@ class LUDecomposition return $matrix; } - public function getU() + public function getU() : Matrix { $U = [[]]; @@ -150,21 +151,25 @@ class LUDecomposition return $d; } - public function solve(Matrix $B) + public function solve(Matrix $B) : Matrix { if ($B->getM() !== $this->m) { + throw new \Exception(); } if (!$this->isNonsingular()) { + throw new \Exception(); } - $nx = $B->getM(); - $X = $B->getMatrix($this->piv, 0, $nx - 1); + $n = $B->getN(); + $X = $B->getMatrix($this->piv, 0, $n - 1); + // todo: fix get extract + // Solve L*Y = B(piv,:) for ($k = 0; $k < $this->n; ++$k) { for ($i = $k + 1; $i < $this->n; ++$i) { - for ($j = 0; $j < $nx; ++$j) { + for ($j = 0; $j < $n; ++$j) { $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; } } @@ -172,16 +177,19 @@ class LUDecomposition // Solve U*X = Y; for ($k = $this->n - 1; $k >= 0; --$k) { - for ($j = 0; $j < $nx; ++$j) { + for ($j = 0; $j < $n; ++$j) { $X[$k][$j] /= $this->LU[$k][$k]; } for ($i = 0; $i < $k; ++$i) { - for ($j = 0; $j < $nx; ++$j) { + for ($j = 0; $j < $n; ++$j) { $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; } } } - return $X; + $solution = new Matrix(); + $solution->setMatrix($X); + + return $solution; } } \ No newline at end of file diff --git a/Math/Matrix/Vector.php b/Math/Matrix/Vector.php index 839079c39..bd0da5754 100644 --- a/Math/Matrix/Vector.php +++ b/Math/Matrix/Vector.php @@ -33,8 +33,10 @@ class Vector extends Matrix * * @since 1.0.0 */ - public function __construct(int $m) + public function __construct(int $m = 1) { parent::__construct($m); } + + // todo: maybe overwrite setMatrix since only one column (is only a visual improvement) } \ No newline at end of file