PageRenderTime 56ms CodeModel.GetById 25ms RepoModel.GetById 0ms app.codeStats 0ms

/PHPExcel_1.7.8-with_documentation-msoffice_format/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php

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