You are here

class PHPExcel_Shared_JAMA_QRDecomposition in Loft Data Grids 7.2

Same name and namespace in other branches
  1. 6.2 vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php \PHPExcel_Shared_JAMA_QRDecomposition

@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

Hierarchy

Expanded class hierarchy of PHPExcel_Shared_JAMA_QRDecomposition

File

vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php, line 19

View source
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);
    }
  }

}

Members

Namesort descending Modifiers Type Description Overrides
PHPExcel_Shared_JAMA_QRDecomposition::$m private property * Row dimension. *
PHPExcel_Shared_JAMA_QRDecomposition::$n private property * Column dimension. *
PHPExcel_Shared_JAMA_QRDecomposition::$QR private property * Array for internal storage of decomposition. *
PHPExcel_Shared_JAMA_QRDecomposition::$Rdiag private property * Array for internal storage of diagonal of R. *
PHPExcel_Shared_JAMA_QRDecomposition::getH public function * Return the Householder vectors * *
PHPExcel_Shared_JAMA_QRDecomposition::getQ public function * Generate and return the (economy-sized) orthogonal factor * *
PHPExcel_Shared_JAMA_QRDecomposition::getR public function * Return the upper triangular factor * *
PHPExcel_Shared_JAMA_QRDecomposition::isFullRank public function * Is the matrix full rank? * *
PHPExcel_Shared_JAMA_QRDecomposition::MatrixRankException constant
PHPExcel_Shared_JAMA_QRDecomposition::solve public function * Least squares solution of A*X = B * *
PHPExcel_Shared_JAMA_QRDecomposition::__construct public function * QR Decomposition computed by Householder reflections. * *