vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php

   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         /**
  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