<?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);
    }
}

• [F] Decomposition.php
  Delete
• [F] LU.php
  Delete
• [F] QR.php
  Delete

Blog Section

Discover the Latest in Digital Payments and NFC Technology

Dive into our blog to explore the cutting-edge trends in digital payments and NFC technology. Stay updated on the innovations that are revolutionizing transactions, boosting security, and making payments quicker and more convenient. Learn how these advancements are shaping the future of financial interactions and driving the global transition towards a cashless world.

The Rise of Contactless Payments:...

In recent years, contactless payments have surged in popularity, driven...

Visit Blog

March 27, 2024

Enhancing Payment Security: The Role...

As digital transactions proliferate, ensuring robust payment security is more critical than ever. Two foundational...

Visit Blog

March 27, 2024

The Future of Digital Wallets:...

Digital wallets have fundamentally transformed how we manage money, offering a streamlined, secure, and highly...

Visit Blog

March 27, 2024