QRDecomposition.php in Loft Data Grids 7.2

File

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

<?php

/**
 *	@package JAMA
 *
 *	For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
 *	orthogonal matrix Q and an n-by-n upper triangular matrix R so that
 *	A = Q*R.
 *
 *	The QR decompostion always exists, even if the matrix does not have
 *	full rank, so the constructor will never fail.  The primary use of the
 *	QR decomposition is in the least squares solution of nonsquare systems
 *	of simultaneous linear equations.  This will fail if isFullRank()
 *	returns false.
 *
 *	@author  Paul Meagher
 *	@license PHP v3.0
 *	@version 1.1
 */
class PHPExcel_Shared_JAMA_QRDecomposition {
  const MatrixRankException = "Can only perform operation on full-rank matrix.";

  /**
   *	Array for internal storage of decomposition.
   *	@var array
   */
  private $QR = array();

  /**
   *	Row dimension.
   *	@var integer
   */
  private $m;

  /**
   *	Column dimension.
   *	@var integer
   */
  private $n;

  /**
   *	Array for internal storage of diagonal of R.
   *	@var  array
   */
  private $Rdiag = array();

  /**
   *	QR Decomposition computed by Householder reflections.
   *
   *	@param matrix $A Rectangular matrix
   *	@return Structure to access R and the Householder vectors and compute Q.
   */
  public function __construct($A) {
    if ($A instanceof PHPExcel_Shared_JAMA_Matrix) {

      // Initialize.
      $this->QR = $A
        ->getArrayCopy();
      $this->m = $A
        ->getRowDimension();
      $this->n = $A
        ->getColumnDimension();

      // Main loop.
      for ($k = 0; $k < $this->n; ++$k) {

        // Compute 2-norm of k-th column without under/overflow.
        $nrm = 0.0;
        for ($i = $k; $i < $this->m; ++$i) {
          $nrm = hypo($nrm, $this->QR[$i][$k]);
        }
        if ($nrm != 0.0) {

          // Form k-th Householder vector.
          if ($this->QR[$k][$k] < 0) {
            $nrm = -$nrm;
          }
          for ($i = $k; $i < $this->m; ++$i) {
            $this->QR[$i][$k] /= $nrm;
          }
          $this->QR[$k][$k] += 1.0;

          // Apply transformation to remaining columns.
          for ($j = $k + 1; $j < $this->n; ++$j) {
            $s = 0.0;
            for ($i = $k; $i < $this->m; ++$i) {
              $s += $this->QR[$i][$k] * $this->QR[$i][$j];
            }
            $s = -$s / $this->QR[$k][$k];
            for ($i = $k; $i < $this->m; ++$i) {
              $this->QR[$i][$j] += $s * $this->QR[$i][$k];
            }
          }
        }
        $this->Rdiag[$k] = -$nrm;
      }
    }
    else {
      throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
    }
  }

  //	function __construct()

  /**
   *	Is the matrix full rank?
   *
   *	@return boolean true if R, and hence A, has full rank, else false.
   */
  public function isFullRank() {
    for ($j = 0; $j < $this->n; ++$j) {
      if ($this->Rdiag[$j] == 0) {
        return false;
      }
    }
    return true;
  }

  //	function isFullRank()

  /**
   *	Return the Householder vectors
   *
   *	@return Matrix Lower trapezoidal matrix whose columns define the reflections
   */
  public function getH() {
    for ($i = 0; $i < $this->m; ++$i) {
      for ($j = 0; $j < $this->n; ++$j) {
        if ($i >= $j) {
          $H[$i][$j] = $this->QR[$i][$j];
        }
        else {
          $H[$i][$j] = 0.0;
        }
      }
    }
    return new PHPExcel_Shared_JAMA_Matrix($H);
  }

  //	function getH()

  /**
   *	Return the upper triangular factor
   *
   *	@return Matrix upper triangular factor
   */
  public function getR() {
    for ($i = 0; $i < $this->n; ++$i) {
      for ($j = 0; $j < $this->n; ++$j) {
        if ($i < $j) {
          $R[$i][$j] = $this->QR[$i][$j];
        }
        elseif ($i == $j) {
          $R[$i][$j] = $this->Rdiag[$i];
        }
        else {
          $R[$i][$j] = 0.0;
        }
      }
    }
    return new PHPExcel_Shared_JAMA_Matrix($R);
  }

  //	function getR()

  /**
   *	Generate and return the (economy-sized) orthogonal factor
   *
   *	@return Matrix orthogonal factor
   */
  public function getQ() {
    for ($k = $this->n - 1; $k >= 0; --$k) {
      for ($i = 0; $i < $this->m; ++$i) {
        $Q[$i][$k] = 0.0;
      }
      $Q[$k][$k] = 1.0;
      for ($j = $k; $j < $this->n; ++$j) {
        if ($this->QR[$k][$k] != 0) {
          $s = 0.0;
          for ($i = $k; $i < $this->m; ++$i) {
            $s += $this->QR[$i][$k] * $Q[$i][$j];
          }
          $s = -$s / $this->QR[$k][$k];
          for ($i = $k; $i < $this->m; ++$i) {
            $Q[$i][$j] += $s * $this->QR[$i][$k];
          }
        }
      }
    }

    /*
    for($i = 0; $i < count($Q); ++$i) {
    	for($j = 0; $j < count($Q); ++$j) {
    		if(! isset($Q[$i][$j]) ) {
    			$Q[$i][$j] = 0;
    		}
    	}
    }
    */
    return new PHPExcel_Shared_JAMA_Matrix($Q);
  }

  //	function getQ()

  /**
   *	Least squares solution of A*X = B
   *
   *	@param Matrix $B A Matrix with as many rows as A and any number of columns.
   *	@return Matrix Matrix that minimizes the two norm of Q*R*X-B.
   */
  public function solve($B) {
    if ($B
      ->getRowDimension() == $this->m) {
      if ($this
        ->isFullRank()) {

        // Copy right hand side
        $nx = $B
          ->getColumnDimension();
        $X = $B
          ->getArrayCopy();

        // Compute Y = transpose(Q)*B
        for ($k = 0; $k < $this->n; ++$k) {
          for ($j = 0; $j < $nx; ++$j) {
            $s = 0.0;
            for ($i = $k; $i < $this->m; ++$i) {
              $s += $this->QR[$i][$k] * $X[$i][$j];
            }
            $s = -$s / $this->QR[$k][$k];
            for ($i = $k; $i < $this->m; ++$i) {
              $X[$i][$j] += $s * $this->QR[$i][$k];
            }
          }
        }

        // Solve R*X = Y;
        for ($k = $this->n - 1; $k >= 0; --$k) {
          for ($j = 0; $j < $nx; ++$j) {
            $X[$k][$j] /= $this->Rdiag[$k];
          }
          for ($i = 0; $i < $k; ++$i) {
            for ($j = 0; $j < $nx; ++$j) {
              $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
            }
          }
        }
        $X = new PHPExcel_Shared_JAMA_Matrix($X);
        return $X
          ->getMatrix(0, $this->n - 1, 0, $nx);
      }
      else {
        throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
      }
    }
    else {
      throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
    }
  }

}

//	PHPExcel_Shared_JAMA_class QRDecomposition

Classes

Name	Description
PHPExcel_Shared_JAMA_QRDecomposition	@package JAMA