Jingga/phpOMS
1
0
Fork 0
mirror of https://github.com/Karaka-Management/phpOMS.git synced 2026-02-10 05:58:42 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Improve matrix implementations

Browse Source
This commit is contained in:
Dennis Eichhorn 2017-11-15 13:04:43 +01:00
parent 102075cd91
commit 39bb8c83d4
3 changed files with 62 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

66
Math/Matrix/CholeskyDecomposition.php
View File

@ -21,8 +21,12 @@ class CholeskyDecomposition
private $m = 0; private $m = 0;
/**
* Is symmetric positive definite
*/
private $isSpd = true; private $isSpd = true;
// see http://www.aip.de/groups/soe/local/numres/bookfpdf/f2-9.pdf
public function __construct(Matrix $M) public function __construct(Matrix $M)
{ {
$this->L = $M->toArray(); $this->L = $M->toArray();
@ -33,14 +37,15 @@ class CholeskyDecomposition
for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) { for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
$sum -= $this->L[$i][$k] * $this->L[$j][$k]; $sum -= $this->L[$i][$k] * $this->L[$j][$k];
} }
if ($i == $j) {
if ($i === $j) {
if ($sum >= 0) { if ($sum >= 0) {
$this->L[$i][$i] = sqrt($sum); $this->L[$i][$i] = sqrt($sum);
} else { } else {
$this->isSpd = false; $this->isSpd = false;
} }
} else { } else {
if ($this->L[$i][$i] != 0) { if ($this->L[$i][$i] !== 0) {
$this->L[$j][$i] = $sum / $this->L[$i][$i]; $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; return $this->isSpd;
} }
public function getL() public function getL() : Matrix
{ {
$matrix = new Matrix(); $matrix = new Matrix();
$matrix->setMatrix($this->L); $matrix->setMatrix($this->L);
@ -65,41 +70,44 @@ class CholeskyDecomposition
return $matrix; return $matrix;
} }
public function solve(Matrix $B) public function solve(Matrix $B) : Matrix
{ {
if ($B->getM() !== $this->m) { if ($B->getM() !== $this->m) {
// invalid dimension throw new \Exception();
} }
if (!$this->isSpd) { if (!$this->isSpd) {
// is not positive definite throw new \Exception();
} }
$X = $B->toArray(); $X = $B->toArray();
$nx = $B->getN(); $n = $B->getN();
for ($k = 0; $k < $this->m; ++$k) { // Solve L*Y = B;
for ($i = $k + 1; $i < $this->m; ++$i) { for ($k = 0; $k < $this->m; $k++) {
for ($j = 0; $j < $nx; ++$j) { for ($j = 0; $j < $n; $j++) {
$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k]; 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];
$X[$k][$j] /= $this->L[$k][$k]; }
} }
}
// Solve L'*X = Y;
for ($k = $this->m - 1; $k >= 0; --$k) { for ($k = $this->m - 1; $k >= 0; $k--) {
for ($j = 0; $j < $nx; ++$j) { for ($j = 0; $j < $n; $j++) {
$X[$k][$j] /= $this->L[$k][$k]; for ($i = $k + 1; $i < $this->m ; $i++) {
} $X[$k][$j] -= $X[$i][$j] * $this->L[$i][$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];
}
}
return new Matrix($X, $this->m, $nx); $solution = new Matrix();
$solution->setMatrix($X);
return $solution;
} }
} }

36
Math/Matrix/LUDecomposition.php
View File

@ -40,40 +40,41 @@ class LUDecomposition
$LUrowi = $LUcolj = []; $LUrowi = $LUcolj = [];
for ($j = 0; $j < $this->n; ++$j) { 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) { for ($i = 0; $i < $this->m; ++$i) {
$LUcolj[$i] = &$this->LU[$i][$j]; $LUcolj[$i] = &$this->LU[$i][$j];
} }
// Apply previous transformations.
for ($i = 0; $i < $this->m; ++$i) { for ($i = 0; $i < $this->m; ++$i) {
$LUrowi = $this->LU[$i]; $LUrowi = $this->LU[$i];
// Most of the time is spent in the following dot product.
$kmax = min($i, $j); $kmax = min($i, $j);
$s = 0.0; $s = 0.0;
for ($k = 0; $k < $kmax; ++$k) { for ($k = 0; $k < $kmax; ++$k) {
$s += $LUrowi[$k] * $LUcolj[$k]; $s += $LUrowi[$k] * $LUcolj[$k];
} }
$LUrowi[$j] = $LUcolj[$i] -= $s; $LUrowi[$j] = $LUcolj[$i] -= $s;
} }
// Find pivot and exchange if necessary.
$p = $j; $p = $j;
for ($i = $j + 1; $i < $this->m; ++$i) { for ($i = $j + 1; $i < $this->m; ++$i) {
if (abs($LUcolj[$i]) > abs($LUcolj[$p])) { if (abs($LUcolj[$i]) > abs($LUcolj[$p])) {
$p = $i; $p = $i;
} }
} }
if ($p != $j) { if ($p != $j) {
for ($k = 0; $k < $this->n; ++$k) { for ($k = 0; $k < $this->n; ++$k) {
$t = $this->LU[$p][$k]; $t = $this->LU[$p][$k];
$this->LU[$p][$k] = $this->LU[$j][$k]; $this->LU[$p][$k] = $this->LU[$j][$k];
$this->LU[$j][$k] = $t; $this->LU[$j][$k] = $t;
} }
$k = $this->piv[$p]; $k = $this->piv[$p];
$this->piv[$p] = $this->piv[$j]; $this->piv[$p] = $this->piv[$j];
$this->piv[$j] = $k; $this->piv[$j] = $k;
$this->pivSign = $this->pivSign * -1; $this->pivSign = $this->pivSign * -1;
} }
// Compute multipliers.
if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
for ($i = $j + 1; $i < $this->m; ++$i) { for ($i = $j + 1; $i < $this->m; ++$i) {
$this->LU[$i][$j] /= $this->LU[$j][$j]; $this->LU[$i][$j] /= $this->LU[$j][$j];
@ -82,7 +83,7 @@ class LUDecomposition
} }
} }
public function getL() public function getL() : Matrix
{ {
$L = [[]]; $L = [[]];
@ -104,7 +105,7 @@ class LUDecomposition
return $matrix; return $matrix;
} }
public function getU() public function getU() : Matrix
{ {
$U = [[]]; $U = [[]];
@ -150,21 +151,25 @@ class LUDecomposition
return $d; return $d;
} }
public function solve(Matrix $B) public function solve(Matrix $B) : Matrix
{ {
if ($B->getM() !== $this->m) { if ($B->getM() !== $this->m) {
throw new \Exception();
} }
if (!$this->isNonsingular()) { if (!$this->isNonsingular()) {
throw new \Exception();
} }
$nx = $B->getM(); $n = $B->getN();
$X = $B->getMatrix($this->piv, 0, $nx - 1); $X = $B->getMatrix($this->piv, 0, $n - 1);
// todo: fix get extract
// Solve L*Y = B(piv,:) // Solve L*Y = B(piv,:)
for ($k = 0; $k < $this->n; ++$k) { for ($k = 0; $k < $this->n; ++$k) {
for ($i = $k + 1; $i < $this->n; ++$i) { 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]; $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
} }
} }
@ -172,16 +177,19 @@ class LUDecomposition
// Solve U*X = Y; // Solve U*X = Y;
for ($k = $this->n - 1; $k >= 0; --$k) { 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]; $X[$k][$j] /= $this->LU[$k][$k];
} }
for ($i = 0; $i < $k; ++$i) { 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]; $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
} }
} }
} }
return $X; $solution = new Matrix();
$solution->setMatrix($X);
return $solution;
} }
} }

4
Math/Matrix/Vector.php
View File

@ -33,8 +33,10 @@ class Vector extends Matrix
* *
* @since 1.0.0 * @since 1.0.0
*/ */
public function __construct(int $m) public function __construct(int $m = 1)
{ {
parent::__construct($m); parent::__construct($m);
} }
// todo: maybe overwrite setMatrix since only one column (is only a visual improvement)
} }