Solving a linear system in PHP | Blog

I would like to show, how to solve a system of linear equations in PHP with NumPHP.

What is NumPHP?

NumPHP is my pretty small library in pre development status. I have build it after using Numpy in the university. I have missed some features like the multiplication of matrices in PHP. Sure, PHP is a language for the web and not for science, but sometimes it can be very helpful to solve problems with the calculation of matrices. Numpy can do a lot more and NumPHP hasn't the claim to be used for science, but maybe it can be helpful to solve little problems on a webpage or to evaluate some user datas. I have called the current status pre developer, cause I just wanted to check if it would be possible and I just have implemented the most important functions for this little example. NumPHP is completely tested. Enough excuses.

System of linear equations

I have a pretty small problem as example. We have matrix A

| 11.0 , 44.0 ,  1.0 |
|  0.1 ,  0.4 ,  3.0 |
|  0.0 ,  1.0 , -1.0 |

and we have vector b

| 1.0 |
| 1.0 |
| 1.0 |

We are searching for a vecor x, so that

| 11.0 , 44.0 ,  1.0 |   | x1 |   | 1.0 |
|  0.1 ,  0.4 ,  3.0 | * | x2 | = | 1.0 |
|  0.0 ,  1.0 , -1.0 |   | x3 |   | 1.0 |

Or as linear system in the school

(1) 11.0*x1 + 44.0*x2 + 1.0*x3 = 1.0
(2)  0.1*x1 +  0.4*x2 + 3.0*x3 = 1.0
(3)                x2      -x3 = 1.0

Gaussian Elimination

Maybe not the best, but a comprehensible for this problem is the Gaussian elimination. We will use it to solve your problem.

Import NumArray

use NumPHP\Core\NumArray;

Create a matrix

$matrixA = new NumArray(
    [
        [11.0, 44.0, 1.0],
        [0.1, 0.4, 3.0],
        [0.0, 1.0, -1.0],
    ]
);

Create a vector

$vectorB = new NumArray(
    [1.0, 1.0, 1.0]
);

Gaussian elimination with pivoting

/**
 * @param NumArray $matrix
 * @param NumArray $vector
 * @return array
 */
protected static function gaussianEliminationPivoting(NumArray $matrix, NumArray $vector)
{
    $shape = $matrix->getShape();

    for ($i = 0; $i < $shape[0]; $i++) {
        // find pivo element
        $max = abs($matrix->get($i, $i)->getData());
        $maxIndex = $i;
        for ($j = $i+1; $j < $shape[0]; $j++) {
            $abs = abs($matrix->get($j, $i)->getData());
            if ($abs > $max) {
                $max = $abs;
                $maxIndex = $j;
            }
        }
        // pivoting
        if ($maxIndex !== $i) {
            // change maxIndex row with i row in $matrix
            $temp = $matrix->get($i);
            $matrix->set($i, $matrix->get($maxIndex));
            $matrix->set($maxIndex, $temp);
            // change maxIndex row with i row in $vector
            $temp = $vector->get($i);
            $vector->set($i, $vector->get($maxIndex));
            $vector->set($maxIndex, $temp);
        }
        for ($j = $i+1; $j < $shape[0]; $j++) {
            $fac = -$matrix->get($j, $i)->getData()/$matrix->get($i, $i)->getData();
            $matrix->set($j, $matrix->get($j)->add($matrix->get($i)->dot($fac)));
            $vector->set($j, $vector->get($j)->add($vector->get($i)->dot($fac)));
        }
    }

    return [
        'M' => $matrix,
        'b' => $vector
    ];
}

Back substitution

/**
 * @param NumArray $matrix
 * @param NumArray $vector
 * @return NumArray
 */
protected static function backSubstitution(NumArray $matrix, NumArray $vector)
{
    $shape = $matrix->getShape();
    $xVector = \NumPHP\Core\NumPHP::zeros($shape[0]);

    for ($i = $shape[0]-1; $i >= 0; $i--) {
        $sum = 0;
        for ($j = $i+1; $j < $shape[0]; $j++) {
            $sum += $matrix->get($i, $j)->dot($xVector->get($j))->getData();
        }
        $xVector->set($i, ($vector->get($i)->getData()-$sum)/$matrix->get($i, $i)->getData());
    }

    return $xVector;
}

Solution

We only have to combine these two functions and put matrix A and vector b inside.

$gaussianResult = self::gaussianEliminationPivoting($matrixA, $vectorB);
$xVector = self::backSubstitution($gaussianResult['M'], $gaussianResult['b']);

The complete example can be found under NumPHP-examples as Symfony Command. The command has the following output

Matrix A:
NumArray([
  [11, 44, 1],
  [0.1, 0.4, 3],
  [0, 1, -1]
])
Vector b:
NumArray([1, 1, 1])
Solution of A*x=b?
Vector x:
NumArray([-5.2644376899696, 1.3313069908815, 0.33130699088146])
Time for calculation: 0.0018529891967773 sec

Checking the solution with Numpy

Just to be sure, that this is the right solution, I have created a small Numpy script with the same example.

import time
import numpy

A = numpy.array([[11., 44., 1.], [0.1, 0.4, 3.], [0., 1., -1.]])
b = numpy.array([1., 1., 1.])

time1 = time.time()
x = numpy.linalg.solve(A, b)
timeDiff = time.time()-time1

print(x)
print "%.16f" % timeDiff

And it has the following output

[-5.26443769  1.33130699  0.33130699]
0.0000879764556885

As you can see, NumPHP can be exact, but it is incredible slow. I have detected some problems and maybe we can save some milliseconds in the future.