diff --git a/Math/Matrix/Cholesky.php b/Math/Matrix/Cholesky.php deleted file mode 100644 index 3d4237aad..000000000 --- a/Math/Matrix/Cholesky.php +++ /dev/null @@ -1,21 +0,0 @@ -L = $M->toArray(); + $this->m = $M->getM(); + + for($i = 0; $i < $this->m; ++$i) { + for($j = $i; $j < $this->m; ++$j) { + 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 ($sum >= 0) { + $this->L[$i][$i] = sqrt($sum); + } else { + $this->isSpd = false; + } + } else { + if ($this->L[$i][$i] != 0) { + $this->L[$j][$i] = $sum / $this->L[$i][$i]; + } + } + } + + for ($k = $i+1; $k < $this->m; ++$k) { + $this->L[$i][$k] = 0.0; + } + } + } + + public function isSpd() + { + return $this->isSpd; + } + + public function getL() + { + $matrix = new Matrix(); + $matrix->setMatrix($this->L); + + return $matrix; + } + + public function solve(Matrix $B) + { + if ($B->getM() !== $this->m) { + // invalid dimension + } + + if (!$this->isSpd) { + // is not positive definite + } + + $X = $B->toArray(); + $nx = $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]; + } + } + 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]; + } + } + } + + return new Matrix($X, $this->m, $nx); + } +} \ No newline at end of file diff --git a/Math/Matrix/EigenvalueDecomposition.php b/Math/Matrix/EigenvalueDecomposition.php new file mode 100644 index 000000000..e69de29bb diff --git a/Math/Matrix/LUDecomposition.php b/Math/Matrix/LUDecomposition.php new file mode 100644 index 000000000..6517f2c4d --- /dev/null +++ b/Math/Matrix/LUDecomposition.php @@ -0,0 +1,185 @@ +LU = $M->toArray(); + $this->m = $M->getM(); + $this->n = $M->getN(); + + for ($i = 0; $i < $this->m; ++$i) { + $this->piv[$i] = $i; + } + + $this->pivSign = 1; + $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]; + } + } + } + } + + public function getL() + { + for ($i = 0; $i < $this->m; ++$i) { + for ($j = 0; $j < $this->n; ++$j) { + if ($i > $j) { + $L[$i][$j] = $this->LU[$i][$j]; + } elseif ($i == $j) { + $L[$i][$j] = 1.0; + } else { + $L[$i][$j] = 0.0; + } + } + } + + $matrix = new Matrix(); + $matrix->setMatrix($L); + + return $matrix; + } + + public function getU() + { + for ($i = 0; $i < $this->n; ++$i) { + for ($j = 0; $j < $this->n; ++$j) { + if ($i <= $j) { + $U[$i][$j] = $this->LU[$i][$j]; + } else { + $U[$i][$j] = 0.0; + } + } + } + + $matrix = new Matrix(); + $matrix->setMatrix($U); + + return $matrix; + } + + public function getPivot() + { + return $this->piv; + } + + public function isNonsingular() : bool + { + for ($j = 0; $j < $this->n; ++$j) { + if ($this->LU[$j][$j] == 0) { + return false; + } + } + + return true; + } + + public function det() + { + $d = $this->pivSign; + for ($j = 0; $j < $this->n; ++$j) { + $d *= $this->LU[$j][$j]; + } + + return $d; + } + + public function solve(Matrix $B) + { + if ($B->getM() !== $this->m) { + } + + if (!$this->isNonsingular()) { + } + + var_dump($this->piv); + + $nx = $B->getM(); + $X = $B->getMatrix($this->piv, 0, $nx-1); + + // 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) { + $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; + } + } + } + + // Solve U*X = Y; + for ($k = $this->n-1; $k >= 0; --$k) { + for ($j = 0; $j < $nx; ++$j) { + $X->A[$k][$j] /= $this->LU[$k][$k]; + } + for ($i = 0; $i < $k; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X->A[$i][$j] -= $X->A[$k][$j] * $this->LU[$i][$k]; + } + } + } + + return $X; + } +} \ No newline at end of file diff --git a/Math/Matrix/Matrix.php b/Math/Matrix/Matrix.php index c040dc0c9..d802dc60a 100644 --- a/Math/Matrix/Matrix.php +++ b/Math/Matrix/Matrix.php @@ -68,7 +68,7 @@ class Matrix implements \ArrayAccess, \Iterator * * @since 1.0.0 */ - public function __construct(int $m, int $n = 1) + public function __construct(int $m = 1, int $n = 1) { $this->n = $n; $this->m = $m; @@ -146,6 +146,18 @@ class Matrix implements \ArrayAccess, \Iterator return $this->matrix; } + /** + * Get matrix array. + * + * @return array + * + * @since 1.0.0 + */ + public function toArray() : array + { + return $this->matrix; + } + /** * Get matrix rank. * @@ -171,10 +183,8 @@ class Matrix implements \ArrayAccess, \Iterator */ public function setMatrix(array $matrix) : Matrix { - if ($this->m !== count($matrix) || $this->n !== count($matrix[0])) { - throw new InvalidDimensionException(count($matrix) . 'x' . count($matrix[0])); - } - + $this->m = count($matrix); + $this->n = count($matrix[0] ?? 1); $this->matrix = $matrix; return $this; @@ -491,16 +501,7 @@ class Matrix implements \ArrayAccess, \Iterator */ public function inverse(int $algorithm = InverseType::GAUSS_JORDAN) : Matrix { - if ($this->n !== $this->m) { - throw new InvalidDimensionException($this->m . 'x' . $this->n); - } - - switch ($algorithm) { - case InverseType::GAUSS_JORDAN: - return $this->inverseGaussJordan(); - default: - throw new \Exception('Inversion algorithm'); - } + return $this->solve(new IdentityMatrix($this->m, $this->m)); } /** @@ -561,9 +562,11 @@ class Matrix implements \ArrayAccess, \Iterator return $newMatrix; } - public function solve($b, int $algorithm) : Matrix + public function solve(Matrix $B) : Matrix { - return $this->gaussElimination($b); + $M = $this->m === $this->n ? new LUDecomposition($this) : new QRDecomposition($this); + + return $M->solve($B); } private function gaussElimination($b) : Matrix @@ -671,18 +674,8 @@ class Matrix implements \ArrayAccess, \Iterator */ public function det() : float { - if ($this->n === 1) { - return $this->matrix[0][0]; - } - - $trianglize = $this->matrix; - $prod = $this->upperTrianglize($trianglize); - - for ($i = 0; $i < $this->n; $i++) { - $prod *= $trianglize[$i][$i]; - } - - return $prod; + $L = new LUDecomposition($this); + return $L->det(); } /** diff --git a/Math/Matrix/QRDecomposition.php b/Math/Matrix/QRDecomposition.php new file mode 100644 index 000000000..fb735ad08 --- /dev/null +++ b/Math/Matrix/QRDecomposition.php @@ -0,0 +1,187 @@ +QR = $M->toArray(); + $this->m = $M->getRowDimension(); + $this->n = $M->getColumnDimension(); + + // Main loop. + for ($k = 0; $k < $this->n; ++$k) { + // Compute 2-norm of k-th column without under/overflow. + $nrm = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $nrm = hypo($nrm, $this->QR[$i][$k]); + } + + if ($nrm != 0.0) { + // Form k-th Householder vector. + if ($this->QR[$k][$k] < 0) { + $nrm = -$nrm; + } + + for ($i = $k; $i < $this->m; ++$i) { + $this->QR[$i][$k] /= $nrm; + } + + $this->QR[$k][$k] += 1.0; + // Apply transformation to remaining columns. + for ($j = $k+1; $j < $this->n; ++$j) { + $s = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $s += $this->QR[$i][$k] * $this->QR[$i][$j]; + } + + $s = -$s/$this->QR[$k][$k]; + for ($i = $k; $i < $this->m; ++$i) { + $this->QR[$i][$j] += $s * $this->QR[$i][$k]; + } + } + } + + $this->Rdiag[$k] = -$nrm; + } + } + + public function isFullRank() : bool + { + for ($j = 0; $j < $this->n; ++$j) { + if ($this->Rdiag[$j] == 0) { + return false; + } + } + + return true; + } + + public function getH() + { + for ($i = 0; $i < $this->m; ++$i) { + for ($j = 0; $j < $this->n; ++$j) { + if ($i >= $j) { + $H[$i][$j] = $this->QR[$i][$j]; + } else { + $H[$i][$j] = 0.0; + } + } + } + + $matrix = new Matrix(); + $matrix->setArray($H); + + return $this->matrix; + } + + public function getR() + { + for ($i = 0; $i < $this->n; ++$i) { + for ($j = 0; $j < $this->n; ++$j) { + if ($i < $j) { + $R[$i][$j] = $this->QR[$i][$j]; + } elseif ($i == $j) { + $R[$i][$j] = $this->Rdiag[$i]; + } else { + $R[$i][$j] = 0.0; + } + } + } + + $matrix = new Matrix(); + $matrix->setArray($R); + + return $this->matrix; + } + + public function getQ() + { + for ($k = $this->n-1; $k >= 0; --$k) { + for ($i = 0; $i < $this->m; ++$i) { + $Q[$i][$k] = 0.0; + } + $Q[$k][$k] = 1.0; + for ($j = $k; $j < $this->n; ++$j) { + if ($this->QR[$k][$k] != 0) { + $s = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $s += $this->QR[$i][$k] * $Q[$i][$j]; + } + $s = -$s/$this->QR[$k][$k]; + for ($i = $k; $i < $this->m; ++$i) { + $Q[$i][$j] += $s * $this->QR[$i][$k]; + } + } + } + } + + $matrix = new Matrix(); + $matrix->setArray($Q); + + return $this->matrix; + } + + public function solve(Matrix $B) + { + if ($B->getRowDimension() !== $this->m) { + } + + if (!$this->isFullRank()) { + } + + $nx = $B->getColumnDimension(); + $X = $B->getArrayCopy(); + // Compute Y = transpose(Q)*B + for ($k = 0; $k < $this->n; ++$k) { + for ($j = 0; $j < $nx; ++$j) { + $s = 0.0; + for ($i = $k; $i < $this->m; ++$i) { + $s += $this->QR[$i][$k] * $X[$i][$j]; + } + $s = -$s/$this->QR[$k][$k]; + for ($i = $k; $i < $this->m; ++$i) { + $X[$i][$j] += $s * $this->QR[$i][$k]; + } + } + } + // Solve R*X = Y; + for ($k = $this->n-1; $k >= 0; --$k) { + for ($j = 0; $j < $nx; ++$j) { + $X[$k][$j] /= $this->Rdiag[$k]; + } + for ($i = 0; $i < $k; ++$i) { + for ($j = 0; $j < $nx; ++$j) { + $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k]; + } + } + } + + $matrix = new Matrix(); + $matrix->setArray($X); + + return $matrix->getMatrix(0, $this->n-1, 0, $nx); + } +} \ No newline at end of file diff --git a/Math/Matrix/SingularValueDecomposition.php b/Math/Matrix/SingularValueDecomposition.php new file mode 100644 index 000000000..e69de29bb diff --git a/Math/Statistic/Forecast/Regression/MultipleLinearRegression.php b/Math/Statistic/Forecast/Regression/MultipleLinearRegression.php index 56fcc32dc..ad4781be4 100644 --- a/Math/Statistic/Forecast/Regression/MultipleLinearRegression.php +++ b/Math/Statistic/Forecast/Regression/MultipleLinearRegression.php @@ -30,7 +30,6 @@ class MultipleLinearRegression $Y = new Matrix(count($y)); $Y->setMatrix($y); - return $XT->mult($X)->inverse()->mult($XT)->mult($Y)->getMatrix(); }