QRDecomposition.php

Go to the documentation of this file.

<?php
class QRDecomposition {

        private $QR = array();

        private $m;

        private $n;

        private $Rdiag = array();


        public function __construct($A) {
                if($A instanceof Matrix) {
                        // Initialize.
                        $this->QR = $A->getArrayCopy();
                        $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 Exception(JAMAError(ArgumentTypeException));
                }
        }       //      function __construct()


        public function isFullRank() {
                for ($j = 0; $j < $this->n; ++$j) {
                        if ($this->Rdiag[$j] == 0) {
                                return false;
                        }
                }
                return true;
        }       //      function isFullRank()


        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;
                                }
                        }
                }
                return new Matrix($H);
        }       //      function getH()


        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;
                                }
                        }
                }
                return new Matrix($R);
        }       //      function getR()


        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];
                                        }
                                }
                        }
                }
                /*
                for($i = 0; $i < count($Q); ++$i) {
                        for($j = 0; $j < count($Q); ++$j) {
                                if(! isset($Q[$i][$j]) ) {
                                        $Q[$i][$j] = 0;
                                }
                        }
                }
                */
                return new Matrix($Q);
        }       //      function getQ()


        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));
                        } else {
                                throw new Exception(JAMAError(MatrixRankException));
                        }
                } else {
                        throw new Exception(JAMAError(MatrixDimensionException));
                }
        }       //      function solve()

}       //      class QRDecomposition