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 |