[Summary view]

Differences Between: [Versions 310 and 400] [Versions 311 and 400] [Versions 39 and 400]

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