Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add QR algorithm #472

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Add QR algorithm #472

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
cd45a97
Replace array_search with floating-point filter
Aweptimum Nov 6, 2023
7c75688
Make 8.2 linter happy
Aweptimum Nov 6, 2023
22ffbbf
Add NumericaMatrix->upperHessenberg
Aweptimum Oct 31, 2023
9fb768a
Test uppserHessenberg
Aweptimum Oct 31, 2023
44637c3
Add qrAlgorithm
Aweptimum Oct 30, 2023
51ba07e
test QRAlgorithm eigenvalues
Aweptimum Oct 31, 2023
edb1bb3
Shift QR by Wilkinson, refactor a bit
Aweptimum Oct 31, 2023
51a36e9
Add qr method to Eigenvector
Aweptimum Nov 1, 2023
6b445d9
Test Eigenvector::qrAlgorithm
Aweptimum Nov 1, 2023
2f08e84
Add failing eigenvector cases to eigenvalues
Aweptimum Nov 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 69 additions & 0 deletions src/LinearAlgebra/Eigenvalue.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,20 +2,24 @@

namespace MathPHP\LinearAlgebra;

use MathPHP\Arithmetic;
use MathPHP\Exception;
use MathPHP\Expression\Polynomial;
use MathPHP\Functions\Support;
use MathPHP\LinearAlgebra\Decomposition\QR;

class Eigenvalue
{
public const CLOSED_FORM_POLYNOMIAL_ROOT_METHOD = 'closedFormPolynomialRootMethod';
public const POWER_ITERATION = 'powerIteration';
public const JACOBI_METHOD = 'jacobiMethod';
public const QR_ALGORITHM = 'qrAlgorithm';

private const METHODS = [
self::CLOSED_FORM_POLYNOMIAL_ROOT_METHOD,
self::POWER_ITERATION,
self::JACOBI_METHOD,
self::QR_ALGORITHM,
];

/**
Expand Down Expand Up @@ -261,4 +265,69 @@

return [$max_ev];
}

public static function qrAlgorithm(NumericMatrix $A): array

Check failure on line 269 in src/LinearAlgebra/Eigenvalue.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Method MathPHP\LinearAlgebra\Eigenvalue::qrAlgorithm() return type has no value type specified in iterable type array.
{
self::checkMatrix($A);

if ($A->isUpperHessenberg()) {
$H = $A;
} else {
$H = $A->upperHessenberg();
}

return self::qrIteration($H);
}

private static function qrIteration(NumericMatrix $A, array &$values = null): array

Check failure on line 282 in src/LinearAlgebra/Eigenvalue.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Method MathPHP\LinearAlgebra\Eigenvalue::qrIteration() has parameter $values with no value type specified in iterable type array.

Check failure on line 282 in src/LinearAlgebra/Eigenvalue.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Method MathPHP\LinearAlgebra\Eigenvalue::qrIteration() return type has no value type specified in iterable type array.
{
$values = $values ?? [];
$e = $A->getError();
$n = $A->getM();

if ($A->getM() === 1) {
$values[] = $A[0][0];
return $values;
}

$ident = MatrixFactory::identity($n);

do {
// Pick shift
$shift = self::calcShift($A);
$A = $A->subtract($ident->scalarMultiply($shift));

// decompose
$qr = QR::decompose($A);
$A = $qr->R->multiply($qr->Q);

// shift back
$A = $A->add($ident->scalarMultiply($shift));
} while (!Arithmetic::almostEqual($A[$n-1][$n-2], 0, $e)); // subdiagonal entry needs to be 0

// Check if we can deflate
$eig = $A[$n-1][$n-1];
$values[] = $eig;
self::qrIteration($A->submatrix(0, 0, $n-2, $n-2), $values);

return array_reverse($values);
}

private static function calcShift(NumericMatrix $A): float
{
$n = $A->getM() - 1;
$sigma = .5 * ($A[$n-1][$n-1] - $A[$n][$n]);
$sign = $sigma >= 0 ? 1 : -1;

$numerator = ($sign * ($A[$n][$n-1])**2);
$denominator = abs($sigma) + sqrt($sigma**2 + ($A[$n-1][$n-1])**2);

try {
$fraction = $numerator / $denominator;
} catch (\DivisionByZeroError $error) {

Check failure on line 327 in src/LinearAlgebra/Eigenvalue.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Dead catch - DivisionByZeroError is never thrown in the try block.
$fraction = 0;
}

return $A[$n][$n] - $fraction;
}
}
18 changes: 16 additions & 2 deletions src/LinearAlgebra/Eigenvector.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,13 +2,21 @@

namespace MathPHP\LinearAlgebra;

use MathPHP\Arithmetic;
use MathPHP\Exception;
use MathPHP\Exception\MatrixException;
use MathPHP\Functions\Map\Single;
use MathPHP\Functions\Special;

class Eigenvector
{
public static function qrAlgorithm(NumericMatrix $A)

Check failure on line 13 in src/LinearAlgebra/Eigenvector.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Method MathPHP\LinearAlgebra\Eigenvector::qrAlgorithm() has no return type specified.
{
$eigenvalues = Eigenvalue::qrAlgorithm($A);

return self::eigenvectors($A, $eigenvalues);
}

/**
* Calculate the Eigenvectors for a matrix
*
Expand Down Expand Up @@ -66,8 +74,14 @@
foreach ($eigenvalues as $eigenvalue) {
// If this is a duplicate eigenvalue, and this is the second instance, the first
// pass already found all the vectors.
$key = \array_search($eigenvalue, \array_column($solution_array, 'eigenvalue'));
if (!$key) {
$key = false;
foreach (\array_column($solution_array, 'eigenvalue') as $i => $v) {
if (Arithmetic::almostEqual($v, $eigenvalue, $A->getError())) {
$key = $i;
break;
}
}
if ($key === false) {
$Iλ = MatrixFactory::identity($number)->scalarMultiply($eigenvalue);
$T = $A->subtract($Iλ);

Expand Down
60 changes: 60 additions & 0 deletions src/LinearAlgebra/NumericMatrix.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

namespace MathPHP\LinearAlgebra;

use MathPHP\Arithmetic;
use MathPHP\Functions\Map;
use MathPHP\Functions\Support;
use MathPHP\Exception;
Expand Down Expand Up @@ -1478,6 +1479,7 @@
* - covarianceMatrix
* - adjugate
* - householder
* - upperHessenberg
**************************************************************************/

/**
Expand Down Expand Up @@ -1886,6 +1888,64 @@
return Householder::transform($this);
}

/**
* Uses householder to convert the matrix to upper hessenberg form
*
* @return NumericMatrix
*
* @throws Exception\MathException
*/
public function upperHessenberg(): NumericMatrix
{
$n = $this->getM();

$hessenberg = $this;
$identity = MatrixFactory::identity($n);

for ($i = 0; $i < $n-2; $i++)
{
$slice = $this->submatrix($i+1, $i, $n-1, $i);
$x = $slice->asVectors()[0];
$sign = $x[0] >= 0 ? 1 : -1;
$x1 = $sign * $x->l2norm();
if (Arithmetic::almostEqual($x1, 0, $this->getError())) {
continue;
}
$u = \array_fill(0, $n-$i-1, 0);
$u[0] = $x1;

Check failure on line 1915 in src/LinearAlgebra/NumericMatrix.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Cannot access offset 0 on array<int, 0>|false.
$u = new Vector($u);

$v = $u->subtract($x);
$v = MatrixFactory::createFromVectors([$v]);
$vt = $v->transpose();

$vvt = $v->multiply($vt);
$vtv = $vt->multiply($v)[0][0];
$P = $vvt->scalarDivide($vtv);

$H = MatrixFactory::identity($P->getM());
$H = $H->subtract($P->scalarMultiply(2));

// augment
$offset = $n - $H->getM();
$elems = $identity->getMatrix();
for ($i = 0; $i < $H->getM(); $i++)
{
for ($j = 0; $j < $H->getN(); $j++)
{
$elems[$i + $offset][$j + $offset] = $H[$i][$j];
}
}

$H = MatrixFactory::create($elems);

// Multiply
$hessenberg = $H->multiply($hessenberg->multiply($H));

Check failure on line 1943 in src/LinearAlgebra/NumericMatrix.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Call to an undefined method MathPHP\LinearAlgebra\Matrix<float|int>|MathPHP\LinearAlgebra\ObjectMatrix::multiply().
}

return $hessenberg;

Check failure on line 1946 in src/LinearAlgebra/NumericMatrix.php

View workflow job for this annotation

GitHub Actions / Static Analysis (7.2) 

Method MathPHP\LinearAlgebra\NumericMatrix::upperHessenberg() should return MathPHP\LinearAlgebra\NumericMatrix but returns $this(MathPHP\LinearAlgebra\NumericMatrix)|MathPHP\LinearAlgebra\ObjectMatrix.
}

/**************************************************************************
* MATRIX VECTOR OPERATIONS - Return a Vector
* - vectorMultiply
Expand Down
55 changes: 55 additions & 0 deletions tests/LinearAlgebra/Eigen/EigenvalueTest.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,31 @@ public function testPowerIteration(array $A, array $S, float $max_abs_eigenvalue
$this->assertEqualsWithDelta([$max_abs_eigenvalue], $eigenvalues, 0.0001);
}

/**
* @test qrAlgorithm returns the expected eigenvalues
* @dataProvider dataProviderForEigenvalues
* @dataProvider dataProviderForLargeMatrixEigenvalues
* @dataProvider dataProviderForSymmetricEigenvalues
* @param array $A
* @param array $S
* @param float $max_abs_eigenvalue maximum absolute eigenvalue
* @throws \Exception
*/
public function testQRAlgorithm(array $A, array $S)
{
// Given
$A = MatrixFactory::create($A);

// When
$eigenvalues = Eigenvalue::qrAlgorithm($A);

sort($S);
sort($eigenvalues);

// Then
$this->assertEqualsWithDelta($S, $eigenvalues, 0.0001);
}

/**
* @test Matrix eigenvalues using powerIterationMethod returns the expected eigenvalues
* @dataProvider dataProviderForEigenvalues
Expand Down Expand Up @@ -195,6 +220,24 @@ public function dataProviderForEigenvalues(): array
[2, 2, 1],
2,
],
[
[
[2, 0, 1],
[2, 1, 2],
[3, 0, 4]
],
[5, 1, 1],
5,
],
[ // Matrix has duplicate eigenvalues. no solution on the axis
[
[2, 2, -3],
[2, 5, -6],
[3, 6, -8],
],
[-3, 1, 1],
-3
],
[
[
[1, 2, 1],
Expand Down Expand Up @@ -251,6 +294,18 @@ public function dataProviderForLargeMatrixEigenvalues(): array
],
[4, 3, 2, -2, 1, -1],
4,
],
[ // Failing case
[
[2,0,0,0,0,0],
[0,2,0,0,0,1729.7],
[0,0,2,0,-1729.7,0],
[0,0,0,0,0,0],
[0,0,-1729.7,0,2.82879*10**6,0],
[0,1729.7,0,0,0,2.82879*10**6]
],
[2828791.05765, 0.94235235527, 0.94235235527, 2828791.05765, 2, 0],
2828791.05765
]
];
}
Expand Down
21 changes: 21 additions & 0 deletions tests/LinearAlgebra/Eigen/EigenvectorTest.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,4 +234,25 @@ public function testMatrixEigenvectorInvalidMethodException()
// When
$A->eigenvectors($invalidMethod);
}

/**
* @test qrAlgorithm
* @dataProvider dataProviderForEigenvector
* @param array $A
* @param array $B
*/
public function testQRAlgorithm(array $A, array $S): void
{
// Given
$A = MatrixFactory::create($A);
$S = MatrixFactory::create($S);

// When
$eigenvectors = Eigenvector::qrAlgorithm($A);

// Then
$this->assertEqualsWithDelta($S, $eigenvectors, 0.0001, sprintf(
"Eigenvectors unequal:\nExpected:\n" . (string) $S . "\nActual:\n" . (string) $eigenvectors
));
}
}
38 changes: 38 additions & 0 deletions tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1485,4 +1485,42 @@ public function dataProviderForRank(): array
],
];
}

/**
* @test upperHessenberg reduction
* @dataProvider dataProviderForUpperHessenberg
* @param array $A
* @throws \Exception
*/
public function testUpperHessenberg(array $A, array $H)
{
// Given
$A = MatrixFactory::create($A);
$H = MatrixFactory::create($H);

// When
$actual = $A->upperHessenberg();

// Then
$this->assertTrue($actual->isUpperHessenberg());
$this->assertEqualsWithDelta($H, $actual, 0.0000001);
}

public static function dataProviderForUpperHessenberg(): \Iterator
{
yield [
'A' => [
[3, 0, 0, 0],
[0, 1, 0, 1],
[0, 0, 2, 0],
[0, 1, 0, 1],
],
'H' => [
[3, 0, 0, 0],
[0, 1, 1, 0],
[0, 1, 1, 0],
[0, 0, 0, 2],
]
Comment on lines +1511 to +1523
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Question about Hessenberg matrices: Is there only one upper Hessenberg matrix for any matrix, or are there multiple possible Hessenberg matrices that satisfy the critiera?

I tried this test case in R and SciPy and it got a slightly different answer.

R

> A <- rbind(c(3, 0, 0, 0), c(0, 1, 0, 1), c(0, 0, 2, 0), c(0, 1, 0, 1))
> hb <- hessenberg(A)
> hb
$H
     [,1]          [,2]          [,3]          [,4]
[1,]    3  0.000000e+00  0.000000e+00  0.000000e+00
[2,]    0  1.000000e+00 -1.000000e+00  2.220446e-16
[3,]    0 -1.000000e+00  1.000000e+00 -5.551115e-16
[4,]    0  2.220446e-16 -4.440892e-16  2.000000e+00

$P
     [,1] [,2]          [,3]          [,4]
[1,]    1    0  0.000000e+00  0.000000e+00
[2,]    0    1  0.000000e+00  0.000000e+00
[3,]    0    0  2.220446e-16 -1.000000e+00
[4,]    0    0 -1.000000e+00  2.220446e-16

Python

In [1]: import numpy as np
In [2]: from scipy.linalg import hessenberg

In [10]: A = np.array([[3, 0, 0, 0], [0, 1, 0, 1], [0, 0, 2, 0], [0, 1, 0, 1]])

In [11]: A
Out[11]:
array([[3, 0, 0, 0],
       [0, 1, 0, 1],
       [0, 0, 2, 0],
       [0, 1, 0, 1]])

In [12]: H = hessenberg(A)

In [13]: H
Out[13]:
array([[ 3.,  0.,  0.,  0.],
       [ 0.,  1., -1.,  0.],
       [ 0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  2.]])

This online calculator however gets the same answer as you have here.

I'd like to help with generating test cases, but I don't know enough about it so I want to make sure I understand how it works before offering help.

Thanks,
Mark

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the values are unique, the signs are not. In addition to the upper hessenberg form, the algorithm also computes a permutation matrix, P, such that A = PHPT. Wolfram returns the hessenberg form with the added context of the permutation matrix, and multiplying it accordingly gives the original matrix.

In the NumericMatrix::upperHessenberg method, the P matrix calculated is exactly the same as R, Python, and Wolfram, it just has different signs:

[1, 0, 0, 0]
[0, 1, 0, 0]
[0, 0, 0, 1]
[0, 0, 1, 0]

Multiplying A = PHPT gives the original matrix as well.

So it might be better to return both matrices for context rather than just the upper hessenberg matrix

Also, I had originally tried the other method of getting the matrix in Upper-Hessenberg form using Givens rotations. It seems more widely used than the householder method, and I think that's because it allows for parallelization + works better on sparse matrices. This particular matrix needed (I think) a planar rotation of pi/2 to zero out the 1 at (4, 1), so the givens rotation block was just:

[cos(π/2), -sin(π/2)]    =    [0, -1]
[sin(π/2),  cos(π/2)]         [1,  0]

All that to say, I wouldn't be surprised if the Givens method introduced negative signs to both P and H

];
}
}
Loading