PageRenderTime 33ms CodeModel.GetById 19ms app.highlight 11ms RepoModel.GetById 1ms app.codeStats 0ms

/PHPExcel/Shared/JAMA/QRDecomposition.php

https://bitbucket.org/nfredricks/wp-employee-time
PHP | 234 lines | 134 code | 18 blank | 82 comment | 43 complexity | 95f548ec797ed8f8c06d6d134d33ebbd MD5 | raw file
  1<?php
  2/**
  3 *	@package JAMA
  4 *
  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
  7 *	A = Q*R.
  8 *
  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()
 13 *	returns false.
 14 *
 15 *	@author  Paul Meagher
 16 *	@license PHP v3.0
 17 *	@version 1.1
 18 */
 19class PHPExcel_Shared_JAMA_QRDecomposition {
 20
 21	const MatrixRankException	= "Can only perform operation on full-rank matrix.";
 22
 23	/**
 24	 *	Array for internal storage of decomposition.
 25	 *	@var array
 26	 */
 27	private $QR = array();
 28
 29	/**
 30	 *	Row dimension.
 31	 *	@var integer
 32	 */
 33	private $m;
 34
 35	/**
 36	*	Column dimension.
 37	*	@var integer
 38	*/
 39	private $n;
 40
 41	/**
 42	 *	Array for internal storage of diagonal of R.
 43	 *	@var  array
 44	 */
 45	private $Rdiag = array();
 46
 47
 48	/**
 49	 *	QR Decomposition computed by Householder reflections.
 50	 *
 51	 *	@param matrix $A Rectangular matrix
 52	 *	@return Structure to access R and the Householder vectors and compute Q.
 53	 */
 54	public function __construct($A) {
 55		if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
 56			// Initialize.
 57			$this->QR = $A->getArrayCopy();
 58			$this->m  = $A->getRowDimension();
 59			$this->n  = $A->getColumnDimension();
 60			// Main loop.
 61			for ($k = 0; $k < $this->n; ++$k) {
 62				// Compute 2-norm of k-th column without under/overflow.
 63				$nrm = 0.0;
 64				for ($i = $k; $i < $this->m; ++$i) {
 65					$nrm = hypo($nrm, $this->QR[$i][$k]);
 66				}
 67				if ($nrm != 0.0) {
 68					// Form k-th Householder vector.
 69					if ($this->QR[$k][$k] < 0) {
 70						$nrm = -$nrm;
 71					}
 72					for ($i = $k; $i < $this->m; ++$i) {
 73						$this->QR[$i][$k] /= $nrm;
 74					}
 75					$this->QR[$k][$k] += 1.0;
 76					// Apply transformation to remaining columns.
 77					for ($j = $k+1; $j < $this->n; ++$j) {
 78						$s = 0.0;
 79						for ($i = $k; $i < $this->m; ++$i) {
 80							$s += $this->QR[$i][$k] * $this->QR[$i][$j];
 81						}
 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];
 85						}
 86					}
 87				}
 88				$this->Rdiag[$k] = -$nrm;
 89			}
 90		} else {
 91			throw new Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
 92		}
 93	}	//	function __construct()
 94
 95
 96	/**
 97	 *	Is the matrix full rank?
 98	 *
 99	 *	@return boolean true if R, and hence A, has full rank, else false.
100	 */
101	public function isFullRank() {
102		for ($j = 0; $j < $this->n; ++$j) {
103			if ($this->Rdiag[$j] == 0) {
104				return false;
105			}
106		}
107		return true;
108	}	//	function isFullRank()
109
110
111	/**
112	 *	Return the Householder vectors
113	 *
114	 *	@return Matrix Lower trapezoidal matrix whose columns define the reflections
115	 */
116	public function getH() {
117		for ($i = 0; $i < $this->m; ++$i) {
118			for ($j = 0; $j < $this->n; ++$j) {
119				if ($i >= $j) {
120					$H[$i][$j] = $this->QR[$i][$j];
121				} else {
122					$H[$i][$j] = 0.0;
123				}
124			}
125		}
126		return new PHPExcel_Shared_JAMA_Matrix($H);
127	}	//	function getH()
128
129
130	/**
131	 *	Return the upper triangular factor
132	 *
133	 *	@return Matrix upper triangular factor
134	 */
135	public function getR() {
136		for ($i = 0; $i < $this->n; ++$i) {
137			for ($j = 0; $j < $this->n; ++$j) {
138				if ($i < $j) {
139					$R[$i][$j] = $this->QR[$i][$j];
140				} elseif ($i == $j) {
141					$R[$i][$j] = $this->Rdiag[$i];
142				} else {
143					$R[$i][$j] = 0.0;
144				}
145			}
146		}
147		return new PHPExcel_Shared_JAMA_Matrix($R);
148	}	//	function getR()
149
150
151	/**
152	 *	Generate and return the (economy-sized) orthogonal factor
153	 *
154	 *	@return Matrix orthogonal factor
155	 */
156	public function getQ() {
157		for ($k = $this->n-1; $k >= 0; --$k) {
158			for ($i = 0; $i < $this->m; ++$i) {
159				$Q[$i][$k] = 0.0;
160			}
161			$Q[$k][$k] = 1.0;
162			for ($j = $k; $j < $this->n; ++$j) {
163				if ($this->QR[$k][$k] != 0) {
164					$s = 0.0;
165					for ($i = $k; $i < $this->m; ++$i) {
166						$s += $this->QR[$i][$k] * $Q[$i][$j];
167					}
168					$s = -$s/$this->QR[$k][$k];
169					for ($i = $k; $i < $this->m; ++$i) {
170						$Q[$i][$j] += $s * $this->QR[$i][$k];
171					}
172				}
173			}
174		}
175		/*
176		for($i = 0; $i < count($Q); ++$i) {
177			for($j = 0; $j < count($Q); ++$j) {
178				if(! isset($Q[$i][$j]) ) {
179					$Q[$i][$j] = 0;
180				}
181			}
182		}
183		*/
184		return new PHPExcel_Shared_JAMA_Matrix($Q);
185	}	//	function getQ()
186
187
188	/**
189	 *	Least squares solution of A*X = B
190	 *
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.
193	 */
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) {
203						$s = 0.0;
204						for ($i = $k; $i < $this->m; ++$i) {
205							$s += $this->QR[$i][$k] * $X[$i][$j];
206						}
207						$s = -$s/$this->QR[$k][$k];
208						for ($i = $k; $i < $this->m; ++$i) {
209							$X[$i][$j] += $s * $this->QR[$i][$k];
210						}
211					}
212				}
213				// Solve R*X = Y;
214				for ($k = $this->n-1; $k >= 0; --$k) {
215					for ($j = 0; $j < $nx; ++$j) {
216						$X[$k][$j] /= $this->Rdiag[$k];
217					}
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];
221						}
222					}
223				}
224				$X = new PHPExcel_Shared_JAMA_Matrix($X);
225				return ($X->getMatrix(0, $this->n-1, 0, $nx));
226			} else {
227				throw new Exception(self::MatrixRankException);
228			}
229		} else {
230			throw new Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
231		}
232	}	//	function solve()
233
234}	//	PHPExcel_Shared_JAMA_class QRDecomposition