diff --git a/README.md b/README.md index 66a1de4..3d59942 100644 --- a/README.md +++ b/README.md @@ -17,34 +17,44 @@ Matrix Transform This library currently provides the following operations: - - addition - - direct sum - - subtraction - - multiplication - - division (using [A].[B]-1) - - division by - - division into - -together with functions for - - - adjoint - - antidiagonal - - cofactors - - determinant - - diagonal - - identity - - inverse - - minors - - trace - - transpose +- addition +- direct sum +- subtraction +- multiplication +- division (using [A].[B]-1) + - division by + - division into +together with functions for + +- adjoint +- antidiagonal +- cofactors +- determinant +- diagonal +- identity +- inverse +- minors +- trace +- transpose + +and classes for + +- Decomposition + - LU Decomposition with partial row pivoting, + + such that [P].[A] = [L].[U] and [A] = [P]|.[L].[U] + - QR Decomposition + + such that [A] = [Q].[R] ## TO DO - - power() - - EigenValues - - EigenVectors - - Decomposition +- power() function +- EigenValues +- EigenVectors +- Decomposition + - Cholesky Decomposition --- @@ -52,7 +62,7 @@ together with functions for To create a new Matrix object, provide an array as the constructor argument -``` +```php $grid = [ [16, 3, 2, 13], [ 5, 10, 11, 8], @@ -63,12 +73,12 @@ $grid = [ $matrix = new Matrix\Matrix($grid); ``` The `Builder` class provides helper methods for creating specific matrices, specifically an identity matrix of a specified size; or a matrix of a specified dimensions, with every cell containing a set value. -``` -$matrix = new Matrix\Builder::createFilledMatrix(1, 5, 3); +```php +$matrix = Matrix\Builder::createFilledMatrix(1, 5, 3); ``` Will create a matrix of 5 rows and 3 columns, filled with a `1` in every cell; while -``` -$matrix = new Matrix\Builder::createIdentityMatrix(3); +```php +$matrix = Matrix\Builder::createIdentityMatrix(3); ``` will create a 3x3 identity matrix. @@ -79,34 +89,34 @@ Matrix objects are immutable: whenever you call a method or pass a grid to a fun To perform mathematical operations with Matrices, you can call the appropriate method against a matrix value, passing other values as arguments -``` -$matrix1 = new Matrix([ +```php +$matrix1 = new Matrix\Matrix([ [2, 7, 6], [9, 5, 1], [4, 3, 8], ]); -$matrix2 = new Matrix([ +$matrix2 = new Matrix\Matrix([ [1, 2, 3], [4, 5, 6], [7, 8, 9], ]); -echo $matrix1->multiply($matrix2); +var_dump($matrix1->multiply($matrix2)->toArray()); ``` or pass all values to the appropriate function -``` -$matrix1 = new Matrix([ +```php +$matrix1 = new Matrix\Matrix([ [2, 7, 6], [9, 5, 1], [4, 3, 8], ]); -$matrix2 = new Matrix([ +$matrix2 = new Matrix\Matrix([ [1, 2, 3], [4, 5, 6], [7, 8, 9], ]); -echo Matrix\multiply($matrix1, $matrix2); +var_dump(Matrix\multiply($matrix1, $matrix2)->toArray()); ``` You can pass in the arguments as Matrix objects, or as arrays. @@ -115,7 +125,7 @@ If you want to perform the same operation against multiple values (e.g. to add t ## Using functions When calling any of the available functions for a matrix value, you can either call the relevant method for the Matrix object -``` +```php $grid = [ [16, 3, 2, 13], [ 5, 10, 11, 8], @@ -127,8 +137,8 @@ $matrix = new Matrix\Matrix($grid); echo $matrix->trace(); ``` -or you can call the function as you would in procedural code, passing the Matrix object as an argument -``` +or you can call the function as you would in procedural code, passing the Matrix object as an argument +```php $grid = [ [16, 3, 2, 13], [ 5, 10, 11, 8], @@ -140,7 +150,7 @@ $matrix = new Matrix\Matrix($grid); echo Matrix\trace($matrix); ``` When called procedurally using the function, you can pass in the argument as a Matrix object, or as an array. -``` +```php $grid = [ [16, 3, 2, 13], [ 5, 10, 11, 8], @@ -151,7 +161,7 @@ $grid = [ echo Matrix\trace($grid); ``` As an alternative, it is also possible to call the method directly from the `Functions` class. -``` +```php $grid = [ [16, 3, 2, 13], [ 5, 10, 11, 8], @@ -163,3 +173,33 @@ $matrix = new Matrix\Matrix($grid); echo Matrix\Functions::trace($matrix); ``` Used this way, methods must be called statically, and the argument must be the Matrix object, and cannot be an array. + +## Decomposition + +The library also provides classes for matrix decomposition. You can access these using +```php +$grid = [ + [1, 2], + [3, 4], +]; + +$matrix = new Matrix\Matrix($grid); + +$decomposition = new Matrix\Decomposition\QR($matrix); +$Q = $decomposition->getQ(); +$R = $decomposition->getR(); +``` + +or alternatively us the `Decomposition` factory, identifying which form of decomposition you want to use +```php +$grid = [ + [1, 2], + [3, 4], +]; + +$matrix = new Matrix\Matrix($grid); + +$decomposition = Matrix\Decomposition\Decomposition::decomposition(Matrix\Decomposition\Decomposition::QR, $matrix); +$Q = $decomposition->getQ(); +$R = $decomposition->getR(); +``` diff --git a/classes/src/Decomposition/Decomposition.php b/classes/src/Decomposition/Decomposition.php new file mode 100644 index 0000000..f014448 --- /dev/null +++ b/classes/src/Decomposition/Decomposition.php @@ -0,0 +1,27 @@ +luMatrix = $matrix->toArray(); + $this->rows = $matrix->rows; + $this->columns = $matrix->columns; + + $this->buildPivot(); + } + + /** + * Get lower triangular factor. + * + * @return Matrix Lower triangular factor + */ + public function getL() + { + $lower = []; + + $columns = min($this->rows, $this->columns); + for ($row = 0; $row < $this->rows; ++$row) { + for ($column = 0; $column < $columns; ++$column) { + if ($row > $column) { + $lower[$row][$column] = $this->luMatrix[$row][$column]; + } elseif ($row === $column) { + $lower[$row][$column] = 1.0; + } else { + $lower[$row][$column] = 0.0; + } + } + } + + return new Matrix($lower); + } + + /** + * Get upper triangular factor. + * + * @return Matrix Upper triangular factor + */ + public function getU() + { + $upper = []; + + $rows = min($this->rows, $this->columns); + for ($row = 0; $row < $rows; ++$row) { + for ($column = 0; $column < $this->columns; ++$column) { + if ($row <= $column) { + $upper[$row][$column] = $this->luMatrix[$row][$column]; + } else { + $upper[$row][$column] = 0.0; + } + } + } + + return new Matrix($upper); + } + + /** + * Return pivot permutation vector. + * + * @return Matrix Pivot matrix + */ + public function getP() + { + $pMatrix = []; + + $pivots = $this->pivot; + $pivotCount = count($pivots); + foreach ($pivots as $row => $pivot) { + $pMatrix[$row] = array_fill(0, $pivotCount, 0); + $pMatrix[$row][$pivot] = 1; + } + + return new Matrix($pMatrix); + } + + /** + * Return pivot permutation vector. + * + * @return array Pivot vector + */ + public function getPivot() + { + return $this->pivot; + } + + /** + * Is the matrix nonsingular? + * + * @return bool true if U, and hence A, is nonsingular + */ + public function isNonsingular() + { + for ($diagonal = 0; $diagonal < $this->columns; ++$diagonal) { + if ($this->luMatrix[$diagonal][$diagonal] === 0.0) { + return false; + } + } + + return true; + } + + private function buildPivot() + { + for ($row = 0; $row < $this->rows; ++$row) { + $this->pivot[$row] = $row; + } + + for ($column = 0; $column < $this->columns; ++$column) { + $luColumn = $this->localisedReferenceColumn($column); + + $this->applyTransformations($column, $luColumn); + + $pivot = $this->findPivot($column, $luColumn); + if ($pivot !== $column) { + $this->pivotExchange($pivot, $column); + } + + $this->computeMultipliers($column); + + unset($luColumn); + } + } + + private function localisedReferenceColumn($column) + { + $luColumn = []; + + for ($row = 0; $row < $this->rows; ++$row) { + $luColumn[$row] = &$this->luMatrix[$row][$column]; + } + + return $luColumn; + } + + private function applyTransformations($column, array $luColumn) + { + for ($row = 0; $row < $this->rows; ++$row) { + $luRow = $this->luMatrix[$row]; + // Most of the time is spent in the following dot product. + $kmax = min($row, $column); + $sValue = 0.0; + for ($kValue = 0; $kValue < $kmax; ++$kValue) { + $sValue += $luRow[$kValue] * $luColumn[$kValue]; + } + $luRow[$column] = $luColumn[$row] -= $sValue; + } + } + + private function findPivot($column, array $luColumn) + { + $pivot = $column; + for ($row = $column + 1; $row < $this->rows; ++$row) { + if (abs($luColumn[$row]) > abs($luColumn[$pivot])) { + $pivot = $row; + } + } + + return $pivot; + } + + private function pivotExchange($pivot, $column) + { + for ($kValue = 0; $kValue < $this->columns; ++$kValue) { + $tValue = $this->luMatrix[$pivot][$kValue]; + $this->luMatrix[$pivot][$kValue] = $this->luMatrix[$column][$kValue]; + $this->luMatrix[$column][$kValue] = $tValue; + } + + $lValue = $this->pivot[$pivot]; + $this->pivot[$pivot] = $this->pivot[$column]; + $this->pivot[$column] = $lValue; + } + + private function computeMultipliers($diagonal) + { + if (($diagonal < $this->rows) && ($this->luMatrix[$diagonal][$diagonal] != 0.0)) { + for ($row = $diagonal + 1; $row < $this->rows; ++$row) { + $this->luMatrix[$row][$diagonal] /= $this->luMatrix[$diagonal][$diagonal]; + } + } + } +} diff --git a/classes/src/Decomposition/QR.php b/classes/src/Decomposition/QR.php new file mode 100644 index 0000000..e267506 --- /dev/null +++ b/classes/src/Decomposition/QR.php @@ -0,0 +1,152 @@ +qrMatrix = $matrix->toArray(); + $this->rows = $matrix->rows; + $this->columns = $matrix->columns; + + $this->decompose(); + } + + public function getHouseholdVectors() + { + $householdVectors = []; + for ($row = 0; $row < $this->rows; ++$row) { + for ($column = 0; $column < $this->columns; ++$column) { + if ($row >= $column) { + $householdVectors[$row][$column] = $this->qrMatrix[$row][$column]; + } else { + $householdVectors[$row][$column] = 0.0; + } + } + } + + return new Matrix($householdVectors); + } + + public function getQ() + { + $qGrid = []; + + $rowCount = $this->rows; + for ($k = $this->columns - 1; $k >= 0; --$k) { + for ($i = 0; $i < $this->rows; ++$i) { + $qGrid[$i][$k] = 0.0; + } + $qGrid[$k][$k] = 1.0; + if ($this->columns > $this->rows) { + $qGrid = array_slice($qGrid, 0, $this->rows); + } + + for ($j = $k; $j < $this->columns; ++$j) { + if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) { + $s = 0.0; + for ($i = $k; $i < $this->rows; ++$i) { + $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j]; + } + $s = -$s / $this->qrMatrix[$k][$k]; + for ($i = $k; $i < $this->rows; ++$i) { + $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k]; + } + } + } + } + + array_walk( + $qGrid, + function (&$row) use ($rowCount) { + $row = array_reverse($row); + $row = array_slice($row, 0, $rowCount); + } + ); + + return new Matrix($qGrid); + } + + public function getR() + { + $rGrid = []; + + for ($row = 0; $row < $this->columns; ++$row) { + for ($column = 0; $column < $this->columns; ++$column) { + if ($row < $column) { + $rGrid[$row][$column] = isset($this->qrMatrix[$row][$column]) ? $this->qrMatrix[$row][$column] : 0.0; + } elseif ($row === $column) { + $rGrid[$row][$column] = isset($this->rDiagonal[$row]) ? $this->rDiagonal[$row] : 0.0; + } else { + $rGrid[$row][$column] = 0.0; + } + } + } + + if ($this->columns > $this->rows) { + $rGrid = array_slice($rGrid, 0, $this->rows); + } + + return new Matrix($rGrid); + } + + private function hypo($a, $b) + { + if (abs($a) > abs($b)) { + $r = $b / $a; + $r = abs($a) * sqrt(1 + $r * $r); + } elseif ($b != 0.0) { + $r = $a / $b; + $r = abs($b) * sqrt(1 + $r * $r); + } else { + $r = 0.0; + } + + return $r; + } + + /** + * QR Decomposition computed by Householder reflections. + */ + private function decompose() + { + for ($k = 0; $k < $this->columns; ++$k) { + // Compute 2-norm of k-th column without under/overflow. + $norm = 0.0; + for ($i = $k; $i < $this->rows; ++$i) { + $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]); + } + if ($norm != 0.0) { + // Form k-th Householder vector. + if ($this->qrMatrix[$k][$k] < 0.0) { + $norm = -$norm; + } + for ($i = $k; $i < $this->rows; ++$i) { + $this->qrMatrix[$i][$k] /= $norm; + } + $this->qrMatrix[$k][$k] += 1.0; + // Apply transformation to remaining columns. + for ($j = $k + 1; $j < $this->columns; ++$j) { + $s = 0.0; + for ($i = $k; $i < $this->rows; ++$i) { + $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j]; + } + $s = -$s / $this->qrMatrix[$k][$k]; + for ($i = $k; $i < $this->rows; ++$i) { + $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k]; + } + } + } + $this->rDiagonal[$k] = -$norm; + } + } +} diff --git a/unitTests/classes/src/Decomposition/LuTest.php b/unitTests/classes/src/Decomposition/LuTest.php new file mode 100644 index 0000000..92328ba --- /dev/null +++ b/unitTests/classes/src/Decomposition/LuTest.php @@ -0,0 +1,426 @@ +assertOriginalMatrixIsUnchanged($grid, $matrix); + } + + /** + * @dataProvider dataProvider + */ + public function testLUDecompositionL($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new LU($matrix); + + $L = $decomposition->getL(); + // Must return an object of the correct type... + $this->assertIsMatrixObject($L); + // ... containing the correct data + $this->assertMatrixValues($L, count($expected['L']), count($expected['L'][0]), $expected['L']); + } + + /** + * @dataProvider dataProvider + */ + public function testLUDecompositionU($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new LU($matrix); + + $U = $decomposition->getU(); + // Must return an object of the correct type... + $this->assertIsMatrixObject($U); + // ... containing the correct data + $this->assertMatrixValues($U, count($expected['U']), count($expected['U'][0]), $expected['U']); + } + + /** + * @dataProvider dataProvider + */ + public function testLUDecompositionPivots($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new LU($matrix); + + $pivots = $decomposition->getPivot(); + $this->assertEquals($expected['pivots'], $pivots); + } + + /** + * @dataProvider dataProvider + */ + public function testLUDecompositionP($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new LU($matrix); + + $P = $decomposition->getP(); + // Must return an object of the correct type... + $this->assertIsMatrixObject($P); + // ... containing the correct data + $this->assertMatrixValues($P, count($expected['P']), count($expected['P'][0]), $expected['P']); + } + + /** + * @dataProvider dataProvider + */ + public function testLUDecompositionResolve($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new LU($matrix); + + $L = $decomposition->getL(); + $U = $decomposition->getU(); + $P = $decomposition->getP(); + + $result1 = $P->inverse()->multiply($L->multiply($U)); + $this->assertEquals($matrix, $result1); + + $matrix2 = $P->multiply($matrix); + $result2 = $L->multiply($U); + $this->assertEquals($matrix2, $result2); + } + + public function dataProvider() + { + return [ + 'Simple 2x2' => [ + [ + 'U' => [ + [3, 4], + [0, 0.666666666667], + ], + 'L' => [ + [1, 0], + [0.333333333333, 1], + ], + 'P' => [ + [0, 1], + [1, 0], + ], + 'pivots' => [1, 0], + ], + [ + [1, 2], + [3, 4], + ], + ], + 'All the Ones' => [ + [ + 'U' => [ + [1, 1, 1, 1], + [0, -2, -2, 0], + [0, 0, -2, -2], + [0, 0, 0, -4], + ], + 'L' => [ + [1, 0, 0, 0], + [1, 1, 0, 0], + [1, 0, 1, 0], + [1, 1, -1, 1], + ], + 'P' => [ + [1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + ], + 'pivots' => [0, 2, 1, 3], + ], + [ + [1, 1, 1, 1], + [1, 1, -1, -1], + [1, -1, -1, 1], + [1, -1, 1, -1], + ], + ], + 'Simple 3x3 Matrix' => [ + [ + 'U' => [ + [7, 8, 9], + [0, 0.85714285714286, 1.7142857142857], + [0, 0, 5e-15], + ], + 'L' => [ + [1, 0, 0], + [0.14285714285714, 1, 0], + [0.57142857142857, 0.5, 1], + ], + 'P' => [ + [0, 0, 1], + [1, 0, 0], + [0, 1, 0], + ], + 'pivots' => [2, 0, 1], + ], + [ + [1, 2, 3], + [4, 5, 6], + [7, 8, 9], + ], + ], + '3x3 Magic Square' => [ + [ + 'U' => [ + [9, 5, 1], + [0, 5.888888888888890, 5.777777777777780], + [0, 0, 6.792452830188680], + ], + 'L' => [ + [1, 0, 0], + [0.222222222222222, 1, 0], + [0.444444444444444, 0.132075471698113, 1], + ], + 'P' => [ + [0, 1, 0], + [1, 0, 0], + [0, 0, 1], + ], + 'pivots' => [1, 0, 2], + ], + [ + [2, 7, 6], + [9, 5, 1], + [4, 3, 8], + ], + ], + 'Another Simple 3x3' => [ + [ + 'U' => [ + [2, 4, 7], + [0, 1, 1.5], + [0, 0, -2], + ], + 'L' => [ + [1, 0, 0], + [0.5, 1, 0], + [0.5, -1, 1], + ], + 'P' => [ + [0, 1, 0], + [1, 0, 0], + [0, 0, 1], + ], + 'pivots' => [1, 0, 2], + ], + [ + [1, 3, 5], + [2, 4, 7], + [1, 1, 0], + ], + ], + '4x4 with positive and negative values' => [ + [ + 'U' => [ + [5, -4, -3, 1], + [0, 5.6, 2.2, -3.4], + [0, 0, -2.6428571428571, 7.3571428571429], + [0, 0, 0, 9.4594594594595], + ], + 'L' => [ + [1, 0, 0, 0], + [0.4, 1, 0, 0], + [0.8, 0.92857142857143, 1, 0], + [-0.2, -0.5, -0.94594594594595, 1], + ], + 'P' => [ + [0, 0, 0, 1], + [1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + ], + 'pivots' => [3, 0, 2, 1], + ], + [ + [2, 4, 1, -3], + [-1, -2, 2, 4], + [4, 2, -3, 5], + [5, -4, -3, 1], + ], + ], + '4x4 with positive and negative float values' => [ + [ + 'U' => [ + [-2.5, -1, 0, 1], + [0, 2.1, 5, 12.4], + [0, 0, 0.11904761904762, -0.9047619047619], + [0, 0, 0, 0.10000000000002], + ], + 'L' => [ + [1, 0, 0, 0], + [-0.4, 1, 0, 0], + [0, 0.47619047619048, 1, 0], + [0.4, 0.19047619047619, 0.40000000000002, 1], + ], + 'P' => [ + [1, 0, 0, 0], + [0, 0, 0, 1], + [0, 0, 1, 0], + [0, 1, 0, 0], + ], + 'pivots' => [0, 3, 2, 1], + ], + [ + [-2.5, -1, 0, 1], + [-1, 0, 1, 2.5], + [0, 1, 2.5, 5], + [1, 2.5, 5, 12], + ], + ], + '4x4 Magic square' =>[ + [ + 'U' => [ + [14, 1, 8, 11], + [0, 14.714285714286, 7.714285714286, 1.857142857143], + [0, 0, -4.951456310680, 8.252427184466], + [0, 0, 0, 0], + ], + 'L' => [ + [1, 0, 0, 0], + [0.285714285714, 1, 0, 0], + [0.642857142857, 0.364077669903, 1, 0], + [0.5, 0.781553398058, -0.6, 1], + ], + 'P' => [ + [0, 0, 1, 0], + [0, 1, 0, 0], + [1, 0, 0, 0], + [0, 0, 0, 1], + ], + 'pivots' => [2, 1, 0, 3], + ], + [ + [9, 6, 3, 16], + [4, 15, 10, 5], + [14, 1, 8, 11], + [7, 12, 13, 2], + ], + ], + 'Asymetric - 2 rows 3 columns' => [ + [ + 'U' => [ + [4, 5, 6], + [0, 0.75, 1.5], + ], + 'L' => [ + [1, 0], + [0.25, 1], + ], + 'P' => [ + [0, 1], + [1, 0], + ], + 'pivots' => [1, 0], + ], + [ + [1, 2, 3], + [4, 5, 6], + ], + ], + 'Asymetric - 3 rows 2 columns' => [ + [ + 'U' => [ + [5, 6], + [0, 0.8], + ], + 'L' => [ + [1, 0], + [0.2, 1], + [0.6, 0.5], + ], + 'P' => [ + [0, 0, 1], + [1, 0, 0], + [0, 1, 0], + ], + 'pivots' => [2, 0, 1], + ], + [ + [1, 2], + [3, 4], + [5, 6], + ], + ], + 'Asymetric - 3 rows 7 columns' => [ + [ + 'U' => [ + [9, 8, 7, 6, 5, 4, 3], + [0, 1.111111111111, 2.222222222222, 3.333333333333, 4.444444444444, 5.555555555556, 6.666666666667], + [0, 0, 0, 0, 0, 0, 0], + ], + 'L' => [ + [1, 0, 0], + [0.111111111111, 1, 0], + [-0.222222222222, 0.700, 1.000], + ], + 'P' => [ + [0, 1, 0], + [1, 0, 0], + [0, 0, 1], + ], + 'pivots' => [1, 0, 2], + ], + [ + [1, 2, 3, 4, 5, 6, 7], + [9, 8, 7, 6, 5, 4, 3], + [-2, -1, 0, 1, 2, 3, 4], + ], + ], + 'Asymetric - 7 rows 3 columns' => [ + [ + 'U' => [ + [7, 9, 8], + [0, -3, -3], + [0, 0, -4.428571428571429], + ], + 'L' => [ + [1, 0, 0], + [1, 1, 0], + [0.5714285714285714, 0.7142857142857141, 1], + [0.5714285714285714, 0.04761904761904745, -0.35483870967741926], + [0.14285714285714285, -0.23809523809523814, -0.25806451612903225], + [-0.14285714285714285, -0.42857142857142855, -0.1935483870967742], + [0.2857142857142857, -0.14285714285714293, -0.29032258064516125], + ], + 'P' => [ + [0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 1, 0, 0], + [0, 1, 0, 0, 0, 0, 0], + [1, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 0, 1], + ], + 'pivots' => [2, 3, 4, 1, 0, 5, 6], + ], + [ + [1, 2, 3], + [4, 5, 6], + [7, 9, 8], + [7, 6, 5], + [4, 3, -2], + [-1, 0, 1], + [2, 3, 4], + ], + ], + ]; + } +} diff --git a/unitTests/classes/src/Decomposition/QrTest.php b/unitTests/classes/src/Decomposition/QrTest.php new file mode 100644 index 0000000..b8059d8 --- /dev/null +++ b/unitTests/classes/src/Decomposition/QrTest.php @@ -0,0 +1,318 @@ +assertOriginalMatrixIsUnchanged($grid, $matrix); + } + + /** + * @dataProvider dataProvider + */ + public function testQRDecompositionQ($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new QR($matrix); + + $Q = $decomposition->getQ(); + // Must return an object of the correct type... + $this->assertIsMatrixObject($Q); + // ... containing the correct data + $this->assertMatrixValues($Q, count($expected['Q']), count($expected['Q'][0]), $expected['Q']); + } + + /** + * @dataProvider dataProvider + */ + public function testQRDecompositionR($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new QR($matrix); + + $R = $decomposition->getR(); + // Must return an object of the correct type... + $this->assertIsMatrixObject($R); + // ... containing the correct data + $this->assertMatrixValues($R, count($expected['R']), count($expected['R'][0]), $expected['R']); + } + + /** + * @dataProvider dataProvider + */ + public function testQRDecompositionResolve($expected, $grid) + { + $matrix = new Matrix($grid); + $decomposition = new QR($matrix); + + $Q = $decomposition->getQ(); + $R = $decomposition->getR(); + $result = $Q->multiply($R); + + $this->assertEquals($matrix, $result); + } + + public function dataProvider() + { + return [ + 'Simple 2x2' => [ + [ + 'Q' => [ + [-0.316227766016838, 0.948683298050514], + [-0.9486832980505138, -0.316227766016838], + ], + 'R' => [ + [-3.162277660168, -4.427188724236], + [0.0, 0.6324555320336751], + ], + ], + [ + [1.0, 2.0], + [3.0, 4.0], + ], + ], + 'All the Ones' => [ + [ + 'Q' => [ + [-0.5, -0.5, 0.5, 0.5], + [-0.5, -0.5, -0.5, -0.5], + [-0.5, 0.5, -0.5, 0.5], + [-0.5, 0.5, 0.5, -0.5], + ], + 'R' => [ + [-2.0, 0.0, 0.0, 0.0], + [0.0, -2.0, 0.0, 0.0], + [0.0, 0.0, 2.0, 0.0], + [0.0, 0.0, 0.0, 2.0], + ], + ], + [ + [1, 1, 1, 1], + [1, 1, -1, -1], + [1, -1, -1, 1], + [1, -1, 1, -1], + ], + ], + 'Simple 3x3 Matrix' => [ + [ + 'Q' => [ + [-0.12309149097933281, 0.9045340337332908, -0.40824829046386363], + [-0.4923659639173309, 0.30151134457776396, 0.8164965809277259], + [-0.8616404368553291, -0.30151134457776396, -0.4082482904638628], + ], + 'R' => [ + [-8.124038404635960, -9.601136296387950, -11.078234188139900], + [0.0, 0.904534033733293, 1.809068067466580], + [0.0, 0.0, -0.000000000000001], + ], + ], + [ + [1, 2, 3], + [4, 5, 6], + [7, 8, 9], + ], + ], + '3x3 Magic Square' => [ + [ + 'Q' => [ + [-0.199007438041998, 0.972488642858958, 0.1210862465943296], + [-0.895533471188990, -0.230643324147080, 0.3805567750107504], + [-0.398014876083996, 0.032703157901452, -0.9167958670713532], + ], + 'R' => [ + [-10.049875621120900, -7.064764050490930, -5.273697108112940], + [0.0, 5.752313352981660, 5.865913796218280], + [0.0, 0.0, -6.227292681994097], + ], + ], + [ + [2, 7, 6], + [9, 5, 1], + [4, 3, 8], + ], + ], + 'Another Simple 3x3' => [ + [ + 'Q' => [ + [-0.408248290463863, 0.707106781186548, -0.5773502691896256], + [-0.816496580927726, 0.0, 0.5773502691896257], + [-0.408248290463863, -0.707106781186547, -0.5773502691896262], + ], + 'R' => [ + [-2.449489742783180, -4.898979485566360, -7.756717518813400], + [0.0, 1.414213562373100, 3.535533905932740], + [0.0, 0.0, 1.1547005383792528], + ], + ], + [ + [1, 3, 5], + [2, 4, 7], + [1, 1, 0], + ], + ], + '4x4 with positive and negative values' => [ + [ + 'Q' => [ + [-0.294883912309794, 0.646908108895345, 0.655396891228784, -0.25496723686380474], + [0.147441956154897, -0.323454054447673, 0.647746343860355, 0.6738419831400553], + [-0.589767824619589, 0.344100057923056, -0.357025543860038, 0.637418092159512], + [-0.737209780774486, -0.598734100786117, 0.153010947368588, -0.27317918235407657], + ], + 'R' => [ + [-6.782329983125270, 0.294883912309793, 3.980932816182220, -2.211629342323460], + [0.0, 6.317677063467310, 0.763902128589184, -2.112774355647560], + [0.0, 0.0, 2.562933368423840, -1.007322070176540], + [0.0, 0.0, 0.0, 6.374180921595119], + ], + ], + [ + [2, 4, 1, -3], + [-1, -2, 2, 4], + [4, 2, -3, 5], + [5, -4, -3, 1], + ], + ], + '4x4 with positive and negative float values' => [ + [ + 'Q' => [ + [-0.870388279778489, -0.225482237707286, 0.265273220832105, -0.34815531191139854], + [-0.348155311911396, -0.265273220832101, -0.225482237707296, 0.8703882797784915], + [0.0, -0.437700814372966, -0.828978815100310, -0.34815531191138704], + [0.348155311911396, -0.828978815100314, 0.437700814372966, -0.000000000000005], + ], + 'R' => [ + [2.872281323269010, 1.740776559556980, 1.392621247645580, 2.437087183379770], + [0.0, -2.284665614416470, -5.504419332266090, -13.024915142856100], + [0.0, 0.0, -0.109425203593242, 0.809083323537904], + [0.0, 0.0, 0.0, 0.0870388279778368], + ], + ], + [ + [-2.5, -1, 0, 1], + [-1, 0, 1, 2.5], + [0, 1, 2.5, 5], + [1, 2.5, 5, 12], + ], + ], + '4x4 Magic square' =>[ + [ + 'Q' => [ + [-0.486664263392288, -0.025409618431111, 0.793825739600044, -0.36380343755449956], + [-0.216295228174350, -0.755583236958464, 0.121092061972888, 0.6063390625908325], + [-0.757033298610225, 0.463372625000684, -0.282548144603405, 0.3638034375544995], + [-0.378516649305113, -0.462313890899387, -0.524732268549182, -0.6063390625908323], + ], + 'R' => [ + [-18.493242008906900, -11.463647093240600, -14.599927901768600, -17.952503938471100], + [0.0, -16.570600330755600, -9.935160806564540, -0.011998986481357], + [0.0, 0.0, -5.489506809437590, 9.149178015729320], + [0.0, 0.0, 0.0, 0.0], + ], + ], + [ + [9, 6, 3, 16], + [4, 15, 10, 5], + [14, 1, 8, 11], + [7, 12, 13, 2], + ], + ], + 'Asymetric - 2 rows 3 columns' => [ + [ + 'Q' => [ + [-0.24253562503633308, 0.9701425001453319], + [-0.9701425001453319, -0.24253562503633308], + ], + 'R' => [ + [-4.123105625617661, -5.335783750799326, -6.5484618759809905], + [0.0, 0.7276068751089992, 1.4552137502179976], + ], + ], + [ + [1, 2, 3], + [4, 5, 6], + ], + ], + 'Asymetric - 3 rows 2 columns' => [ + [ + 'Q' => [ + [-0.169030850945703, 0.897085227145060], + [-0.507092552837110, 0.276026223736942], + [-0.845154254728517, -0.345032779671177], + ], + 'R' => [ + [-5.916079783099620, -7.437357441610940], + [0.0, 0.828078671210823], + ], + ], + [ + [1, 2], + [3, 4], + [5, 6], + ], + ], + 'Asymetric - 3 rows 7 columns' => [ + [ + 'Q' => [ + [-0.10783277320343831, 0.8235566226707017, 0.5568900989230126], + [-0.9704949588309455, 0.0343148592779462, -0.23866718525271913], + [0.2156655464068768, 0.56619517808611, -0.7955572841757288], + ], + 'R' => [ + // phpcs:disable Generic.Files.LineLength + [-9.273618495495706, -8.195290763461317, -7.116963031426934, -6.038635299392549, -4.960307567358166, -3.881979835323781, -2.803652103289398], + [0.0, 1.355436941478862, 2.7108738829577295, 4.066310824436593, 5.421747765915459, 6.7771847073943245, 8.13262164887319], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + // phpcs:enable + ], + ], + [ + [1, 2, 3, 4, 5, 6, 7], + [9, 8, 7, 6, 5, 4, 3], + [-2, -1, 0, 1, 2, 3, 4], + ], + ], + 'Asymetric - 7 rows 3 columns' => [ + [ + 'Q' => [ + [-0.085749292571255, -0.3045082017063, 0.058918029340522], + [-0.34299717028502, -0.23977023756402, 0.27272366158654], + [-0.60024504799878, -0.5011197965088, -0.30268884489375], + [-0.60024504799878, 0.4771427727524, 0.51123032331553], + [-0.34299717028502, 0.41240480861011, -0.73839880757692], + [0.085749292571254, -0.34766684446783, -0.08361905882349], + [-0.17149858514251, -0.28292888032554, 0.13018657342253], + ], + 'R' => [ + [-11.661903789690601, -12.433647422831891, -10.032667230836765], + [0.0, -3.0666613384437933, -6.27958252180164], + [0.0, 0.0, 3.8616617649889395], + ], + ], + [ + [1, 2, 3], + [4, 5, 6], + [7, 9, 8], + [7, 6, 5], + [4, 3, -2], + [-1, 0, 1], + [2, 3, 4], + ], + ], + ]; + } +}