[ index.php ] [ docs.php ] [ package.php ] [ test.php ] [ example.php ] [ download.php ]

Source Listing:


Viewing: 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 BSD
* @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->$Arg->getRowDimension();
    
$this->$Arg->getColumnDimension();
    
$nu      min($this->m$this->n);
    
$e       = array();
    
$work    = array();
    
$wantu   true;
    
$wantv   true;      
    
$nct min($this->1$this->n);
    
$nrt max(0min($this->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+$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->1);
    if (
$nct $this->n)
      
$this->s[$nct] = $A[$nct][$nct];
    if (
$this->$p)
      
$this->s[$p-1] = 0.0;
    if (
$nrt $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->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 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->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->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->mmin($this->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;
  }
}
?>