JAMA

# Source for file QRDecomposition.php Documentation is available at QRDecomposition.php <?php /** * @package JAMA * * 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 * @license PHP v3.0 * @version 1.1 */ class QRDecomposition { /** * Array for internal storage of decomposition. * @var array */ private $QR = array(); /** * Row dimension. * @var integer */ private $m; /** * Column dimension. * @var integer */ private $n; /** * Array for internal storage of diagonal of R. * @var array */ private $Rdiag = array(); /** * QR Decomposition computed by Householder reflections. * * @param matrix $A Rectangular matrix * @return Structure to access R and the Householder vectors and compute Q. */ 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() /** * Is the matrix full rank? * * @return boolean 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() { 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() { 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() { 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() /** * 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)); } else { throw new Exception(JAMAError(MatrixRankException)); } } else { throw new Exception(JAMAError(MatrixDimensionException)); } } // function solve() } // class QRDecomposition Documentation generated on Mon, 11 Jan 2010 08:13:08 +0100 by phpDocumentor 1.4.1