MATRIX_RANK_EXCEPTION
MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.'
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.
$QR : array
Array for internal storage of decomposition.
$m : int
Row dimension.
$n : int
Column dimension.
$Rdiag : array
Array for internal storage of diagonal of R.
solve(\PhpOffice\PhpSpreadsheet\Shared\JAMA\Matrix $B) : \PhpOffice\PhpSpreadsheet\Shared\JAMA\Matrix
Least squares solution of A*X = B.
| \PhpOffice\PhpSpreadsheet\Shared\JAMA\Matrix | $B | a Matrix with as many rows as A and any number of columns |
matrix that minimizes the two norm of QRX-B
<?php
namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
/**
* 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
*
* @version 1.1
*/
class QRDecomposition
{
const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';
/**
* Array for internal storage of decomposition.
*
* @var array
*/
private $QR = [];
/**
* Row dimension.
*
* @var int
*/
private $m;
/**
* Column dimension.
*
* @var int
*/
private $n;
/**
* Array for internal storage of diagonal of R.
*
* @var array
*/
private $Rdiag = [];
/**
* QR Decomposition computed by Householder reflections.
*
* @param matrix $A Rectangular matrix
*/
public function __construct($A)
{
if ($A instanceof Matrix) {
// Initialize.
$this->QR = $A->getArray();
$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 CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
}
}
// function __construct()
/**
* Is the matrix full rank?
*
* @return bool 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()
{
$H = [];
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()
{
$R = [];
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()
{
$Q = [];
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];
}
}
}
}
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(Matrix $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);
}
throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
}
throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
}
}