[Summary view] [Print] [Text view]

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