| New file |
|---|
| 0,0 → 1,234 |
|---|
| <?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 PHPExcel_Shared_JAMA_QRDecomposition { |
| |
| const MatrixRankException = "Can only perform operation on full-rank matrix."; |
| |
| /** |
| * 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 PHPExcel_Shared_JAMA_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 PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::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 PHPExcel_Shared_JAMA_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 PHPExcel_Shared_JAMA_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 PHPExcel_Shared_JAMA_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 PHPExcel_Shared_JAMA_Matrix($X); |
| return ($X->getMatrix(0, $this->n-1, 0, $nx)); |
| } else { |
| throw new PHPExcel_Calculation_Exception(self::MatrixRankException); |
| } |
| } else { |
| throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException); |
| } |
| } // function solve() |
| |
| } // PHPExcel_Shared_JAMA_class QRDecomposition |