Source for file Matrix.php
Documentation is available at Matrix.php
// +----------------------------------------------------------------------+
// +----------------------------------------------------------------------+
// | Copyright (c) 1997-2001 The PHP Group |
// +----------------------------------------------------------------------+
// | This source file is subject to version 2.02 of the PHP license, |
// | that is bundled with this package in the file LICENSE, and is |
// | available at through the world-wide-web at |
// | http://www.php.net/license/2_02.txt. |
// | If you did not receive a copy of the PHP license and are unable to |
// | obtain it through the world-wide-web, please send a note to |
// | license@php.net so we can mail you a copy immediately. |
// +----------------------------------------------------------------------+
// | Authors: Jesus M. Castagnetto <jmcastagnetto@php.net> |
// +----------------------------------------------------------------------+
// Matrix definition and manipulation package
// $Id: Matrix.php 304049 2010-10-05 00:30:02Z clockwerx $
require_once 'Math/Vector.php';
* Defines a matrix object.
* A matrix is implemented as an array of arrays such that:
* [0][0] [0][1] [0][2] ... [0][M]
* [1][0] [1][1] [1][2] ... [1][M]
* [N][0] [n][1] [n][2] ... [n][M]
* Originally this class was part of NumPHP (Numeric PHP package)
* @author Jesus M. Castagnetto <jmcastagnetto@php.net>
* Contains the array of arrays defining the matrix
* The number of rows in the matrix
* The number of columns in the matrix
* A flag indicating if the matrix is square
* i.e. if $this->_num_cols == $this->_num_rows
* The smallest value of all matrix cells
* The biggest value of all matrix cells
* The Euclidean norm for the matrix: sqrt(sum(e[i][j]^2))
* Cutoff error used to test for singular or ill-conditioned matrices
* Constructor for the matrix object
* @param array|Math_Matrix $data a numeric array of arrays of a Math_Matrix object
* @return object Math_Matrix
* Validates the data and initializes the internal variables (except for the determinant).
* The validation is performed by by checking that
* each row (first dimension in the array of arrays)
* contains the same number of colums (e.g. arrays of the
* @param array $data array of arrays of numbers or a valid Math_Matrix object
* @return boolean|PEAR_Error true on success, a PEAR_Error object otherwise
$this->_data = $data->getData ();
// check that we got a numeric bidimensional array
// and that all rows are of the same size
for ($i=0; $i < $nr; $i++ ) {
if (count($data[$i]) != $nc) {
return PEAR ::raiseError ('Invalid data, cannot create/modify matrix.'.
' Expecting an array of arrays or an initialized Math_Matrix object');
for ($j=0; $j < $nc; $j++ ) {
return PEAR ::raiseError ('Invalid data, cannot create/modify matrix.'.
' Expecting an array of arrays or an initialized Math_Matrix object');
$data[$i][$j] = (float) $data[$i][$j];
$eucnorm += $data[$i][$j] * $data[$i][$j];
$this->_square = ($nr == $nc);
$this->_min = !empty ($tmp)? min($tmp) : null;
$this->_max = !empty ($tmp)? max($tmp) : null;
$this->_norm = sqrt($eucnorm);
$this->_det = null; // lazy initialization ;-)
return PEAR ::raiseError ('Invalid data, cannot create/modify matrix.'.
' Expecting an array of arrays or an initialized Math_Matrix object');
* Returns the array of arrays.
* @return array|PEAR_Erroran array of array of numbers on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Sets the threshold to consider a numeric value as zero:
* if number <= epsilon then number = 0
* @param number $epsilon the upper bound value
* @return boolean|PEAR_Errortrue if successful, a PEAR_Error otherwise
return PEAR ::raisError ('Expection a number for threshold, using the old value: '. $this->_epsilon);
$this->_epsilon = $epsilon;
* Returns the value of the upper bound used to minimize round off errors
* Checks if the matrix has been initialized.
* @return boolean TRUE on success, FALSE otherwise
return ( empty ($this->_data) || is_null($this->_data) );
* Returns an array with the number of rows and columns in the matrix
* @return array|PEAR_Error an array of integers on success, a PEAR_Error object otherwise
return PEAR ::raiseError ('Matrix has not been populated');
return array ($this->_num_rows, $this->_num_cols);
* Checks if it is a square matrix (i.e. num rows == num cols)
* @return boolean TRUE on success, FALSE otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Returns the Euclidean norm of the matrix.
* Euclidean norm = sqrt( sum( e[i][j]^2 ) )
* @return float|PEAR_Errora number on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Uninitialized Math_Matrix object');
* Returns a new Math_Matrix object with the same data as the current one
* @return object Math_Matrix|PEAR_Errora Math_Matrix objects on succes, a
return PEAR ::raiseError ('Matrix has not been populated');
* Sets the value of the element at (row,col)
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
if ($row >= $this->_num_rows && $col >= $this->_num_cols) {
return PEAR ::raiseError ('Incorrect row and column values');
return PEAR ::raiseError ('Incorrect value, expecting a number');
$this->_data[$row][$col] = $value;
* Returns the value of the element at (row,col)
* @return number|PEAR_Errora number on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
if ($row >= $this->_num_rows && $col >= $this->_num_cols) {
return PEAR ::raiseError ('Incorrect row and column values');
return $this->_data[$row][$col];
* Returns the row with the given index
* This method checks that matrix has been initialized and that the
* row requested is not outside the range of rows.
* @param optional boolean $asVector whether to return a Math_Vector or a simple array. Default = false.
* @return array|Math_Vector|PEAR_Erroran array of numbers or a Math_Vector on success, a PEAR_Error otherwise
function getRow ($row, $asVector = false ) {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
if (is_integer($row) && $row >= $this->_num_rows) {
return PEAR ::raiseError ('Incorrect row value');
if (!in_array("math_vector", $classes) || !in_array("math_vectopop", $classes)) {
return PEAR ::raiseError ("Classes Math_Vector and Math_VectorOp undefined".
" add \"require_once 'Math/Vector/Vector.php'\" to your script");
return new Math_Vector ($this->_data[$row]);
return $this->_data[$row];
* Sets the row with the given index to the array
* This method checks that the row is less than the size of the matrix
* rows, and that the array size equals the number of columns in the
* @param integer $row index of the row
* @param array $arr array of numbers
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
function setRow ($row, $arr) {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
if ($row >= $this->_num_rows) {
return PEAR ::raiseError ('Row index out of bounds');
if (count($arr) != $this->_num_cols) {
return PEAR ::raiseError ('Incorrect size for matrix row: expecting '. $this->_num_cols
. ' columns, got '. count($arr). ' columns');
for ($i=0; $i < $this->_num_cols; $i++ ) {
return PEAR ::raiseError ('Incorrect values, expecting numbers');
$this->_data[$row] = $arr;
* Returns the column with the given index
* This method checks that matrix has been initialized and that the
* column requested is not outside the range of column.
* @param optional boolean $asVector whether to return a Math_Vector or a simple array. Default = false.
* @return array|Math_Vector|PEAR_Erroran array of numbers or a Math_Vector on success, a PEAR_Error otherwise
function getCol ($col, $asVector=false ) {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
if (is_integer($col) && $col >= $this->_num_cols) {
return PEAR ::raiseError ('Incorrect column value');
for ($i=0; $i < $this->_num_rows; $i++ ) {
if (!in_array("math_vector", $classes) || !in_array("math_vectopop", $classes)) {
return PEAR ::raiseError ("Classes Math_Vector and Math_VectorOp undefined".
" add \"require_once 'Math/Vector/Vector.php'\" to your script");
return new Math_Vector ($ret);
* Sets the column with the given index to the array
* This method checks that the column is less than the size of the matrix
* columns, and that the array size equals the number of rows in the
* @param integer $col index of the column
* @param array $arr array of numbers
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
function setCol ($col, $arr) {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
if ($col >= $this->_num_cols) {
return PEAR ::raiseError ('Incorrect column value');
if (count($arr) != $this->_num_cols) {
return PEAR ::raiseError ('Incorrect size for matrix column');
for ($i=0; $i < $this->_num_rows; $i++ ) {
return PEAR ::raiseError ('Incorrect values, expecting numbers');
if (PEAR ::isError ($err)) {
* Swaps the rows with the given indices
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
if (PEAR ::isError ($r1)) {
if (PEAR ::isError ($r2)) {
* Swaps the columns with the given indices
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
if (PEAR ::isError ($r1)) {
if (PEAR ::isError ($r2)) {
* Swaps a given row with a given column. Only valid for square matrices.
* @param integer $row index of row
* @param integer $col index of column
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
return PEAR ::raiseError ("Parameters must be row and a column indices");
* Returns the minimum value of the elements in the matrix
* @return number|PEAR_Errora number on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Returns the maximum value of the elements in the matrix
* @return number|PEAR_Errora number on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Gets the position of the first element with the given value
* @return array|PEAR_Erroran array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
for ($i=0; $i < $this->_num_rows; $i++ ) {
for ($j=0; $j < $this->_num_cols; $j++ ) {
if ($this->_data[$i][$j] == $val) {
* Gets the position of the element with the minimum value
* @return array|PEAR_Erroran array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Gets the position of the element with the maximum value
* @return array|PEAR_Erroran array of two numbers on success, FALSE if value is not found, and PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
* Transpose the matrix rows and columns
* @return boolean|PEAR_ErrorTRUE on success, PEAR_Error otherwise
/* John Pye noted that this operation is defined for
if (!$this->isSquare()) {
return PEAR::raiseError("Transpose is undefined for non-sqaure matrices");
for ($i=0; $i < $nc; $i++ ) {
if (PEAR ::isError ($col)) {
* Returns the trace of the matrix. Trace = sum(e[i][j]), for all i == j
* @return number|PEAR_Errora number on success, PEAR_Error otherwise
function trace() {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
return PEAR ::raiseError ('Trace undefined for non-square matrices');
for ($i=0; $i < $this->_num_rows; $i++ ) {
* Calculates the matrix determinant using Gaussian elimination with partial pivoting.
* At each step of the pivoting proccess, it checks that the normalized
* determinant calculated so far is less than 10^-18, trying to detect
* singular or ill-conditioned matrices
* @return number|PEAR_Errora number on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
return PEAR ::raiseError ('Determinant undefined for non-square matrices');
if (PEAR ::isError ($norm)) {
list ($nr, $nc) = $m->getSize ();
for ($r=0; $r< $nr; $r++ ) {
// find the maximum element in the column under the current diagonal element
$ridx = $m->_maxElementIndex ($r);
if (PEAR ::isError ($ridx)) {
$e = $m->swapRows ($r, $ridx);
$pelement = $m->getElement ($r, $r);
if (PEAR ::isError ($pelement)) {
// Is this an singular or ill-conditioned matrix?
// i.e. is the normalized determinant << 1 and -> 0?
if ((abs($det)/ $norm) < $this->_epsilon) {
return PEAR ::raiseError ('Probable singular or ill-conditioned matrix, normalized determinant = '
return PEAR ::raiseError ('Cannot continue, pivoting element is zero');
// zero all elements in column below the pivoting element
for ($i = $r + 1; $i < $nr; $i++ ) {
$factor = $m->getElement ($i, $r) / $pelement;
for ($j= $r; $j < $nc; $j++ ) {
$val = $m->getElement ($i, $j) - $factor* $m->getElement ($r, $j);
$e = $m->setElement ($i, $j, $val);
//echo $m->toString()."\n";
* Returns the normalized determinant = abs(determinant)/(euclidean norm)
* @return number|PEAR_Errora positive number on success, a PEAR_Error otherwise
if (PEAR ::isError ($det)) {
if (PEAR ::isError ($norm)) {
return PEAR ::raiseError ('Undefined normalized determinant, euclidean norm is zero');
return abs($det / $norm);
* Inverts a matrix using Gauss-Jordan elimination with partial pivoting
* @return number|PEAR_Errorthe value of the matrix determinant on success, PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
return PEAR ::raiseError ('Determinant undefined for non-square matrices');
// work on a copy to be safe
list ($nr, $nc) = $m->getSize ();
// Unit matrix to use as target
for ($i=0; $i< $nr; $i++ ) {
$e = $m->swapRows ($i, $ridx);
$e = $q->swapRows ($i, $ridx);
$pelement = $m->getElement ($i, $i);
if (PEAR ::isError ($pelement)) {
return PEAR ::raiseError ('Cannot continue inversion, pivoting element is zero');
if ((abs($det)/ $norm) < $this->_epsilon) {
return PEAR ::raiseError ('Probable singular or ill-conditioned matrix, normalized determinant = '
$e = $m->scaleRow ($i, 1/ $pelement);
$e = $q->scaleRow ($i, 1/ $pelement);
// zero all column elements execpt for the current one
$tpelement = $m->getElement ($i, $i);
for ($j=0; $j< $nr; $j++ ) {
$factor = $m->getElement ($j, $i) / $tpelement;
for ($k=0; $k< $nc; $k++ ) {
$vm = $m->getElement ($j, $k) - $factor * $m->getElement ($i, $k);
$vq = $q->getElement ($j, $k) - $factor * $q->getElement ($i, $k);
$m->setElement ($j, $k, $vm);
$q->setElement ($j, $k, $vq);
echo $m->toString()."\n";
echo $q->toString()."\n";
echo $m->toString()."\n";
echo $q->toString()."\n";
* Returns a submatrix from the position (row, col), with nrows and ncols
* @return object Math_Matrix|PEAR_ErrorMath_Matrix on success, PEAR_Error otherwise
function &getSubMatrix ($row, $col, $nrows, $ncols) {/*{{{*/
return PEAR ::raiseError ('Parameters must be a initial row and column, and number of rows and columns in submatrix');
if ($row + $nrows > $nr) {
return PEAR ::raiseError ('Rows in submatrix more than in original matrix');
if ($col + $ncols > $nc) {
return PEAR ::raiseError ('Columns in submatrix more than in original matrix');
for ($i=0; $i < $nrows; $i++ ) {
for ($j=0; $j < $ncols; $j++ ) {
$data[$i][$j] = $this->getElement($i + $row, $j + $col);
* Returns the diagonal of a square matrix as a Math_Vector
* @return object Math_Vector|PEAR_ErrorMath_Vector on success, PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
return PEAR ::raiseError ('Cannot get diagonal vector of a non-square matrix');
for ($i=0; $i< $n; $i++ ) {
return new Math_Vector ($vals);
* Returns a simple string representation of the matrix
* @param optional string $format a sprintf() format used to print the matrix elements (default = '%6.2f')
* @return string|PEAR_Errora string on success, PEAR_Error otherwise
function toString ($format= '%6.2f') {/*{{{*/
return PEAR ::raiseError ('Matrix has not been populated');
for ($i=0; $i < $this->_num_rows; $i++ ) {
for ($j=0; $j < $this->_num_cols; $j++ ) {
// remove the -0.0 output
$entry = $this->_data[$i][$j];
if (sprintf('%2.1f',$entry) == '-0.0') {
* Returns an HTML table representation of the matrix elements
* @return a string on success, PEAR_Error otherwise
return PEAR ::raiseError ('Matrix has not been populated');
$out = "<table border>\n\t<caption align=\"top\"><b>Matrix</b>";
$out .= "</caption>\n\t<tr align=\"center\">\n\t\t<th>";
$out .= $this->_num_rows. "x". $this->_num_cols. "</th>";
for ($i=0; $i < $this->_num_cols; $i++ ) {
$out .= "<th>". $i. "</th>";
for ($i=0; $i < $this->_num_rows; $i++ ) {
$out .= "\t<tr align=\"center\">\n\t\t<th>". $i. "</th>";
for ($j=0; $j < $this->_num_cols; $j++ ) {
$out .= "<td bgcolor=\"#ffffdd\">". $this->_data[$i][$j]. "</td>";
return $out. "\n</table>\n";
* Returns the index of the row with the maximum value under column of the element e[i][i]
for ($i= $r; $i< $nr; $i++ ) {
$val = abs($this->_data[$i][$r]);
* Adds a matrix to this one
* @param object Math_Matrix $m1
* @return boolean|PEAR_ErrorTRUE on success, PEAR_Error otherwise
function add ($m1) {/*{{{*/
return PEAR ::raiseError ("Parameter must be a Math_Matrix object");
if ($this->getSize() != $m1->getSize ()) {
return PEAR ::raiseError ("Matrices must have the same dimensions");
list ($nr, $nc) = $m1->getSize ();
for ($i=0; $i < $nr; $i++ ) {
for ($j=0; $j < $nc; $j++ ) {
$el1 = $m1->getElement ($i,$j);
if (PEAR ::isError ($el1)) {
if (PEAR ::isError ($el)) {
$data[$i][$j] = $el + $el1;
return PEAR ::raiseError ('Undefined error');
* Substracts a matrix from this one
* @param object Math_Matrix $m1
* @return boolean|PEAR_ErrorTRUE on success, PEAR_Error otherwise
function sub (&$m1) {/*{{{*/
return PEAR ::raiseError ("Parameter must be a Math_Matrix object");
if ($this->getSize() != $m1->getSize ()) {
return PEAR ::raiseError ("Matrices must have the same dimensions");
list ($nr, $nc) = $m1->getSize ();
for ($i=0; $i < $nr; $i++ ) {
for ($j=0; $j < $nc; $j++ ) {
$el1 = $m1->getElement ($i,$j);
if (PEAR ::isError ($el1)) {
if (PEAR ::isError ($el)) {
$data[$i][$j] = $el - $el1;
return PEAR ::raiseError ('Undefined error');
* Scales the matrix by a given factor
* @param numeric $scale the scaling factor
* @return boolean|PEAR_ErrorTRUE on success, PEAR_Error otherwise
function scale ($scale) {/*{{{*/
return PEAR ::raiseError ("Parameter must be a number");
for ($i=0; $i < $nr; $i++ ) {
for ($j=0; $j < $nc; $j++ ) {
$data[$i][$j] = $scale * $this->getElement($i,$j);
return PEAR ::raiseError ('Undefined error');
* Multiplies (scales) a row by the given factor
* @param integer $row the row index
* @param numeric $factor the scaling factor
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
function scaleRow($row, $factor) {/*{{{*/
return PEAR ::raiseError ('Uninitialized Math_Matrix object');
return PEAR ::raiseError ('Row index must be an integer, and factor a valid number');
if ($row >= $this->_num_rows) {
return PEAR ::raiseError ('Row index out of bounds');
for ($i=0; $i< $nr; $i++ ) {
return $this->setRow($row, $r);
* Multiplies this matrix (A) by another one (B), and stores
* @param object Math_Matrix $m1
* @return boolean|PEAR_ErrorTRUE on success, PEAR_Error otherwise
* @see setZeroThreshold()
return PEAR ::raiseError ('Wrong parameter, expected a Math_Matrix object');
list ($nrA, $ncA) = $this->getSize();
list ($nrB, $ncB) = $B->getSize ();
return PEAR ::raiseError ('Incompatible sizes columns in matrix must be the same as rows in parameter matrix');
for ($i=0; $i < $nrA; $i++ ) {
for ($j=0; $j < $ncB; $j++ ) {
for ($k=0; $k < $ncA; $k++ ) {
$rctot += $this->getElement($i,$k) * $B->getElement ($k, $j);
// take care of some round-off errors
if (abs($rctot) <= $this->_epsilon) {
return PEAR ::raiseError ('Undefined error');
* Multiplies a vector by this matrix
* @param object Math_Vector $v1
* @return object Math_Vector|PEAR_ErrorMath_Vector on success, PEAR_Error otherwise
* @see Math_Vector::get()
if (!Math_VectorOp ::isVector ($v1)) {
return PEAR ::raiseError ("Wrong parameter, a Math_Vector object");
return PEAR ::raiseError (" Incompatible number of columns in matrix ($nc) must ".
" be the same as the number of elements ($nv) in the vector" );
for ($i=0; $i < $nr; $i++ ) {
for ($j=0; $j < $nv; $j++ ) {
$data[$i] += $this->getElement($i,$j) * $v1->get ($j);
$obj = new Math_Vector ($data);
* Create a matrix from a file, using data stored in the given format
function &readFromFile ($filename, $format= 'serialized') {/*{{{*/
return PEAR ::raiseError ('File cannot be opened for reading');
return PEAR ::raiseError ('File is empty');
if ($format == 'serialized') {
return PEAR ::raiseError ('File did not contain a Math_Matrix object');
} else { // assume CSV data
$lines = file($filename);
foreach ($lines as $line) {
* Writes matrix object to a file using the given format
* @param object Math_Matrix $matrix the matrix object to store
* @param string $filename name of file to contain the matrix data
* @param optional string $format one of 'serialized' (default) or 'csv'
* @return boolean|PEAR_ErrorTRUE on success, a PEAR_Error otherwise
function writeToFile($matrix, $filename, $format= 'serialized') {/*{{{*/
return PEAR ::raiseError ("Parameter must be a Math_Matrix object");
if ($matrix->isEmpty ()) {
return PEAR ::raiseError ("Math_Matrix object is empty");
if ($format == 'serialized') {
list ($nr, $nc) = $matrix->getSize ();
for ($i=0; $i< $nr; $i++ ) {
$row = $matrix->getRow ($i);
if (PEAR ::isError ($row)) {
$fp = fopen($filename, "w");
return PEAR ::raiseError (" Cannot write matrix to file $filename" );
* Checks if the object is a Math_Matrix instance
* @param object Math_Matrix $matrix
* @return boolean TRUE on success, FALSE otherwise
* Returns a Math_Matrix object of size (nrows, ncols) filled with a value
* @param integer $nrows number of rows in the generated matrix
* @param integer $ncols number of columns in the generated matrix
* @param numeric $value the fill value
* @return object Math_Matrix|PEAR_ErrorMath_Matrix instance on success, PEAR_Error otherwise
function &makeMatrix ($nrows, $ncols, $value) {/*{{{*/
return PEAR ::raiseError ('Number of rows, columns, and a numeric fill value expected');
for ($i=0; $i< $nrows; $i++ ) {
* Returns the Math_Matrix object of size (nrows, ncols), filled with the value 1 (one)
* @param integer $nrows number of rows in the generated matrix
* @param integer $ncols number of columns in the generated matrix
* @return object Math_Matrix|PEAR_ErrorMath_Matrix instance on success, PEAR_Error otherwise
* @see Math_Matrix::makeMatrix()
function &makeOne ($nrows, $ncols) {/*{{{*/
* Returns the Math_Matrix object of size (nrows, ncols), filled with the value 0 (zero)
* @param integer $nrows number of rows in the generated matrix
* @param integer $ncols number of columns in the generated matrix
* @return object Math_Matrix|PEAR_ErrorMath_Matrix instance on success, PEAR_Error otherwise
* @see Math_Matrix::makeMatrix()
function &makeZero ($nrows, $ncols) {/*{{{*/
* Returns a square unit Math_Matrix object of the given size
* A unit matrix is one in which the elements follow the rules:
* Such a matrix is also called an 'identity matrix'
* @param integer $size number of rows and columns in the generated matrix
* @return object Math_Matrix|PEAR_Errora square unit Math_Matrix instance on success, PEAR_Error otherwise
* @see Math_Matrix::makeIdentity()
return PEAR ::raiseError ('An integer expected for the size of the Identity matrix');
for ($i=0; $i< $size; $i++ ) {
for ($j=0; $j< $size; $j++ ) {
$data[$i][$j] = (float) 1.0;
$data[$i][$j] = (float) 0.0;
* Returns the identity matrix of the given size. An alias of Math_Matrix::makeUnit()
* @param integer $size number of rows and columns in the generated matrix
* @return object Math_Matrix|PEAR_Errora square unit Math_Matrix instance on success, PEAR_Error otherwise
* @see Math_Matrix::makeUnit()
* Returns a Hilbert matrix of the given size: H(i,j) = 1 / (i + j - 1) where {i,j = 1..n}
* @param integer $size number of rows and columns in the Hilbert matrix
* @return object Math_Matrix|PEAR_Errora Hilber matrix on success, a PEAR_Error otherwise
return PEAR ::raiseError ('An integer expected for the size of the Hilbert matrix');
for ($i=1; $i <= $size; $i++ ) {
for ($j=1; $j <= $size; $j++ ) {
$data[$i - 1 ][$j - 1 ] = 1 / ($i + $j - 1 );
* Returns a Hankel matrix from a array of size m (C), and (optionally) of
* an array if size n (R). C will define the first column and R the last
* row. If R is not defined, C will be used. Also, if the last element of C
* is not the same to the first element of R, the last element of C is
* H(i,j) = C(i+j-1), i+j-1 <= m
* H(i,j) = R(i+j-m), otherwise
* @param array $c first column of Hankel matrix
* @param optional array $r last row of Hankel matrix
* @return object Math_Matrix|PEAR_Errora Hankel matrix on success, a PEAR_Error otherwise
return PEAR ::raiseError ('Expecting an array of values for the first column of the Hankel matrix');
return PEAR ::raiseError ('Expecting an array of values for the last row of the Hankel matrix');
// make sure that the first element of r is the same as the last element of c
for ($i=1; $i <= $nc; $i++ ) {
for ($j=1; $j <= $nr; $j++ ) {
if (($i + $j - 1 ) <= $nc) {
$val = $c[($i + $j - 1 ) - 1 ];
$val = $r[($i + $j - $nc) - 1 ];
$data[($i - 1 )][($j - 1 )] = $val;
// methods for solving linear equations
* Solves a system of linear equations: Ax = b
* a11*x1 + a12*x2 + ... + a1n*xn = b1
* a21*x1 + a22*x2 + ... + a2n*xn = b2
* ak1*x1 + ak2*x2 + ... + akn*xn = bk
* - A is matrix of coefficients (aij, i=1..k, j=1..n),
* - b a vector of values (bi, i=1..k),
* - x the vector of unkowns (xi, i=1..n)
* - Ainv is the inverse of A
* @param object Math_Matrix $a the matrix of coefficients
* @param object Math_Vector $b the vector of values
* @return object Math_Vector|PEAR_Errora Math_Vector object on succcess, PEAR_Error otherwise
function solve($a, $b) {/*{{{*/
// check that the vector classes are defined
return PEAR ::raiseError ('Incorrect parameters, expecting a Math_Matrix and a Math_Vector');
return $a->vectorMultiply ($b);
* Solves a system of linear equations: Ax = b, using an iterative error correction algorithm
* a11*x1 + a12*x2 + ... + a1n*xn = b1
* a21*x1 + a22*x2 + ... + a2n*xn = b2
* ak1*x1 + ak2*x2 + ... + akn*xn = bk
* - A is matrix of coefficients (aij, i=1..k, j=1..n),
* - b a vector of values (bi, i=1..k),
* - x the vector of unkowns (xi, i=1..n)
* - Ainv is the inverse of A
* The error correction algorithm uses the approach that if:
* - xp is the approximate solution
* - bp the values obtained from pluging xp into the original equation
* - xadj is the adjusted value (= Ainv*(b - bp))
* therefore, we calculate iteratively new values of x using the estimated
* xadj and testing to check if we have decreased the error.
* @param object Math_Matrix $a the matrix of coefficients
* @param object Math_Vector $b the vector of values
* @return object Math_Vector|PEAR_Errora Math_Vector object on succcess, PEAR_Error otherwise
* @see Math_VectorOp::add()
* @see Math_VectorOp::substract()
* @see Math_VectorOp::length()
$ainv = $a->cloneMatrix ();
$x = $ainv->vectorMultiply ($b);
$bprime = $a->vectorMultiply ($x);
if (PEAR ::isError ($bprime)) {
$err = Math_VectorOp ::substract ($b, $bprime);
$adj = $ainv->vectorMultiply ($err);
if (PEAR ::isError ($adj)) {
$adjnorm = $adj->length ();
// compute new solutions and test for accuracy
// iterate no more than 10 times
for ($i=0; $i<10; $i++ ) {
$xnew = Math_VectorOp ::add ($x, $adj);
$bprime = $a->vectorMultiply ($xnew);
$err = Math_VectorOp ::substract ($b, $bprime);
$newadj = $ainv->vectorMultiply ($err);
$newadjnorm = $newadj->length ();
// did we improve the accuracy?
if ($newadjnorm < $adjnorm) {
} else { // we did improve the accuracy, so break;
} // end of Math_Matrix class /*}}}*/
Documentation generated on Mon, 11 Mar 2019 15:39:23 -0400 by phpDocumentor 1.4.4. PEAR Logo Copyright © PHP Group 2004.
|