Skip to content

Commit d5a82cc

Browse files
markrogoyskimarkrogoyski
committed
Add Hessenberge matrix decomposition.
1 parent 20295fa commit d5a82cc

File tree

7 files changed

+3198
-4
lines changed

7 files changed

+3198
-4
lines changed

‎src/LinearAlgebra/Decomposition/Hessenberg.php‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
<?php
2+
3+
namespace MathPHP\LinearAlgebra\Decomposition;
4+
5+
use MathPHP\Exception;
6+
use MathPHP\Functions\Support;
7+
use MathPHP\LinearAlgebra\Householder;
8+
use MathPHP\LinearAlgebra\NumericMatrix;
9+
use MathPHP\LinearAlgebra\MatrixFactory;
10+
11+
/**
12+
* Hessenberg Decomposition
13+
*
14+
* A = QHQ*
15+
*
16+
* Where:
17+
* - Q is an orthogonal matrix
18+
* - H is an upper Hessenberg matrix (zeros below the first subdiagonal)
19+
* - Q* is the conjugate transpose of Q
20+
*
21+
* The Hessenberg decomposition is useful as a preprocessing step for eigenvalue algorithms,
22+
* as it preserves eigenvalues while reducing the matrix to a form that is easier to work with.
23+
*
24+
* @property-read NumericMatrix $Q Orthogonal transformation matrix
25+
* @property-read NumericMatrix $H Upper Hessenberg matrix
26+
*/
27+
class Hessenberg extends Decomposition
28+
{
29+
/** @var NumericMatrix Orthogonal transformation matrix */
30+
private $Q;
31+
32+
/** @var NumericMatrix Upper Hessenberg matrix */
33+
private $H;
34+
35+
/**
36+
* Hessenberg constructor
37+
*
38+
* @param NumericMatrix $Q Orthogonal transformation matrix
39+
* @param NumericMatrix $H Upper Hessenberg matrix
40+
*/
41+
private function __construct(NumericMatrix $Q, NumericMatrix $H)
42+
{
43+
$this->Q = $Q;
44+
$this->H = $H;
45+
}
46+
47+
/**
48+
* Decompose a matrix into Hessenberg form using Householder reflections
49+
* Factory method to create Hessenberg objects.
50+
*
51+
* A = QHQ*
52+
*
53+
* Where:
54+
* - Q is an orthogonal matrix
55+
* - H is an upper Hessenberg matrix (zeros below the first subdiagonal)
56+
* - Q* is the conjugate transpose of Q
57+
*
58+
* Algorithm:
59+
* - Uses Householder reflections to systematically zero out elements below the first subdiagonal
60+
* - Each reflection preserves eigenvalues through similarity transformations
61+
* - Requires n-2 Householder transformations for an n×n matrix
62+
*
63+
* Numerical Stability:
64+
* - Uses consistent error tolerance throughout the decomposition process
65+
* - Householder transformations use the input matrix's error tolerance for consistent behavior
66+
* - The algorithm maintains numerical stability through orthogonal similarity transformations
67+
*
68+
* @param NumericMatrix $A Matrix to decompose
69+
*
70+
* @return Hessenberg
71+
*
72+
* @throws Exception\BadDataException
73+
* @throws Exception\IncorrectTypeException
74+
* @throws Exception\MathException
75+
* @throws Exception\MatrixException
76+
* @throws Exception\OutOfBoundsException
77+
* @throws Exception\VectorException
78+
*/
79+
public static function decompose(NumericMatrix $A): Hessenberg
80+
{
81+
if (!$A->isSquare()) {
82+
throw new Exception\MatrixException('Hessenberg decomposition only works on square matrices');
83+
}
84+
85+
$n = $A->getN();
86+
$ε = $A->getError();
87+
88+
// For 1x1 matrices, return trivial decomposition
89+
if ($n < 2) {
90+
$Q = MatrixFactory::identity(1);
91+
$H = clone $A;
92+
$Q->setError($ε);
93+
$H->setError($ε);
94+
return new Hessenberg($Q, $H);
95+
}
96+
97+
// Start with copies of the input matrix and identity matrix
98+
$H = MatrixFactory::createNumeric($A->getMatrix());
99+
$Q = MatrixFactory::identity($n);
100+
101+
// Apply Householder transformations to columns 0 through n-3
102+
// (we need n-2 transformations total for an n×n matrix)
103+
for ($k = 0; $k < $n - 2; $k++) {
104+
// Extract the submatrix starting from row k+1, column k
105+
$subH = $H->submatrix($k + 1, $k, $n - 1, $k);
106+
107+
// Skip if the column below the diagonal is already effectively zero
108+
if (Support::isZero($subH->frobeniusNorm(), $ε)) {
109+
continue;
110+
}
111+
112+
// Create Householder transformation for this column
113+
$reflector = Householder::transform($subH, $ε);
114+
115+
// Embed the Householder matrix in a full-size identity matrix
116+
$fullReflector = self::embedHouseholderMatrix($reflector, $n, $k + 1);
117+
118+
// Apply the transformation: H = QHQ'
119+
$H = $fullReflector->multiply($H)->multiply($fullReflector->transpose());
120+
121+
// Accumulate the orthogonal transformations: Q = QQ'
122+
$Q = $Q->multiply($fullReflector->transpose());
123+
}
124+
125+
// Propagate error tolerance to result matrices
126+
$Q->setError($ε);
127+
$H->setError($ε);
128+
129+
return new Hessenberg($Q, $H);
130+
}
131+
132+
/**
133+
* Embed a Householder transformation matrix into a full-size identity matrix
134+
*
135+
* The Householder matrix is placed starting at the specified offset in both row and column dimensions.
136+
*
137+
* The implementation used is not the most efficient,
138+
* but it helps us visualize what it is, using matrix operations.
139+
* We prefer clarity and simplicity over performance.
140+
*
141+
* An alternative implementation would copy Householder matrix elements
142+
* to the appropriate position in identity matrix in a nested loop.
143+
* It is more efficient, but harder to visualize what the transformation is.
144+
*
145+
* $I = MatrixFactory::identity($n)->getMatrix();
146+
*
147+
* for ($i = 0; $i < $householder->getM(); $i++) {
148+
* for ($j = 0; $j < $householder->getN(); $j++) {
149+
* $I[$offset + $i][$offset + $j] = $householder[$i][$j];
150+
*. }
151+
* }
152+
* return MatrixFactory::createNumeric($I);
153+
*
154+
* @param NumericMatrix $householder The Householder transformation matrix
155+
* @param int $n The size of the full identity matrix
156+
* @param int $offset The starting position to embed the Householder matrix
157+
*
158+
* @return NumericMatrix Full-size matrix with embedded Householder transformation
159+
*
160+
* @throws Exception\MatrixException
161+
* @throws Exception\IncorrectTypeException
162+
*/
163+
private static function embedHouseholderMatrix(NumericMatrix $householder, int $n, int $offset): NumericMatrix
164+
{
165+
$I = MatrixFactory::identity($offset);
166+
167+
// Create the top block: [ I | 0 ]
168+
$topBlock = $I->augment(MatrixFactory::zero($offset, $n - $offset));
169+
170+
// Create the bottom block: [ 0 | H ]
171+
$bottomBlock = MatrixFactory::zero($n - $offset, $offset)->augment($householder);
172+
173+
// Stack them vertically: [ topBlock ]
174+
// [ bottomBlock ]
175+
return $topBlock->augmentBelow($bottomBlock);
176+
}
177+
178+
/**
179+
* Get Q or H matrix
180+
*
181+
* @param string $name
182+
*
183+
* @return NumericMatrix
184+
*
185+
* @throws Exception\MatrixException
186+
*/
187+
public function __get(string $name): NumericMatrix
188+
{
189+
switch ($name) {
190+
case 'Q':
191+
case 'H':
192+
return $this->$name;
193+
194+
default:
195+
throw new Exception\MatrixException("Hessenberg class does not have a gettable property: $name");
196+
}
197+
}
198+
}

‎src/LinearAlgebra/Householder.php‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
namespace MathPHP\LinearAlgebra;
44

55
use MathPHP\Exception;
6+
use MathPHP\Functions\Support;
67

78
/**
89
* A Householder transformation (also known as a Householder reflection or elementary reflector) is a linear
@@ -22,12 +23,13 @@ class Householder
2223
* https://en.wikipedia.org/wiki/Householder_transformation
2324
*
2425
* @param NumericMatrix $A source Matrix
26+
* @param float $ε optional tolerance for zero checking (defaults to PHP_FLOAT_EPSILON)
2527
*
2628
* @return NumericMatrix
2729
*
2830
* @throws Exception\MathException
2931
*/
30-
public static function transform(NumericMatrix $A): NumericMatrix
32+
public static function transform(NumericMatrix $A, float $ε = \PHP_FLOAT_EPSILON): NumericMatrix
3133
{
3234
$m = $A->getM();
3335
$I = MatrixFactory::identity($m);
@@ -47,7 +49,7 @@ public static function transform(NumericMatrix $A): NumericMatrix
4749
$uᵀ = $u->transpose();
4850
$uᵀu = $uᵀ->multiply($u)->get(0, 0);
4951
$uuᵀ = $u->multiply($uᵀ);
50-
if ($uᵀu == 0) {
52+
if (Support::isZero($uᵀu, $ε)) {
5153
return $I;
5254
}
5355

‎src/LinearAlgebra/MatrixCatalog.php‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,9 @@ class MatrixCatalog
3636
/** @var Decomposition\SVD */
3737
private $SVD;
3838

39+
/** @var Decomposition\Hessenberg */
40+
private $hessenberg;
41+
3942
/** @var int|float|ObjectArithmetic determinant */
4043
private $det;
4144

@@ -166,6 +169,7 @@ public function getReducedRowEchelonForm(): Reduction\ReducedRowEchelonForm
166169
* - Cholesky decomposition
167170
* - Crout decomposition
168171
* - SVD
172+
* - Hessenberg decomposition
169173
**************************************************************************/
170174

171175

@@ -303,6 +307,33 @@ public function getSVD(): Decomposition\SVD
303307
{
304308
return $this->SVD;
305309
}
310+
311+
// HESSENBERG DECOMPOSITION
312+
313+
/**
314+
* @param Decomposition\Hessenberg $hessenberg
315+
*/
316+
public function addHessenbergDecomposition(Decomposition\Hessenberg $hessenberg): void
317+
{
318+
$this->hessenberg = $hessenberg;
319+
}
320+
321+
/**
322+
* @return bool
323+
*/
324+
public function hasHessenbergDecomposition(): bool
325+
{
326+
// @phpstan-ignore-next-line
327+
return isset($this->hessenberg);
328+
}
329+
330+
/**
331+
* @return Decomposition\Hessenberg
332+
*/
333+
public function getHessenbergDecomposition(): Decomposition\Hessenberg
334+
{
335+
return $this->hessenberg;
336+
}
306337
/**************************************************************************
307338
* DERIVED DATA
308339
* - determinant

‎src/LinearAlgebra/NumericMatrix.php‎

Lines changed: 27 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -886,7 +886,7 @@ public function isUpperHessenberg(): bool
886886
// Elements below lower diagonal are zero
887887
for ($i = 2; $i < $this->m; $i++) {
888888
for ($j = 0; $j < $i - 1; $j++) {
889-
if ($this->A[$i][$j] != 0) {
889+
if (Support::isNotZero($this->A[$i][$j], $this->ε)) {
890890
return false;
891891
}
892892
}
@@ -913,7 +913,7 @@ public function isLowerHessenberg(): bool
913913
// Elements above upper diagonal are zero
914914
for ($i = 0; $i < $this->m; $i++) {
915915
for ($j = $i + 2; $j < $this->n; $j++) {
916-
if ($this->A[$i][$j] != 0) {
916+
if (Support::isNotZero($this->A[$i][$j], $this->ε)) {
917917
return false;
918918
}
919919
}
@@ -2741,6 +2741,7 @@ public function rref(): Reduction\ReducedRowEchelonForm
27412741
* - Cholesky decomposition
27422742
* - Crout decomposition
27432743
* - SVD (Singular Value Decomposition)
2744+
* - Hessenberg decomposition
27442745
********************************************************************************/
27452746

27462747
/**
@@ -2864,6 +2865,30 @@ public function svd(): Decomposition\SVD
28642865
return $this->catalog->getSVD();
28652866
}
28662867

2868+
/**
2869+
* Hessenberg Decomposition
2870+
*
2871+
* A = QHQ*
2872+
*
2873+
* Where:
2874+
* Q is an orthogonal matrix
2875+
* H is an upper Hessenberg matrix (zeros below the first subdiagonal)
2876+
* Q* is the conjugate transpose of Q
2877+
*
2878+
* @return Decomposition\Hessenberg
2879+
*
2880+
* @throws Exception\MatrixException if matrix is not square
2881+
* @throws Exception\MathException
2882+
*/
2883+
public function hessenbergDecomposition(): Decomposition\Hessenberg
2884+
{
2885+
if (!$this->catalog->hasHessenbergDecomposition()) {
2886+
$this->catalog->addHessenbergDecomposition(Decomposition\Hessenberg::decompose($this));
2887+
}
2888+
2889+
return $this->catalog->getHessenbergDecomposition();
2890+
}
2891+
28672892
/**************************************************************************
28682893
* SOLVE LINEAR SYSTEM OF EQUATIONS
28692894
* - solve

0 commit comments

Comments
 (0)