diff options
Diffstat (limited to 'vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php')
| -rw-r--r-- | vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php | 249 |
1 files changed, 249 insertions, 0 deletions
diff --git a/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php new file mode 100644 index 0000000..373c074 --- /dev/null +++ b/vendor/phpoffice/phpspreadsheet/src/PhpSpreadsheet/Shared/JAMA/QRDecomposition.php @@ -0,0 +1,249 @@ +<?php
+
+namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
+
+use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
+
+/**
+ * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
+ * orthogonal matrix Q and an n-by-n upper triangular matrix R so that
+ * A = Q*R.
+ *
+ * The QR decompostion always exists, even if the matrix does not have
+ * full rank, so the constructor will never fail. The primary use of the
+ * QR decomposition is in the least squares solution of nonsquare systems
+ * of simultaneous linear equations. This will fail if isFullRank()
+ * returns false.
+ *
+ * @author Paul Meagher
+ *
+ * @version 1.1
+ */
+class QRDecomposition
+{
+ const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';
+
+ /**
+ * Array for internal storage of decomposition.
+ *
+ * @var array
+ */
+ private $QR = [];
+
+ /**
+ * Row dimension.
+ *
+ * @var int
+ */
+ private $m;
+
+ /**
+ * Column dimension.
+ *
+ * @var int
+ */
+ private $n;
+
+ /**
+ * Array for internal storage of diagonal of R.
+ *
+ * @var array
+ */
+ private $Rdiag = [];
+
+ /**
+ * QR Decomposition computed by Householder reflections.
+ *
+ * @param matrix $A Rectangular matrix
+ */
+ public function __construct($A)
+ {
+ if ($A instanceof Matrix) {
+ // Initialize.
+ $this->QR = $A->getArray();
+ $this->m = $A->getRowDimension();
+ $this->n = $A->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;
+ }
+ } else {
+ throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
+ }
+ }
+
+ // function __construct()
+
+ /**
+ * Is the matrix full rank?
+ *
+ * @return bool true if R, and hence A, has full rank, else false
+ */
+ public function isFullRank()
+ {
+ for ($j = 0; $j < $this->n; ++$j) {
+ if ($this->Rdiag[$j] == 0) {
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+ // function isFullRank()
+
+ /**
+ * Return the Householder vectors.
+ *
+ * @return Matrix Lower trapezoidal matrix whose columns define the reflections
+ */
+ public function getH()
+ {
+ $H = [];
+ 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;
+ }
+ }
+ }
+
+ return new Matrix($H);
+ }
+
+ // function getH()
+
+ /**
+ * Return the upper triangular factor.
+ *
+ * @return Matrix upper triangular factor
+ */
+ public function getR()
+ {
+ $R = [];
+ 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;
+ }
+ }
+ }
+
+ return new Matrix($R);
+ }
+
+ // function getR()
+
+ /**
+ * Generate and return the (economy-sized) orthogonal factor.
+ *
+ * @return Matrix orthogonal factor
+ */
+ public function getQ()
+ {
+ $Q = [];
+ 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];
+ }
+ }
+ }
+ }
+
+ return new Matrix($Q);
+ }
+
+ // function getQ()
+
+ /**
+ * Least squares solution of A*X = B.
+ *
+ * @param Matrix $B a Matrix with as many rows as A and any number of columns
+ *
+ * @return Matrix matrix that minimizes the two norm of Q*R*X-B
+ */
+ public function solve($B)
+ {
+ if ($B->getRowDimension() == $this->m) {
+ if ($this->isFullRank()) {
+ // Copy right hand side
+ $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];
+ }
+ }
+ }
+ $X = new Matrix($X);
+
+ return $X->getMatrix(0, $this->n - 1, 0, $nx);
+ }
+
+ throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
+ }
+
+ throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
+ }
+}
|