5 * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
6 * orthogonal matrix Q and an n-by-n upper triangular matrix R so that
9 * The QR decompostion always exists, even if the matrix does not have
10 * full rank, so the constructor will never fail. The primary use of the
11 * QR decomposition is in the least squares solution of nonsquare systems
12 * of simultaneous linear equations. This will fail if isFullRank()
15 * @author Paul Meagher
19 class PHPExcel_Shared_JAMA_QRDecomposition
{
21 const MatrixRankException
= "Can only perform operation on full-rank matrix.";
24 * Array for internal storage of decomposition.
27 private $QR = array();
42 * Array for internal storage of diagonal of R.
45 private $Rdiag = array();
49 * QR Decomposition computed by Householder reflections.
51 * @param matrix $A Rectangular matrix
52 * @return Structure to access R and the Householder vectors and compute Q.
54 public function __construct($A) {
55 if($A instanceof PHPExcel_Shared_JAMA_Matrix
) {
57 $this->QR
= $A->getArrayCopy();
58 $this->m
= $A->getRowDimension();
59 $this->n
= $A->getColumnDimension();
61 for ($k = 0; $k < $this->n
; ++
$k) {
62 // Compute 2-norm of k-th column without under/overflow.
64 for ($i = $k; $i < $this->m
; ++
$i) {
65 $nrm = hypo($nrm, $this->QR
[$i][$k]);
68 // Form k-th Householder vector.
69 if ($this->QR
[$k][$k] < 0) {
72 for ($i = $k; $i < $this->m
; ++
$i) {
73 $this->QR
[$i][$k] /= $nrm;
75 $this->QR
[$k][$k] +
= 1.0;
76 // Apply transformation to remaining columns.
77 for ($j = $k+
1; $j < $this->n
; ++
$j) {
79 for ($i = $k; $i < $this->m
; ++
$i) {
80 $s +
= $this->QR
[$i][$k] * $this->QR
[$i][$j];
82 $s = -$s/$this->QR
[$k][$k];
83 for ($i = $k; $i < $this->m
; ++
$i) {
84 $this->QR
[$i][$j] +
= $s * $this->QR
[$i][$k];
88 $this->Rdiag
[$k] = -$nrm;
91 throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix
::ArgumentTypeException
);
93 } // function __construct()
97 * Is the matrix full rank?
99 * @return boolean true if R, and hence A, has full rank, else false.
101 public function isFullRank() {
102 for ($j = 0; $j < $this->n
; ++
$j) {
103 if ($this->Rdiag
[$j] == 0) {
108 } // function isFullRank()
112 * Return the Householder vectors
114 * @return Matrix Lower trapezoidal matrix whose columns define the reflections
116 public function getH() {
117 for ($i = 0; $i < $this->m
; ++
$i) {
118 for ($j = 0; $j < $this->n
; ++
$j) {
120 $H[$i][$j] = $this->QR
[$i][$j];
126 return new PHPExcel_Shared_JAMA_Matrix($H);
131 * Return the upper triangular factor
133 * @return Matrix upper triangular factor
135 public function getR() {
136 for ($i = 0; $i < $this->n
; ++
$i) {
137 for ($j = 0; $j < $this->n
; ++
$j) {
139 $R[$i][$j] = $this->QR
[$i][$j];
140 } elseif ($i == $j) {
141 $R[$i][$j] = $this->Rdiag
[$i];
147 return new PHPExcel_Shared_JAMA_Matrix($R);
152 * Generate and return the (economy-sized) orthogonal factor
154 * @return Matrix orthogonal factor
156 public function getQ() {
157 for ($k = $this->n
-1; $k >= 0; --$k) {
158 for ($i = 0; $i < $this->m
; ++
$i) {
162 for ($j = $k; $j < $this->n
; ++
$j) {
163 if ($this->QR
[$k][$k] != 0) {
165 for ($i = $k; $i < $this->m
; ++
$i) {
166 $s +
= $this->QR
[$i][$k] * $Q[$i][$j];
168 $s = -$s/$this->QR
[$k][$k];
169 for ($i = $k; $i < $this->m
; ++
$i) {
170 $Q[$i][$j] +
= $s * $this->QR
[$i][$k];
176 for($i = 0; $i < count($Q); ++$i) {
177 for($j = 0; $j < count($Q); ++$j) {
178 if(! isset($Q[$i][$j]) ) {
184 return new PHPExcel_Shared_JAMA_Matrix($Q);
189 * Least squares solution of A*X = B
191 * @param Matrix $B A Matrix with as many rows as A and any number of columns.
192 * @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
194 public function solve($B) {
195 if ($B->getRowDimension() == $this->m
) {
196 if ($this->isFullRank()) {
197 // Copy right hand side
198 $nx = $B->getColumnDimension();
199 $X = $B->getArrayCopy();
200 // Compute Y = transpose(Q)*B
201 for ($k = 0; $k < $this->n
; ++
$k) {
202 for ($j = 0; $j < $nx; ++
$j) {
204 for ($i = $k; $i < $this->m
; ++
$i) {
205 $s +
= $this->QR
[$i][$k] * $X[$i][$j];
207 $s = -$s/$this->QR
[$k][$k];
208 for ($i = $k; $i < $this->m
; ++
$i) {
209 $X[$i][$j] +
= $s * $this->QR
[$i][$k];
214 for ($k = $this->n
-1; $k >= 0; --$k) {
215 for ($j = 0; $j < $nx; ++
$j) {
216 $X[$k][$j] /= $this->Rdiag
[$k];
218 for ($i = 0; $i < $k; ++
$i) {
219 for ($j = 0; $j < $nx; ++
$j) {
220 $X[$i][$j] -= $X[$k][$j]* $this->QR
[$i][$k];
224 $X = new PHPExcel_Shared_JAMA_Matrix($X);
225 return ($X->getMatrix(0, $this->n
-1, 0, $nx));
227 throw new PHPExcel_Calculation_Exception(self
::MatrixRankException
);
230 throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix
::MatrixDimensionException
);
232 } // function solve()
234 } // PHPExcel_Shared_JAMA_class QRDecomposition