Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

SingularValueDecomposition.php @ 1

Last change on this file since 1 was 1, checked in by dungnv, 11 years ago

File size: 12.7 KB

Rev	Line
[1]	1	<?php
	2	/**
	3	* @package JAMA
	4	*
	5	* For an m-by-n matrix A with m >= n, the singular value decomposition is
	6	* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
	7	* an n-by-n orthogonal matrix V so that A = USV'.
	8	*
	9	* The singular values, sigma[$k] = S[$k][$k], are ordered so that
	10	* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
	11	*
	12	* The singular value decompostion always exists, so the constructor will
	13	* never fail. The matrix condition number and the effective numerical
	14	* rank can be computed from this decomposition.
	15	*
	16	* @author Paul Meagher
	17	* @license PHP v3.0
	18	* @version 1.1
	19	*/
	20	class SingularValueDecomposition {
	21
	22	/**
	23	* Internal storage of U.
	24	* @var array
	25	*/
	26	private $U = array();
	27
	28	/**
	29	* Internal storage of V.
	30	* @var array
	31	*/
	32	private $V = array();
	33
	34	/**
	35	* Internal storage of singular values.
	36	* @var array
	37	*/
	38	private $s = array();
	39
	40	/**
	41	* Row dimension.
	42	* @var int
	43	*/
	44	private $m;
	45
	46	/**
	47	* Column dimension.
	48	* @var int
	49	*/
	50	private $n;
	51
	52
	53	/**
	54	* Construct the singular value decomposition
	55	*
	56	* Derived from LINPACK code.
	57	*
	58	* @param $A Rectangular matrix
	59	* @return Structure to access U, S and V.
	60	*/
	61	public function __construct($Arg) {
	62
	63	// Initialize.
	64	$A = $Arg->getArrayCopy();
	65	$this->m = $Arg->getRowDimension();
	66	$this->n = $Arg->getColumnDimension();
	67	$nu = min($this->m, $this->n);
	68	$e = array();
	69	$work = array();
	70	$wantu = true;
	71	$wantv = true;
	72	$nct = min($this->m - 1, $this->n);
	73	$nrt = max(0, min($this->n - 2, $this->m));
	74
	75	// Reduce A to bidiagonal form, storing the diagonal elements
	76	// in s and the super-diagonal elements in e.
	77	for ($k = 0; $k < max($nct,$nrt); ++$k) {
	78
	79	if ($k < $nct) {
	80	// Compute the transformation for the k-th column and
	81	// place the k-th diagonal in s[$k].
	82	// Compute 2-norm of k-th column without under/overflow.
	83	$this->s[$k] = 0;
	84	for ($i = $k; $i < $this->m; ++$i) {
	85	$this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
	86	}
	87	if ($this->s[$k] != 0.0) {
	88	if ($A[$k][$k] < 0.0) {
	89	$this->s[$k] = -$this->s[$k];
	90	}
	91	for ($i = $k; $i < $this->m; ++$i) {
	92	$A[$i][$k] /= $this->s[$k];
	93	}
	94	$A[$k][$k] += 1.0;
	95	}
	96	$this->s[$k] = -$this->s[$k];
	97	}
	98
	99	for ($j = $k + 1; $j < $this->n; ++$j) {
	100	if (($k < $nct) & ($this->s[$k] != 0.0)) {
	101	// Apply the transformation.
	102	$t = 0;
	103	for ($i = $k; $i < $this->m; ++$i) {
	104	$t += $A[$i][$k] * $A[$i][$j];
	105	}
	106	$t = -$t / $A[$k][$k];
	107	for ($i = $k; $i < $this->m; ++$i) {
	108	$A[$i][$j] += $t * $A[$i][$k];
	109	}
	110	// Place the k-th row of A into e for the
	111	// subsequent calculation of the row transformation.
	112	$e[$j] = $A[$k][$j];
	113	}
	114	}
	115
	116	if ($wantu AND ($k < $nct)) {
	117	// Place the transformation in U for subsequent back
	118	// multiplication.
	119	for ($i = $k; $i < $this->m; ++$i) {
	120	$this->U[$i][$k] = $A[$i][$k];
	121	}
	122	}
	123
	124	if ($k < $nrt) {
	125	// Compute the k-th row transformation and place the
	126	// k-th super-diagonal in e[$k].
	127	// Compute 2-norm without under/overflow.
	128	$e[$k] = 0;
	129	for ($i = $k + 1; $i < $this->n; ++$i) {
	130	$e[$k] = hypo($e[$k], $e[$i]);
	131	}
	132	if ($e[$k] != 0.0) {
	133	if ($e[$k+1] < 0.0) {
	134	$e[$k] = -$e[$k];
	135	}
	136	for ($i = $k + 1; $i < $this->n; ++$i) {
	137	$e[$i] /= $e[$k];
	138	}
	139	$e[$k+1] += 1.0;
	140	}
	141	$e[$k] = -$e[$k];
	142	if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
	143	// Apply the transformation.
	144	for ($i = $k+1; $i < $this->m; ++$i) {
	145	$work[$i] = 0.0;
	146	}
	147	for ($j = $k+1; $j < $this->n; ++$j) {
	148	for ($i = $k+1; $i < $this->m; ++$i) {
	149	$work[$i] += $e[$j] * $A[$i][$j];
	150	}
	151	}
	152	for ($j = $k + 1; $j < $this->n; ++$j) {
	153	$t = -$e[$j] / $e[$k+1];
	154	for ($i = $k + 1; $i < $this->m; ++$i) {
	155	$A[$i][$j] += $t * $work[$i];
	156	}
	157	}
	158	}
	159	if ($wantv) {
	160	// Place the transformation in V for subsequent
	161	// back multiplication.
	162	for ($i = $k + 1; $i < $this->n; ++$i) {
	163	$this->V[$i][$k] = $e[$i];
	164	}
	165	}
	166	}
	167	}
	168
	169	// Set up the final bidiagonal matrix or order p.
	170	$p = min($this->n, $this->m + 1);
	171	if ($nct < $this->n) {
	172	$this->s[$nct] = $A[$nct][$nct];
	173	}
	174	if ($this->m < $p) {
	175	$this->s[$p-1] = 0.0;
	176	}
	177	if ($nrt + 1 < $p) {
	178	$e[$nrt] = $A[$nrt][$p-1];
	179	}
	180	$e[$p-1] = 0.0;
	181	// If required, generate U.
	182	if ($wantu) {
	183	for ($j = $nct; $j < $nu; ++$j) {
	184	for ($i = 0; $i < $this->m; ++$i) {
	185	$this->U[$i][$j] = 0.0;
	186	}
	187	$this->U[$j][$j] = 1.0;
	188	}
	189	for ($k = $nct - 1; $k >= 0; --$k) {
	190	if ($this->s[$k] != 0.0) {
	191	for ($j = $k + 1; $j < $nu; ++$j) {
	192	$t = 0;
	193	for ($i = $k; $i < $this->m; ++$i) {
	194	$t += $this->U[$i][$k] * $this->U[$i][$j];
	195	}
	196	$t = -$t / $this->U[$k][$k];
	197	for ($i = $k; $i < $this->m; ++$i) {
	198	$this->U[$i][$j] += $t * $this->U[$i][$k];
	199	}
	200	}
	201	for ($i = $k; $i < $this->m; ++$i ) {
	202	$this->U[$i][$k] = -$this->U[$i][$k];
	203	}
	204	$this->U[$k][$k] = 1.0 + $this->U[$k][$k];
	205	for ($i = 0; $i < $k - 1; ++$i) {
	206	$this->U[$i][$k] = 0.0;
	207	}
	208	} else {
	209	for ($i = 0; $i < $this->m; ++$i) {
	210	$this->U[$i][$k] = 0.0;
	211	}
	212	$this->U[$k][$k] = 1.0;
	213	}
	214	}
	215	}
	216
	217	// If required, generate V.
	218	if ($wantv) {
	219	for ($k = $this->n - 1; $k >= 0; --$k) {
	220	if (($k < $nrt) AND ($e[$k] != 0.0)) {
	221	for ($j = $k + 1; $j < $nu; ++$j) {
	222	$t = 0;
	223	for ($i = $k + 1; $i < $this->n; ++$i) {
	224	$t += $this->V[$i][$k]* $this->V[$i][$j];
	225	}
	226	$t = -$t / $this->V[$k+1][$k];
	227	for ($i = $k + 1; $i < $this->n; ++$i) {
	228	$this->V[$i][$j] += $t * $this->V[$i][$k];
	229	}
	230	}
	231	}
	232	for ($i = 0; $i < $this->n; ++$i) {
	233	$this->V[$i][$k] = 0.0;
	234	}
	235	$this->V[$k][$k] = 1.0;
	236	}
	237	}
	238
	239	// Main iteration loop for the singular values.
	240	$pp = $p - 1;
	241	$iter = 0;
	242	$eps = pow(2.0, -52.0);
	243
	244	while ($p > 0) {
	245	// Here is where a test for too many iterations would go.
	246	// This section of the program inspects for negligible
	247	// elements in the s and e arrays. On completion the
	248	// variables kase and k are set as follows:
	249	// kase = 1 if s(p) and e[k-1] are negligible and k<p
	250	// kase = 2 if s(k) is negligible and k<p
	251	// kase = 3 if e[k-1] is negligible, k<p, and
	252	// s(k), ..., s(p) are not negligible (qr step).
	253	// kase = 4 if e(p-1) is negligible (convergence).
	254	for ($k = $p - 2; $k >= -1; --$k) {
	255	if ($k == -1) {
	256	break;
	257	}
	258	if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
	259	$e[$k] = 0.0;
	260	break;
	261	}
	262	}
	263	if ($k == $p - 2) {
	264	$kase = 4;
	265	} else {
	266	for ($ks = $p - 1; $ks >= $k; --$ks) {
	267	if ($ks == $k) {
	268	break;
	269	}
	270	$t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
	271	if (abs($this->s[$ks]) <= $eps * $t) {
	272	$this->s[$ks] = 0.0;
	273	break;
	274	}
	275	}
	276	if ($ks == $k) {
	277	$kase = 3;
	278	} else if ($ks == $p-1) {
	279	$kase = 1;
	280	} else {
	281	$kase = 2;
	282	$k = $ks;
	283	}
	284	}
	285	++$k;
	286
	287	// Perform the task indicated by kase.
	288	switch ($kase) {
	289	// Deflate negligible s(p).
	290	case 1:
	291	$f = $e[$p-2];
	292	$e[$p-2] = 0.0;
	293	for ($j = $p - 2; $j >= $k; --$j) {
	294	$t = hypo($this->s[$j],$f);
	295	$cs = $this->s[$j] / $t;
	296	$sn = $f / $t;
	297	$this->s[$j] = $t;
	298	if ($j != $k) {
	299	$f = -$sn * $e[$j-1];
	300	$e[$j-1] = $cs * $e[$j-1];
	301	}
	302	if ($wantv) {
	303	for ($i = 0; $i < $this->n; ++$i) {
	304	$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
	305	$this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
	306	$this->V[$i][$j] = $t;
	307	}
	308	}
	309	}
	310	break;
	311	// Split at negligible s(k).
	312	case 2:
	313	$f = $e[$k-1];
	314	$e[$k-1] = 0.0;
	315	for ($j = $k; $j < $p; ++$j) {
	316	$t = hypo($this->s[$j], $f);
	317	$cs = $this->s[$j] / $t;
	318	$sn = $f / $t;
	319	$this->s[$j] = $t;
	320	$f = -$sn * $e[$j];
	321	$e[$j] = $cs * $e[$j];
	322	if ($wantu) {
	323	for ($i = 0; $i < $this->m; ++$i) {
	324	$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
	325	$this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
	326	$this->U[$i][$j] = $t;
	327	}
	328	}
	329	}
	330	break;
	331	// Perform one qr step.
	332	case 3:
	333	// Calculate the shift.
	334	$scale = max(max(max(max(
	335	abs($this->s[$p-1]),abs($this->s[$p-2])),abs($e[$p-2])),
	336	abs($this->s[$k])), abs($e[$k]));
	337	$sp = $this->s[$p-1] / $scale;
	338	$spm1 = $this->s[$p-2] / $scale;
	339	$epm1 = $e[$p-2] / $scale;
	340	$sk = $this->s[$k] / $scale;
	341	$ek = $e[$k] / $scale;
	342	$b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
	343	$c = ($sp * $epm1) * ($sp * $epm1);
	344	$shift = 0.0;
	345	if (($b != 0.0) \|\| ($c != 0.0)) {
	346	$shift = sqrt($b * $b + $c);
	347	if ($b < 0.0) {
	348	$shift = -$shift;
	349	}
	350	$shift = $c / ($b + $shift);
	351	}
	352	$f = ($sk + $sp) * ($sk - $sp) + $shift;
	353	$g = $sk * $ek;
	354	// Chase zeros.
	355	for ($j = $k; $j < $p-1; ++$j) {
	356	$t = hypo($f,$g);
	357	$cs = $f/$t;
	358	$sn = $g/$t;
	359	if ($j != $k) {
	360	$e[$j-1] = $t;
	361	}
	362	$f = $cs * $this->s[$j] + $sn * $e[$j];
	363	$e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
	364	$g = $sn * $this->s[$j+1];
	365	$this->s[$j+1] = $cs * $this->s[$j+1];
	366	if ($wantv) {
	367	for ($i = 0; $i < $this->n; ++$i) {
	368	$t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
	369	$this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
	370	$this->V[$i][$j] = $t;
	371	}
	372	}
	373	$t = hypo($f,$g);
	374	$cs = $f/$t;
	375	$sn = $g/$t;
	376	$this->s[$j] = $t;
	377	$f = $cs * $e[$j] + $sn * $this->s[$j+1];
	378	$this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
	379	$g = $sn * $e[$j+1];
	380	$e[$j+1] = $cs * $e[$j+1];
	381	if ($wantu && ($j < $this->m - 1)) {
	382	for ($i = 0; $i < $this->m; ++$i) {
	383	$t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
	384	$this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
	385	$this->U[$i][$j] = $t;
	386	}
	387	}
	388	}
	389	$e[$p-2] = $f;
	390	$iter = $iter + 1;
	391	break;
	392	// Convergence.
	393	case 4:
	394	// Make the singular values positive.
	395	if ($this->s[$k] <= 0.0) {
	396	$this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
	397	if ($wantv) {
	398	for ($i = 0; $i <= $pp; ++$i) {
	399	$this->V[$i][$k] = -$this->V[$i][$k];
	400	}
	401	}
	402	}
	403	// Order the singular values.
	404	while ($k < $pp) {
	405	if ($this->s[$k] >= $this->s[$k+1]) {
	406	break;
	407	}
	408	$t = $this->s[$k];
	409	$this->s[$k] = $this->s[$k+1];
	410	$this->s[$k+1] = $t;
	411	if ($wantv AND ($k < $this->n - 1)) {
	412	for ($i = 0; $i < $this->n; ++$i) {
	413	$t = $this->V[$i][$k+1];
	414	$this->V[$i][$k+1] = $this->V[$i][$k];
	415	$this->V[$i][$k] = $t;
	416	}
	417	}
	418	if ($wantu AND ($k < $this->m-1)) {
	419	for ($i = 0; $i < $this->m; ++$i) {
	420	$t = $this->U[$i][$k+1];
	421	$this->U[$i][$k+1] = $this->U[$i][$k];
	422	$this->U[$i][$k] = $t;
	423	}
	424	}
	425	++$k;
	426	}
	427	$iter = 0;
	428	--$p;
	429	break;
	430	} // end switch
	431	} // end while
	432
	433	} // end constructor
	434
	435
	436	/**
	437	* Return the left singular vectors
	438	*
	439	* @access public
	440	* @return U
	441	*/
	442	public function getU() {
	443	return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
	444	}
	445
	446
	447	/**
	448	* Return the right singular vectors
	449	*
	450	* @access public
	451	* @return V
	452	*/
	453	public function getV() {
	454	return new Matrix($this->V);
	455	}
	456
	457
	458	/**
	459	* Return the one-dimensional array of singular values
	460	*
	461	* @access public
	462	* @return diagonal of S.
	463	*/
	464	public function getSingularValues() {
	465	return $this->s;
	466	}
	467
	468
	469	/**
	470	* Return the diagonal matrix of singular values
	471	*
	472	* @access public
	473	* @return S
	474	*/
	475	public function getS() {
	476	for ($i = 0; $i < $this->n; ++$i) {
	477	for ($j = 0; $j < $this->n; ++$j) {
	478	$S[$i][$j] = 0.0;
	479	}
	480	$S[$i][$i] = $this->s[$i];
	481	}
	482	return new Matrix($S);
	483	}
	484
	485
	486	/**
	487	* Two norm
	488	*
	489	* @access public
	490	* @return max(S)
	491	*/
	492	public function norm2() {
	493	return $this->s[0];
	494	}
	495
	496
	497	/**
	498	* Two norm condition number
	499	*
	500	* @access public
	501	* @return max(S)/min(S)
	502	*/
	503	public function cond() {
	504	return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
	505	}
	506
	507
	508	/**
	509	* Effective numerical matrix rank
	510	*
	511	* @access public
	512	* @return Number of nonnegligible singular values.
	513	*/
	514	public function rank() {
	515	$eps = pow(2.0, -52.0);
	516	$tol = max($this->m, $this->n) * $this->s[0] * $eps;
	517	$r = 0;
	518	for ($i = 0; $i < count($this->s); ++$i) {
	519	if ($this->s[$i] > $tol) {
	520	++$r;
	521	}
	522	}
	523	return $r;
	524	}
	525
	526	} // class SingularValueDecomposition

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: sourcecode/application/libraries/PHPExcel/Shared/JAMA/SingularValueDecomposition.php @ 1

Download in other formats: