/branches/v1.6.5/Classes/PHPExcel/Shared/JAMA/SingularValueDecomposition.php

<?php
/**
* @package JAMA
*
* For an m-by-n matrix A with m >= n, the singular value decomposition is
* an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
* an n-by-n orthogonal matrix V so that A = U*S*V'.
*
* The singular values, sigma[$k] = S[$k][$k], are ordered so that
* sigma[0] >= sigma[1] >= ... >= sigma[n-1].
*
* The singular value decompostion always exists, so the constructor will
* never fail.  The matrix condition number and the effective numerical
* rank can be computed from this decomposition.
*
* @author  Paul Meagher
* @license PHP v3.0
* @version 1.1
*/
class SingularValueDecomposition  {

  /**
  * Internal storage of U.
  * @var array
  */
  var $U = array();

  /**
  * Internal storage of V.
  * @var array
  */
  var $V = array();

  /**
  * Internal storage of singular values.
  * @var array
  */
  var $s = array();

  /**
  * Row dimension.
  * @var int
  */
  var $m;

  /**
  * Column dimension.
  * @var int
  */
  var $n;

  /**
  * Construct the singular value decomposition
  *
  * Derived from LINPACK code.
  *
  * @param $A Rectangular matrix
  * @return Structure to access U, S and V.
  */
  function SingularValueDecomposition ($Arg) {

    // Initialize.

    $A = $Arg->getArrayCopy();
    $this->m = $Arg->getRowDimension();
    $this->n = $Arg->getColumnDimension();
    $nu      = min($this->m, $this->n);
    $e       = array();
    $work    = array();
    $wantu   = true;
    $wantv   = true;      
    $nct = min($this->m - 1, $this->n);
    $nrt = max(0, min($this->n - 2, $this->m));   

    // Reduce A to bidiagonal form, storing the diagonal elements
    // in s and the super-diagonal elements in e.

    for ($k = 0; $k < max($nct,$nrt); $k++) {

      if ($k < $nct) {
        // Compute the transformation for the k-th column and
        // place the k-th diagonal in s[$k].
        // Compute 2-norm of k-th column without under/overflow.
        $this->s[$k] = 0;
        for ($i = $k; $i < $this->m; $i++)
          $this->s[$k] = hypo($this->s[$k], $A[$i][$k]);
        if ($this->s[$k] != 0.0) {
          if ($A[$k][$k] < 0.0)
            $this->s[$k] = -$this->s[$k];
          for ($i = $k; $i < $this->m; $i++)
            $A[$i][$k] /= $this->s[$k];
          $A[$k][$k] += 1.0;
        }
        $this->s[$k] = -$this->s[$k];
      }           

      for ($j = $k + 1; $j < $this->n; $j++) {
        if (($k < $nct) & ($this->s[$k] != 0.0)) {
          // Apply the transformation.
          $t = 0;
          for ($i = $k; $i < $this->m; $i++)
            $t += $A[$i][$k] * $A[$i][$j];
          $t = -$t / $A[$k][$k];
          for ($i = $k; $i < $this->m; $i++)
            $A[$i][$j] += $t * $A[$i][$k];
          // Place the k-th row of A into e for the
          // subsequent calculation of the row transformation.
          $e[$j] = $A[$k][$j];
        }
      }

      if ($wantu AND ($k < $nct)) {
        // Place the transformation in U for subsequent back
        // multiplication.
        for ($i = $k; $i < $this->m; $i++)
          $this->U[$i][$k] = $A[$i][$k];
      }

      if ($k < $nrt) {
        // Compute the k-th row transformation and place the
        // k-th super-diagonal in e[$k].
        // Compute 2-norm without under/overflow.
        $e[$k] = 0;
        for ($i = $k + 1; $i < $this->n; $i++)
          $e[$k] = hypo($e[$k], $e[$i]);
        if ($e[$k] != 0.0) {
          if ($e[$k+1] < 0.0)
            $e[$k] = -$e[$k];
          for ($i = $k + 1; $i < $this->n; $i++)
            $e[$i] /= $e[$k];
          $e[$k+1] += 1.0;
        }
        $e[$k] = -$e[$k];
        if (($k+1 < $this->m) AND ($e[$k] != 0.0)) {
          // Apply the transformation.
          for ($i = $k+1; $i < $this->m; $i++)
            $work[$i] = 0.0;
          for ($j = $k+1; $j < $this->n; $j++)
            for ($i = $k+1; $i < $this->m; $i++)
              $work[$i] += $e[$j] * $A[$i][$j];
          for ($j = $k + 1; $j < $this->n; $j++) {
            $t = -$e[$j] / $e[$k+1];
            for ($i = $k + 1; $i < $this->m; $i++)
              $A[$i][$j] += $t * $work[$i];
          }
        }
        if ($wantv) {
          // Place the transformation in V for subsequent
          // back multiplication.
          for ($i = $k + 1; $i < $this->n; $i++)
            $this->V[$i][$k] = $e[$i];
        }        
      }
    } 
       
    // Set up the final bidiagonal matrix or order p.
    $p = min($this->n, $this->m + 1);
    if ($nct < $this->n)
      $this->s[$nct] = $A[$nct][$nct];
    if ($this->m < $p)
      $this->s[$p-1] = 0.0;
    if ($nrt + 1 < $p)
      $e[$nrt] = $A[$nrt][$p-1];
    $e[$p-1] = 0.0;
    // If required, generate U.
    if ($wantu) {
      for ($j = $nct; $j < $nu; $j++) {
        for ($i = 0; $i < $this->m; $i++)
          $this->U[$i][$j] = 0.0;
        $this->U[$j][$j] = 1.0;
      }
      for ($k = $nct - 1; $k >= 0; $k--) {
        if ($this->s[$k] != 0.0) {
          for ($j = $k + 1; $j < $nu; $j++) {
            $t = 0;
            for ($i = $k; $i < $this->m; $i++)
              $t += $this->U[$i][$k] * $this->U[$i][$j];
            $t = -$t / $this->U[$k][$k];
            for ($i = $k; $i < $this->m; $i++)
              $this->U[$i][$j] += $t * $this->U[$i][$k];
          }
          for ($i = $k; $i < $this->m; $i++ )
            $this->U[$i][$k] = -$this->U[$i][$k];
          $this->U[$k][$k] = 1.0 + $this->U[$k][$k];
          for ($i = 0; $i < $k - 1; $i++)
             $this->U[$i][$k] = 0.0;
         } else {
           for ($i = 0; $i < $this->m; $i++)
             $this->U[$i][$k] = 0.0;
           $this->U[$k][$k] = 1.0;
         }
      }
    }

    // If required, generate V.
    if ($wantv) {
      for ($k = $this->n - 1; $k >= 0; $k--) {
        if (($k < $nrt) AND ($e[$k] != 0.0)) {
          for ($j = $k + 1; $j < $nu; $j++) {
            $t = 0;
            for ($i = $k + 1; $i < $this->n; $i++)
              $t += $this->V[$i][$k]* $this->V[$i][$j];
            $t = -$t / $this->V[$k+1][$k];
            for ($i = $k + 1; $i < $this->n; $i++)
              $this->V[$i][$j] += $t * $this->V[$i][$k];
          }
        }
        for ($i = 0; $i < $this->n; $i++)
           $this->V[$i][$k] = 0.0;
        $this->V[$k][$k] = 1.0;
      }
    }
        
    // Main iteration loop for the singular values.
    $pp   = $p - 1;
    $iter = 0;
    $eps  = pow(2.0, -52.0);
    while ($p > 0) {      
      
      // Here is where a test for too many iterations would go.
      // This section of the program inspects for negligible
      // elements in the s and e arrays.  On completion the
      // variables kase and k are set as follows:
      // kase = 1  if s(p) and e[k-1] are negligible and k<p
      // kase = 2  if s(k) is negligible and k<p
      // kase = 3  if e[k-1] is negligible, k<p, and
      //           s(k), ..., s(p) are not negligible (qr step).
      // kase = 4  if e(p-1) is negligible (convergence).

      for ($k = $p - 2; $k >= -1; $k--) {
        if ($k == -1)
          break;
        if (abs($e[$k]) <= $eps * (abs($this->s[$k]) + abs($this->s[$k+1]))) {
          $e[$k] = 0.0;
          break;
        }
      }
      if ($k == $p - 2)
        $kase = 4;
      else {
        for ($ks = $p - 1; $ks >= $k; $ks--) {
          if ($ks == $k)
            break;
          $t = ($ks != $p ? abs($e[$ks]) : 0.) + ($ks != $k + 1 ? abs($e[$ks-1]) : 0.);
          if (abs($this->s[$ks]) <= $eps * $t)  {
            $this->s[$ks] = 0.0;
            break;
          }
        }
        if ($ks == $k)
           $kase = 3;
        else if ($ks == $p-1)
           $kase = 1;
        else {
           $kase = 2;
           $k = $ks;
        }
      }
      $k++;

      // Perform the task indicated by kase.
      switch ($kase) {
         // Deflate negligible s(p).
         case 1:
           $f = $e[$p-2];
           $e[$p-2] = 0.0;
           for ($j = $p - 2; $j >= $k; $j--) {
             $t  = hypo($this->s[$j],$f);
             $cs = $this->s[$j] / $t;
             $sn = $f / $t;
             $this->s[$j] = $t;
             if ($j != $k) {
               $f = -$sn * $e[$j-1];
               $e[$j-1] = $cs * $e[$j-1];
             }
             if ($wantv) {
               for ($i = 0; $i < $this->n; $i++) {
                 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$p-1];
                 $this->V[$i][$p-1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$p-1];
                 $this->V[$i][$j] = $t;
               }
             }
           }
           break;
         // Split at negligible s(k).
         case 2:
           $f = $e[$k-1];
           $e[$k-1] = 0.0;
           for ($j = $k; $j < $p; $j++) {
             $t = hypo($this->s[$j], $f);
             $cs = $this->s[$j] / $t;
             $sn = $f / $t;
             $this->s[$j] = $t;
             $f = -$sn * $e[$j];
             $e[$j] = $cs * $e[$j];
             if ($wantu) {
               for ($i = 0; $i < $this->m; $i++) {
                 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$k-1];
                 $this->U[$i][$k-1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$k-1];
                 $this->U[$i][$j] = $t;
               }
             }
           }
           break;           
         // Perform one qr step.
         case 3:                  
           // Calculate the shift.                          
           $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]));
           $sp   = $this->s[$p-1] / $scale;
           $spm1 = $this->s[$p-2] / $scale;
           $epm1 = $e[$p-2] / $scale;
           $sk   = $this->s[$k] / $scale;
           $ek   = $e[$k] / $scale;
           $b    = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
           $c    = ($sp * $epm1) * ($sp * $epm1);
           $shift = 0.0;
           if (($b != 0.0) || ($c != 0.0)) {
             $shift = sqrt($b * $b + $c);
             if ($b < 0.0)
               $shift = -$shift;
             $shift = $c / ($b + $shift);
           }
           $f = ($sk + $sp) * ($sk - $sp) + $shift;
           $g = $sk * $ek;  
           // Chase zeros.  
           for ($j = $k; $j < $p-1; $j++) {
             $t  = hypo($f,$g);
             $cs = $f/$t;
             $sn = $g/$t;
             if ($j != $k)
               $e[$j-1] = $t;
             $f = $cs * $this->s[$j] + $sn * $e[$j];
             $e[$j] = $cs * $e[$j] - $sn * $this->s[$j];
             $g = $sn * $this->s[$j+1];
             $this->s[$j+1] = $cs * $this->s[$j+1];
             if ($wantv) {
               for ($i = 0; $i < $this->n; $i++) {
                 $t = $cs * $this->V[$i][$j] + $sn * $this->V[$i][$j+1];
                 $this->V[$i][$j+1] = -$sn * $this->V[$i][$j] + $cs * $this->V[$i][$j+1];
                 $this->V[$i][$j] = $t;
               }
             }                                                                      
             $t = hypo($f,$g);
             $cs = $f/$t;
             $sn = $g/$t;
             $this->s[$j] = $t;
             $f = $cs * $e[$j] + $sn * $this->s[$j+1];
             $this->s[$j+1] = -$sn * $e[$j] + $cs * $this->s[$j+1];
             $g = $sn * $e[$j+1];
             $e[$j+1] = $cs * $e[$j+1];
             if ($wantu && ($j < $this->m - 1)) {
               for ($i = 0; $i < $this->m; $i++) {
                 $t = $cs * $this->U[$i][$j] + $sn * $this->U[$i][$j+1];
                 $this->U[$i][$j+1] = -$sn * $this->U[$i][$j] + $cs * $this->U[$i][$j+1];
                 $this->U[$i][$j] = $t;
               }
             }                
           }
           $e[$p-2] = $f;
           $iter = $iter + 1;
           break;
         // Convergence.
         case 4:
           // Make the singular values positive.
           if ($this->s[$k] <= 0.0) {
             $this->s[$k] = ($this->s[$k] < 0.0 ? -$this->s[$k] : 0.0);
             if ($wantv) {
               for ($i = 0; $i <= $pp; $i++)
                 $this->V[$i][$k] = -$this->V[$i][$k];
             }
           }
           // Order the singular values.  
           while ($k < $pp) {
            if ($this->s[$k] >= $this->s[$k+1])
              break;
            $t = $this->s[$k];
            $this->s[$k] = $this->s[$k+1];
            $this->s[$k+1] = $t;
            if ($wantv AND ($k < $this->n - 1)) {
              for ($i = 0; $i < $this->n; $i++) {
                $t = $this->V[$i][$k+1];
                $this->V[$i][$k+1] = $this->V[$i][$k];
                $this->V[$i][$k] = $t;
              }
            }
            if ($wantu AND ($k < $this->m-1)) {
              for ($i = 0; $i < $this->m; $i++) {
                $t = $this->U[$i][$k+1];
                $this->U[$i][$k+1] = $this->U[$i][$k];
                $this->U[$i][$k] = $t;
              }
            }
            $k++;
           }
           $iter = 0;
           $p--;
           break;                                   
      } // end switch      
    } // end while

    /*
    echo "<p>Output A</p>";   
    $A = new Matrix($A);        
    $A->toHTML();
    
    echo "<p>Matrix U</p>";    
    echo "<pre>";
    print_r($this->U);        
    echo "</pre>";
    
    echo "<p>Matrix V</p>";    
    echo "<pre>";
    print_r($this->V);        
    echo "</pre>";    
    
    echo "<p>Vector S</p>";    
    echo "<pre>";
    print_r($this->s);        
    echo "</pre>";    
    exit;
    */
    
  } // end constructor
  
  /**
  * Return the left singular vectors
  * @access public
  * @return U
  */
  function getU() {
    return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
  }

  /**
  * Return the right singular vectors
  * @access public   
  * @return V
  */
  function getV() {
    return new Matrix($this->V);
  }

  /**
  * Return the one-dimensional array of singular values
  * @access public   
  * @return diagonal of S.
  */
  function getSingularValues() {
    return $this->s;
  }

  /**
  * Return the diagonal matrix of singular values
  * @access public   
  * @return S
  */
  function getS() {
    for ($i = 0; $i < $this->n; $i++) {
      for ($j = 0; $j < $this->n; $j++)
        $S[$i][$j] = 0.0;
      $S[$i][$i] = $this->s[$i];
    }
    return new Matrix($S);
  }

  /**
  * Two norm
  * @access public   
  * @return max(S)
  */
  function norm2() {
    return $this->s[0];
  }

  /**
  * Two norm condition number
  * @access public   
  * @return max(S)/min(S)
  */
  function cond() {
    return $this->s[0] / $this->s[min($this->m, $this->n) - 1];
  }

  /**
  * Effective numerical matrix rank
  * @access public   
  * @return Number of nonnegligible singular values.
  */
  function rank() {
    $eps = pow(2.0, -52.0);
    $tol = max($this->m, $this->n) * $this->s[0] * $eps;
    $r = 0;
    for ($i = 0; $i < count($this->s); $i++) {
      if ($this->s[$i] > $tol)
        $r++;
    }
    return $r;
  }
}
?>