You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
235 lines
7.5 KiB
235 lines
7.5 KiB
<?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 MATRIX_RANK_EXCEPTION = "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::ARGUMENT_TYPE_EXCEPTION);
|
|
}
|
|
} // 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::MATRIX_RANK_EXCEPTION);
|
|
}
|
|
} else {
|
|
throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MATRIX_DIMENSION_EXCEPTION);
|
|
}
|
|
}
|
|
}
|