1 files changed, 191 insertions, 0 deletions
diff --git a/vendor/markbaker/matrix/classes/src/Decomposition/QR.php b/vendor/markbaker/matrix/classes/src/Decomposition/QR.php
new file mode 100644
index 0000000..93bf414
--- /dev/null
+++ b/vendor/markbaker/matrix/classes/src/Decomposition/QR.php
@@ -0,0 +1,191 @@
+<?php
+
+namespace Matrix\Decomposition;
+
+use Matrix\Exception;
+use Matrix\Matrix;
+
+class QR
+{
+    private $qrMatrix;
+    private $rows;
+    private $columns;
+
+    private $rDiagonal = [];
+
+    public function __construct(Matrix $matrix)
+    {
+        $this->qrMatrix = $matrix->toArray();
+        $this->rows = $matrix->rows;
+        $this->columns = $matrix->columns;
+
+        $this->decompose();
+    }
+
+    public function getHouseholdVectors(): Matrix
+    {
+        $householdVectors = [];
+        for ($row = 0; $row < $this->rows; ++$row) {
+            for ($column = 0; $column < $this->columns; ++$column) {
+                if ($row >= $column) {
+                    $householdVectors[$row][$column] = $this->qrMatrix[$row][$column];
+                } else {
+                    $householdVectors[$row][$column] = 0.0;
+                }
+            }
+        }
+
+        return new Matrix($householdVectors);
+    }
+
+    public function getQ(): Matrix
+    {
+        $qGrid = [];
+
+        $rowCount = $this->rows;
+        for ($k = $this->columns - 1; $k >= 0; --$k) {
+            for ($i = 0; $i < $this->rows; ++$i) {
+                $qGrid[$i][$k] = 0.0;
+            }
+            $qGrid[$k][$k] = 1.0;
+            if ($this->columns > $this->rows) {
+                $qGrid = array_slice($qGrid, 0, $this->rows);
+            }
+
+            for ($j = $k; $j < $this->columns; ++$j) {
+                if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) {
+                    $s = 0.0;
+                    for ($i = $k; $i < $this->rows; ++$i) {
+                        $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j];
+                    }
+                    $s = -$s / $this->qrMatrix[$k][$k];
+                    for ($i = $k; $i < $this->rows; ++$i) {
+                        $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k];
+                    }
+                }
+            }
+        }
+
+        array_walk(
+            $qGrid,
+            function (&$row) use ($rowCount) {
+                $row = array_reverse($row);
+                $row = array_slice($row, 0, $rowCount);
+            }
+        );
+
+        return new Matrix($qGrid);
+    }
+
+    public function getR(): Matrix
+    {
+        $rGrid = [];
+
+        for ($row = 0; $row < $this->columns; ++$row) {
+            for ($column = 0; $column < $this->columns; ++$column) {
+                if ($row < $column) {
+                    $rGrid[$row][$column] = $this->qrMatrix[$row][$column] ?? 0.0;
+                } elseif ($row === $column) {
+                    $rGrid[$row][$column] = $this->rDiagonal[$row] ?? 0.0;
+                } else {
+                    $rGrid[$row][$column] = 0.0;
+                }
+            }
+        }
+
+        if ($this->columns > $this->rows) {
+            $rGrid = array_slice($rGrid, 0, $this->rows);
+        }
+
+        return new Matrix($rGrid);
+    }
+
+    private function hypo($a, $b): float
+    {
+        if (abs($a) > abs($b)) {
+            $r = $b / $a;
+            $r = abs($a) * sqrt(1 + $r * $r);
+        } elseif ($b != 0.0) {
+            $r = $a / $b;
+            $r = abs($b) * sqrt(1 + $r * $r);
+        } else {
+            $r = 0.0;
+        }
+
+        return $r;
+    }
+
+    /**
+     * QR Decomposition computed by Householder reflections.
+     */
+    private function decompose(): void
+    {
+        for ($k = 0; $k < $this->columns; ++$k) {
+            // Compute 2-norm of k-th column without under/overflow.
+            $norm = 0.0;
+            for ($i = $k; $i < $this->rows; ++$i) {
+                $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]);
+            }
+            if ($norm != 0.0) {
+                // Form k-th Householder vector.
+                if ($this->qrMatrix[$k][$k] < 0.0) {
+                    $norm = -$norm;
+                }
+                for ($i = $k; $i < $this->rows; ++$i) {
+                    $this->qrMatrix[$i][$k] /= $norm;
+                }
+                $this->qrMatrix[$k][$k] += 1.0;
+                // Apply transformation to remaining columns.
+                for ($j = $k + 1; $j < $this->columns; ++$j) {
+                    $s = 0.0;
+                    for ($i = $k; $i < $this->rows; ++$i) {
+                        $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j];
+                    }
+                    $s = -$s / $this->qrMatrix[$k][$k];
+                    for ($i = $k; $i < $this->rows; ++$i) {
+                        $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k];
+                    }
+                }
+            }
+            $this->rDiagonal[$k] = -$norm;
+        }
+    }
+
+    public function isFullRank(): bool
+    {
+        for ($j = 0; $j < $this->columns; ++$j) {
+            if ($this->rDiagonal[$j] == 0.0) {
+                return false;
+            }
+        }
+
+        return true;
+    }
+
+    /**
+     * Least squares solution of A*X = B.
+     *
+     * @param Matrix $B a Matrix with as many rows as A and any number of columns
+     *
+     * @throws Exception
+     *
+     * @return Matrix matrix that minimizes the two norm of Q*R*X-B
+     */
+    public function solve(Matrix $B): Matrix
+    {
+        if ($B->rows !== $this->rows) {
+            throw new Exception('Matrix row dimensions are not equal');
+        }
+
+        if (!$this->isFullRank()) {
+            throw new Exception('Can only perform this operation on a full-rank matrix');
+        }
+
+        // Compute Y = transpose(Q)*B
+        $Y = $this->getQ()->transpose()
+            ->multiply($B);
+        // Solve R*X = Y;
+        return $this->getR()->inverse()
+            ->multiply($Y);
+    }
+}