php/Math   
Recreational Mathematics   
   home  |  library  |  contact
 Math Notes
 Math Programming [25]
 Regression [3]
 Data Mining [17]
 Notation [6]
 Linear Algebra [9]
 Stats & Prob [15]
 Math Cognition [5]
 Space & Physics [6]
 Formulas [5]
 Fun & Games [2]
 Haskell [1]
 Bayes Theory [1]
 Site News [0]
 Math Projects [5]
 Polynomials [1]
 Calculus [9]
 Number Theory [3]
 Optimization [2]
 Financial [1]

 Math Links
 PHP/ir
 Andrew Gelman
 Chance Wiki
 Daniel Lemire
 KD Knuggets
 Social Stats
 MySQL Performance
 Hunch.net
 Matthew Hurst
 JMLR
 JSS
 Hal Daume III
 Math Notes >> Linear Algebra

PHP functions to perform linear algebra decompositons [Linear Algebra
Posted on February 11, 2012 @ 09:16:36 PM by Paul Meagher

Here is a session in Octave demonstrating how to compute the QR decomposition of a matrix:

octave-3.2.4.exe:1> m = [1, 2; 3, 4]
m =

   1   2
   3   4

octave-3.2.4.exe:2> [Q R] = qr(m)
Q =

  -0.31623  -0.94868
  -0.94868   0.31623

R =

  -3.16228  -4.42719
   0.00000  -0.63246

To emulate this functional approach (versus object-oriented approach) to decomposing a matrix in PHP, we need to wrap calls to the JAMA Matrix.php object inside a qr( ) function as follows:

<?php

// Set path to JAMA Matrix object.
require_once "../Matrix.php";

// Wrap JAMA matrix functionality inside a qr function.
function qr($m) {
  
$A  = new Matrix($m);
  
$QR $A->qr();
  
$ans['R']  = $QR->getR();
  
$ans['Q']  = $QR->getQ();
  return 
$ans;
}

// Create matrix.
$m = array(array(12),array(34));

// Qet QR decompostion of a matrix.
$ans qr($m);

// Output Q and R components of answer object.
?>

<pre>
Q = 
<?php print_r($ans['Q']); ?>
</pre>

<pre>
R = 
<?php print_r($ans['R']); ?>
</pre>

The output of this script looks like this:

Q = 
Matrix Object
(
    [A] => Array
        (
            [0] => Array
                (
                    [1] => 0.94868329805051
                    [0] => -0.31622776601684
                )

            [1] => Array
                (
                    [1] => -0.31622776601684
                    [0] => -0.94868329805051
                )

        )

    [m] => 2
    [n] => 2
)

R = 
Matrix Object
(
    [A] => Array
        (
            [0] => Array
                (
                    [0] => -3.1622776601684
                    [1] => -4.4271887242357
                )

            [1] => Array
                (
                    [0] => 0
                    [1] => 0.63245553203368
                )

        )

    [m] => 2
    [n] => 2
)

Interestingly, the answers produced by the Octave and PHP scripts differ in the signs associated with a few of the answer values. Does this matter? Well, if we take the values of Q and R output by PHP, enter them into Octave, and then multiply Q and R together, we get the original matrix back. This reconstruction of the original matrix demonstrates the correctness of the decomposition and that the decomposition of a matrix into Q and R components is not unique.

octave-3.2.4.exe:1> Q = [-0.31622, 0.94868; -0.94868, -0.31622]
Q =

  -0.31622   0.94868
  -0.94868  -0.31622

octave-3.2.4.exe:2> R = [-3.1622, -4.42718; 0, 0.63245]
R =

  -3.16220  -4.42718
   0.00000   0.63245

octave-3.2.4.exe:3> Q * R
ans =

   0.99995   1.99996
   2.99992   3.99998

Note that if I didn't reduce the precision of the Q and R matrix values, the output would be [1 2; 3 4] or the exact version of the original matrix.

This proof-of-concept demonstrates how we can develop a matlab/octave-like functional API for linear algrebra decompositions in PHP by incorporating calls to the JAMA library methods inside qr(), svd(), lu(), eig(), and chol() functions. A functional API for linear algebra would make doing linear algrebra in PHP much more direct, interactive and enjoyable. Object-oriented linear algebra libraries like JAMA play a supporting role, hidden behind a set of linear algebra functions like qr(), svd(), lu(), eig() and chol() for operating on matrices.

Permalink 

A compact gaussian elimination solver for PHP [Linear Algebra
Posted on June 12, 2007 @ 03:10:35 AM by Paul Meagher

I have looked long and hard for a compact Gaussian Elimination Solver that I might port to PHP and here is my favorite implementation so far:

<?php
/**
* @class GaussianElimination.php

* Gaussian elimination solver with partial pivoting.  Originally 
* written in Java by Robert Sedgewick and Kevin Wayne:  
*
* @see http://www.cs.princeton.edu/introcs/95linear/GaussianElimination.java.html
*
* Ported to PHP by Paul Meagher

* @modified June  12/2007 reimplemented as a GaussianElimination.php class
* @modified March 26/2007 implemented as a function gauss($A, $b)
* @version 0.3
*/
class GaussianElimination {

  
/** 
  * Smallest deviation allowed in floating point comparisons. 
  */
  
const EPSILON 1e-10;

  
/**
  * Implements gaussian elimination solver with partial pivoting.
  *
  * @param double[][] $A coefficient matrix
  * @param double[]   $b output vector
  * @return double[]  $x solution vector
  */
  
public static function solve($A$b) {    

    
// number of rows
    
$N  count($b);

    
// forward elimination
    
for ($p=0$p<$N$p++) {

      
// find pivot row and swap
      
$max $p;
      for (
$i $p+1$i $N$i++)
        if (
abs($A[$i][$p]) > abs($A[$max][$p]))
          
$max $i;
      
$temp $A[$p]; $A[$p] = $A[$max]; $A[$max] = $temp;
      
$t    $b[$p]; $b[$p] = $b[$max]; $b[$max] = $t;

      
// check if matrix is singular
      
if (abs($A[$p][$p]) <= EPSILON) die("Matrix is singular or nearly singular");

      
// pivot within A and b
      
for ($i $p+1$i $N$i++) {
        
$alpha $A[$i][$p] / $A[$p][$p];
        
$b[$i] -= $alpha $b[$p];
        for (
$j $p$j $N$j++)
          
$A[$i][$j] -= $alpha $A[$p][$j];
      }
    }

    
// zero the solution vector
    
$x array_fill(0$N-10);

    
// back substitution
    
for ($i $N 1$i >= 0$i--) {
      
$sum 0.0;
      for (
$j $i 1$j $N$j++)
        
$sum += $A[$i][$j] * $x[$j];
      
$x[$i] = ($b[$i] - $sum) / $A[$i][$i];
    }

    return 
$x;

  }

}

?>

Here is a script that illustrates usage and performs a unit test of the GaussianElimination::solve method:

<?php
/**
*  Unit test for the GaussianElimination::solve
*
* @see http://www.cs.princeton.edu/introcs/95linear/GaussianElimination.java.html
*
* Ported to PHP by Paul Meagher

* @modified June 25/2007
* @modified March 25/2007
* @version 0.2
*/

// Step 1: Solve the linear equation. 

include "GaussianElimination.php";

$A = array(array(01,  1),
           array(
24, -2),
           array(
0315));
                                 
$b = array(4236);

$g = new GaussianElimination;

$x $g->solve($A$b);

// Step 2: Test whether computed and expected solution vectors are the same.

// use machine independent comparison methods
include "Utils.php";

// define expected solution vector
$expected = array(-1.02.02.0);

// tallys used for unit testing
$expected_success 0;
$expected_failure 0;

$n count($x);

for (
$i=0$i<$n$i++) {

  
// print solution vector component
  
echo $x[$i]."<br />";

  
// compare computed and expected solution vectors
  
if (Utils::eq($x[$i], $expected[$i]))
    
$expected_success++;
  else
    
$expected_failure++;  
    
}  

echo 
"Test of <code>GaussianElimination::solve</code> method ... "

// The solve method passes the test if all elements of the computed 
// solution vector match the expected solution vector

if ($expected_success == $n)
  echo 
"<font color='red'><b>passed</b></font>";
else 
  echo 
"<font color='red'><b>failed</b></font>";

?>

The output of the solve.php unit test script looks like this:

-1
2
2
Test of GaussianElimination::solve method ... passed

Adding GaussianElimination.php class to JAMA

The JAMA package does not offer a standard Gaussian Elimination solver. This is arguably because other decompositions in the package are considered superior for most of the work you might want to use a Gaussian Elimination solver on. It would nevertheless still be a good idea to have a standard Gaussian Elimination solver around for a number of reason:

  • Textbooks and classrooms in linear algebra and numerical analysis still devote a lot of their curricula to discussing Gaussian Elimination so it would be good to have a standard solver to apply to text book examples and classroom projects.
  • Numerical analysis code often uses a basic gaussian solver as an intermediate part of the calculation. Having a standard gaussian solver available makes porting such code easier.
  • A stable guassian elimination solver might be all we need to arrive at a useful and transparent solution to a math programing problem (e.g., a spline method that uses gaussian elimination to calculate curve weights).

A natural place to put a GaussianElimination.php class is among the other JAMA decomposition methods (e.g., LUDecomposition.php, QRDecomposition, CholeskyDecompostion.php, etc...). This might be a longer term goal after the present GaussianElimination.php class has undergone further testing and development if necessary. One change that would be required would be to have the class operate on a standard JAMA matrix object ($this->A) rather than the home grown version the class currently uses. The unit testing for the class might be done a bit differently as JAMA uses vector norm tests to determine if the computed and expected solution vectors are the same.

Arbitrary precision support?

I wonder if it is necessary to add arbitrary precision support to the GaussianElimination solver? One reason for doing so is that the solver might perform better on matrices that are ill-conditioned. I would also simply be interested to see what an arbitrary precision version of GaussianElimination might look like.

Gaussian Elimination for other matrix storage formats.

The Gaussian Elimination algorithm works great with standard arrays as the storage format. If the matrix is stored in a compressed or sparse format than alternative gaussian elimination algorithms might be required that would be better adapted to the particular constraints of the matrix storage format.

Complex gaussian elimination

How would gaussian elimination on a matrix of complex numbers work? Would there be any market for gaussian elimination using matrices of complex numbers?

Visualizing and auralizing the gaussian elimination process

An interesting multimedia project might involve devising new methods to visualize and auralize the gaussian elimination process for 2D, 3D, and 4D systems of equations. One interactive mapping that might be used is the mapping between an elementrary row operation and a visual or aural quality. I have done some research and development on elementary row operations (i.e., an ElementaryRowOperations.php class) and devised a web interface to apply them to a matrix as a way to reduce the matrix to a simpler form and eventually solve it. If might interesting to explore various visualizations and auralizations of what of what each elimination step is doing to the vector space of the matrix.

Permalink 

Vector norms [Linear Algebra
Posted on January 17, 2007 @ 05:33:45 AM by Paul Meagher

The notation for a vector norm consists of a vector enclosed by a pair of double vertical lines ||(a1,..., an)|| with a subscript to indicate what type of norm ||(a1,..., an)||2 you want. The three most common vector norms are the infinity-norm, one-norm, and two-norm:

  • ||(a1,..., an)|| = max(|ai|)
  • ||(a1,..., an)||1 = |a1| + |a2| + ... + |an|
  • ||(a1,..., an)||2 = sqrt( |a1|2 + |a2|2 + ... + |an|2)

The three functions below implement these three norms:

<?php
/**
* Collection of some common vector norms
*/

/**
* Infinity norm
*
* @param array $v 
* @returns maximum of absolute values
*/ 
function normInf($v) {
  
$n count($v);
  
$biggest 0.0;
  for(
$i=0$i<$n$i++)
    if (
abs($v[$i]) > $biggest)
      
$biggest abs($v[$i]);
  return 
$biggest;
}


/**
* Norm 1
*
* @param array $v 
* @returns sum of absolute values
*/ 
function norm1($v) {
  
$n count($v);
  for(
$i=0$i<$n$i++)
    
$sum += abs($v[$i]);
  return 
$sum;
}

/**
* Norm 2
*
* @param array $v 
* @returns square root of sum of squares
*/ 
function norm2($v) {
  
$n count($v);
  for(
$i=0$i<$n$i++)
    
$sum += pow($v[$i], 2);
  return 
sqrt($sum);  
}

$v = array(13, -64);

echo 
"normInf = ".normInf($v)."<br />";
echo 
"norm1 = ".norm1($v)."<br />";
echo 
"norm2 = ".norm2($v);

// Output:

// normInf = 6
// norm1 = 14
// norm2 = 7.87400787401
?>

Exercises

1. What is the difference between a vector norm and a matrix norm?

2. What are some common matrix norms?

3. Implement some common matrix norms.

Permalink 

Barycentric coordinates for collinear and non-collinear points [Linear Algebra
Posted on January 13, 2007 @ 02:14:26 PM by Paul Meagher

Take three collinear points p1, p2, p3 where p2 is a point interior to the two end points p1 and p3. Instead of expressing p2 as the x and y coordinates of the point in a Cartesian plane R2, we can instead express p2 as a ratio of the segment lengths (p1 to p2 and p2 to p3) to the overall length (p1 to p3). These ratios are called the barycentric coordinates of p2.

The first function below computes the barycentric coordinates $bc for the interior point $in of 3 collinear points. The second function takes the computed barycentric coordinates $bc, plus the original two endpoints, and returns a point $pt expressed as x and y values. The demonstration code shows that the second inverse function return the interior point supplied to the first function (i.e., $pt=$in) which verifies that everything is working as it should.

<?php
/**
* Finds barycentric coordinates for interior point $in 
* where $e1 and $e2 are the end points and all three
* points are collinear.
*
* @param array $e1 first 2D endpoint
* @param array $e2 second 2D endpoint
* @param array $in interior 2D point
* @returns barycentric coordinates $c for interior point $in
*/
function getBarycentricCoordinates($e1$e2$in) {

  
// array to hold lengths
  
$l = array(); 

  
// compute segment lengths and overall length
  
$l[0] = sqrt(pow($in[0]-$e1[0],2)+pow($in[1]-$e1[1], 2));
  
$l[1] = sqrt(pow($e2[0]-$in[0],2)+pow($e2[1]-$in[1], 2));
  
$l[2] = $l[0]+$l[1];
  
  
// array to hold barycentric coordinates
  
$bc = array(); 
  
  
// ratios of segment lengths to overall length
  
$bc[0] = $l[0]/$l[2];
  
$bc[1] = $l[1]/$l[2];
  
  return 
$bc;     
}

/**
* Gets the interior point $pt between two endpoints
* $e1 and $e2 that corresponds to the supplied barycentric
* coordinate $bc.
*
* @param array $e1 first 2D endpoint
* @param array $e2 second 2D endpoint
* @param array $bc barycentric coordinates for interior 2D point
* @returns interior point $pt corresponding to barycentric coordinates $bc
*/
function getPointUsingBarycentricCoordinates($e1$e2$bc) {
    
$pt[0] = ($e1[0]*$bc[1])+($e2[0]*$bc[0]);
    
$pt[1] = ($e1[1]*$bc[1])+($e2[1]*$bc[0]);    
    return 
$pt;
}

// specify collinear endpoints and interior point
$e1 = array(24);
$e2 = array(88);
$in = array(6.57);

// get barycentric coordinates for interior point
$bc getBarycentricCoordinates($e1$e2$in);
echo 
"<pre>";
print_r($bc);
echo 
"</pre>";

// get interior point corresponding to barycentric coordinates
$pt getPointUsingBarycentricCoordinates($e1$e2$bc);
echo 
"<pre>";
print_r($pt);
echo 
"</pre>";
?>

We can also compute the barycentric coordinates when our three points are not collinear. Three non-collinear points form a triangle and computing the barycentric coordinates in this case involves expressing any point p interior to these three non-collinear points (p1, p2, p3) as area ratios instead of length ratios. p can be expressed as the area ratios (i.e., barycentric coordinates) multiplied by each point:

p = up1 + vp2 + wp3

Where u + v + w sum to 1 (u, v, and w are ratios of exclusive and exhaustive area patches to overall areas so naturally sum to 1 just like the barycentric coordinates for an interior collinear point does).

Permalink 

Application and implementation of dot product [Linear Algebra
Posted on January 12, 2007 @ 02:41:53 AM by Paul Meagher

A simple but useful application of the dot product is computing a course average.

Let w be the vector of test weights [0.25, 0.50, 0.25].

Let g be the vector of test grades [80, 75, 90].

We compute the course average using the dot product as follows:

w . gT = course average

(0.25 * 80) + (0.50 * 75) + (0.25 * 90) = 80

It is trivial to implement the basic dot product calculation. Here is implementation that is not so trivial:

<?php
/**
* This method calculates the dot product of two vectors. It uses unrolled
* loops for increments equal to one.  It is derived from a translation from
* FORTRAN to Java of the LINPACK function DDOT.  In the LINPACK listing DDOT
* is attributed to Jack Dongarra with a date of 3/11/78.
*
* Java translation by Steve Verrill, February 25, 1997.
* @see http://www1.fpl.fs.fed.us/Blas_f77.html
*
* PHP translation by Paul Meagher, Feb 11, 2005.
*
* @param int   $n     The order of the vectors $dx[] and $dy[]
* @param float $dx[]  vector
* @param int   $incx  The subscript increment for $dx[]
* @param float $dy[]  vector
* @param int   $incy  The subscript increment for $dy[]
*/
function ddot($n$dx$incx$dy$incy) {
  
$ddot 0.0;
  if (
$n <= 0) return $ddot;
  if ((
$incx == 1) && ($incy == 1)) {
    
// both increments equal to 1
    
$m $n 5;
    for (
$i=1$i<=$m$i++)
      
$ddot += $dx[$i] * $dy[$i];
    for (
$i=$m+1$i<=$n$i+=5)
      
$ddot += $dx[$i] * $dy[$i] + $dx[$i+1] * $dy[$i+1] +
               
$dx[$i+2] * $dy[$i+2] + $dx[$i+3] *$dy[$i+3] +
               
$dx[$i+4] * $dy[$i+4];
    return 
$ddot;
  } else {
    
// at least one increment not equal to 1
    
$ix 1;
    
$iy 1;
    if (
$incx 0$ix = (-$n+1) * $incx 1;
    if (
$incy 0$iy = (-$n+1) * $incy 1;
    for (
$i=1$i<=$n$i++) {
      
$ddot += $dx[$ix] * $dy[$iy];
      
$ix += $incx;
      
$iy += $incy;
    }
    return 
$ddot;
  }
}
?>

Exercise

Create two massive vectors and compare the time it takes to calculate a dot product using a trivial dot product implementation and the ddot implementation above.

Permalink 

Discrete dynamical systems in a nutshell [Linear Algebra
Posted on December 27, 2006 @ 03:52:47 PM by Paul Meagher

A discrete dynamical system consists of a sequence of n x 1 matrices:

u0, u1, u2,...

We can model the evolution of the u vectors via the recursive formula:

uk+1 = Auk

Where A is the step matrix and uk is the step vector.

The kth step vector uk can be found via the formula:

uk = Aku0

Where u0 is the initial step vector.

Permalink 

polyfit and printpoly functions [Linear Algebra
Posted on November 15, 2006 @ 02:15:22 PM by Paul Meagher

The polyfit and printpoly functions were largely developed by Mike Bommarito and provide a nice application for the JAMA Matrix library.

<?php
require_once "../Matrix.php";
/*
* @package JAMA
* @author Michael Bommarito
* @author Paul Meagher
* @version 0.1
*
* Function to fit an order n polynomial function through
* a series of x-y data points using least squares.
*
* @param $X array x values
* @param $Y array y values
* @param $n int order of polynomial to be used for fitting
* @returns array $coeffs of polynomial coefficients
* Pre-Conditions: the system is not underdetermined: sizeof($X) > $n+1
*/
function polyfit($X$Y$n) {  
  for (
$i 0$i sizeof($X); $i++) 
      for (
$j 0$j <= $n$j++) 
        
$A[$i][$j] = pow($X[$i], $j);
  for (
$i=0$i sizeof($Y); $i++)
    
$B[$i] = array($Y[$i]);   
  
$matrixA = new Matrix($A);
  
$matrixB = new Matrix($B);
  
$C $matrixA->solve($matrixB);
  return 
$C->getMatrix(0$n01);
}

function 
printpoly$C null ) {
  for(
$i $C->1$i >= 0$i-- ) {
    
$r $C->get($i0);
    if ( 
abs($r) <= pow(10, -9) )
      
$r 0;
    if (
$i == $C->1)
      echo 
$r "x<sup>$i</sup>";
    else if (
$i $C->1)
      echo 
" + " $r "x<sup>$i</sup>";
    else if (
$i == 0)
      echo 
" + " $r;
  }
}

$X = array(0,1,2,3,4,5);
$Y = array(4,3,12,67,228579);
$points = new Matrix(array($X$Y));
$points->toHTML();
printpoly(polyfit($X$Y4));

echo 
'<hr />';

$X = array(0,1,2,3,4,5);
$Y = array(1,2,5,10,1726);
$points = new Matrix(array($X$Y));
$points->toHTML();
printpoly(polyfit($X$Y2));

echo 
'<hr />';

$X = array(0,1,2,3,4,5,6);
$Y = array(-90,-104,-178,-252,-2611604446);
$points = new Matrix(array($X$Y));
$points->toHTML();
printpoly(polyfit($X$Y5));

echo 
'<hr />';

$X = array(0,1,2,3,4);
$Y = array(mt_rand(010), mt_rand(4080), mt_rand(240400), mt_rand(18002215), mt_rand(80009000)); 
$points = new Matrix(array($X$Y));
$points->toHTML();
printpoly(polyfit($X$Y3));
?> 

Exercises

1. Call the polyfit function in a loop where the $n argument ranges from polyfit(..,..,2) to polyfit(..,..,8). You may need to implement destructors for the looping code if you have problems embedding the polyfit function in a loop.

2. Graph the original data points and polynomial function graphs for various values of polyfit(.., .., $n).

3. Develop a Polyfit class that has a constructor which sets $X and $Y arrays and makes looping though $n values easier by sharing the same $this->X and $this->Y instance variables.

Permalink 

Polar Coordinate Functions in PHP, Mathlab, and Octave [Linear Algebra
Posted on October 15, 2006 @ 01:26:14 PM by Paul Meagher

To compute CART2SPH:

rho   = ħsqrt(x2+y2+z2),
theta = arctan(y/x),
phi   = arccos(z/sqrt[x2+y2+z2])

$rho   = sqrt(($x*$x) + ($y*$y) + ($z*$z));
$phi   = acos($z/sqrt(($x*$x) + ($y*$y) + ($z*$z)));
$theta = atan($y/$x);

To compute SPH2CART

x = rho cos(theta)sin(phi),
y = rho sin(theta)sin(phi),
z = rho cos(phi)

The full MatLab API for some of this work consists of:

CART2POL Transform Cartesian to polar coordinates.
CART2SPH Transform Cartesian to spherical coordinates.
POL2CART Transform polar to Cartesian coordinates.
SPH2CART Transform spherical to Cartesian coordinates.

If you do a search on "cart2sph" it takes you to this Octave script.

Documentation for Octave (gnu version of matlab) is quite impressive.

Permalink 

Quiz: Special Matrices [Linear Algebra
Posted on June 12, 2006 @ 01:04:35 AM by Paul Meagher

Define and provide an example of each of these special matrix types:

  1. diagonal
  2. tridiagonal
  3. upper triangular
  4. upper Hessenberg
  5. symmetric
  6. nonsingular
  7. positive-definite
  8. orthogonal

These are the most popular matrix types discussed in linear algebra (Sewell, 2005).

Permalink 

 Archive 
 


php/Math Project
© 2011. All rights reserved.