php-ai / php-ml

src/Phpml/Math/LinearAlgebra/EigenvalueDecomposition.php

<?php

declare(strict_types=1);
/**
 *
 *	Class to obtain eigenvalues and eigenvectors of a real matrix.
 *
 *	If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
 *	is diagonal and the eigenvector matrix V is orthogonal (i.e.
 *	A = V.times(D.times(V.transpose())) and V.times(V.transpose())
 *	equals the identity matrix).
 *
 *	If A is not symmetric, then the eigenvalue matrix D is block diagonal
 *	with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
 *	lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
 *	columns of V represent the eigenvectors in the sense that A*V = V*D,
 *	i.e. A.times(V) equals V.times(D).  The matrix V may be badly
 *	conditioned, or even singular, so the validity of the equation
 *	A = V*D*inverse(V) depends upon V.cond().
 *
 *	@author  Paul Meagher
 *	@license PHP v3.0
 *	@version 1.1
 *
 *  Slightly changed to adapt the original code to PHP-ML library
 *  @date 2017/04/11
 *  @author Mustafa Karabulut
 */

namespace Phpml\Math\LinearAlgebra;

use Phpml\Math\Matrix;

class EigenvalueDecomposition
{
    /**
     *	Row and column dimension (square matrix).
     *	@var int
     */
    private $n;

    /**
     *	Internal symmetry flag.
     *	@var bool
     */
    private $symmetric;

    /**
     *	Arrays for internal storage of eigenvalues.
     *	@var array
     */
    private $d = [];
    private $e = [];

    /**
     *	Array for internal storage of eigenvectors.
     *	@var array
     */
    private $V = [];

    /**
    *	Array for internal storage of nonsymmetric Hessenberg form.
    *	@var array
    */
    private $H = [];

    /**
    *	Working storage for nonsymmetric algorithm.
    *	@var array
    */
    private $ort;

    /**
    *	Used for complex scalar division.
    *	@var float
    */
    private $cdivr;
    private $cdivi;

    /**
     *	Constructor: Check for symmetry, then construct the eigenvalue decomposition
     *
     * @param array $Arg
     */
    public function __construct(array $Arg)
    {
        $this->A = $Arg;

The property A does not exist. Did you maybe forget to declare it?

        $this->n = count($Arg[0]);
        $this->symmetric = true;

        for ($j = 0; ($j < $this->n) && $this->symmetric; ++$j) {
            for ($i = 0; ($i < $this->n) & $this->symmetric; ++$i) {
                $this->symmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
            }
        }

        if ($this->symmetric) {
            $this->V = $this->A;
            // Tridiagonalize.
            $this->tred2();
            // Diagonalize.
            $this->tql2();
        } else {
            $this->H = $this->A;
            $this->ort = [];
            // Reduce to Hessenberg form.
            $this->orthes();
            // Reduce Hessenberg to real Schur form.
            $this->hqr2();
        }
    }

    /**
     *	Symmetric Householder reduction to tridiagonal form.
     */
    private function tred2()
    {
        //  This is derived from the Algol procedures tred2 by
        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
        //  Fortran subroutine in EISPACK.
        $this->d = $this->V[$this->n - 1];
        // Householder reduction to tridiagonal form.
        for ($i = $this->n - 1; $i > 0; --$i) {
            $i_ = $i - 1;
            // Scale to avoid under/overflow.
            $h = $scale = 0.0;
            $scale += array_sum(array_map('abs', $this->d));
            if ($scale == 0.0) {
                $this->e[$i] = $this->d[$i_];
                $this->d = array_slice($this->V[$i_], 0, $i_);
                for ($j = 0; $j < $i; ++$j) {
                    $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
                }
            } else {
                // Generate Householder vector.
                for ($k = 0; $k < $i; ++$k) {
                    $this->d[$k] /= $scale;
                    $h += pow($this->d[$k], 2);
                }

                $f = $this->d[$i_];
                $g = sqrt($h);
                if ($f > 0) {
                    $g = -$g;
                }

                $this->e[$i] = $scale * $g;
                $h = $h - $f * $g;
                $this->d[$i_] = $f - $g;

                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] = 0.0;
                }
                // Apply similarity transformation to remaining columns.
                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $this->V[$j][$i] = $f;
                    $g = $this->e[$j] + $this->V[$j][$j] * $f;

                    for ($k = $j + 1; $k <= $i_; ++$k) {
                        $g += $this->V[$k][$j] * $this->d[$k];
                        $this->e[$k] += $this->V[$k][$j] * $f;
                    }
                    $this->e[$j] = $g;
                }

                $f = 0.0;
                if ($h === 0 || $h < 1e-32) {
                    $h = 1e-32;
                }

                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] /= $h;
                    $f += $this->e[$j] * $this->d[$j];
                }

                $hh = $f / (2 * $h);
                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] -= $hh * $this->d[$j];
                }
                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $g = $this->e[$j];
                    for ($k = $j; $k <= $i_; ++$k) {
                        $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
                    }
                    $this->d[$j] = $this->V[$i - 1][$j];
                    $this->V[$i][$j] = 0.0;
                }
            }
            $this->d[$i] = $h;
        }

        // Accumulate transformations.
        for ($i = 0; $i < $this->n - 1; ++$i) {
            $this->V[$this->n - 1][$i] = $this->V[$i][$i];
            $this->V[$i][$i] = 1.0;
            $h = $this->d[$i + 1];
            if ($h != 0.0) {
                for ($k = 0; $k <= $i; ++$k) {

This code seems to be duplicated across your project.

                    $this->d[$k] = $this->V[$k][$i + 1] / $h;
                }
                for ($j = 0; $j <= $i; ++$j) {
                    $g = 0.0;
                    for ($k = 0; $k <= $i; ++$k) {

This code seems to be duplicated across your project.

                        $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
                    }
                    for ($k = 0; $k <= $i; ++$k) {

This code seems to be duplicated across your project.

                        $this->V[$k][$j] -= $g * $this->d[$k];
                    }
                }
            }
            for ($k = 0; $k <= $i; ++$k) {
                $this->V[$k][$i + 1] = 0.0;
            }
        }

        $this->d = $this->V[$this->n - 1];
        $this->V[$this->n - 1] = array_fill(0, $j, 0.0);

The variable $j does not seem to be defined for all execution paths leading up to this point.

        $this->V[$this->n - 1][$this->n - 1] = 1.0;
        $this->e[0] = 0.0;
    }


    /**
     *	Symmetric tridiagonal QL algorithm.
     *
     *	This is derived from the Algol procedures tql2, by
     *	Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
     *	Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
     *	Fortran subroutine in EISPACK.
     */
    private function tql2()
    {
        for ($i = 1; $i < $this->n; ++$i) {
            $this->e[$i - 1] = $this->e[$i];
        }
        $this->e[$this->n - 1] = 0.0;
        $f = 0.0;
        $tst1 = 0.0;
        $eps  = pow(2.0, -52.0);

        for ($l = 0; $l < $this->n; ++$l) {
            // Find small subdiagonal element
            $tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
            $m = $l;
            while ($m < $this->n) {
                if (abs($this->e[$m]) <= $eps * $tst1) {
                    break;
                }
                ++$m;
            }
            // If m == l, $this->d[l] is an eigenvalue,
            // otherwise, iterate.
            if ($m > $l) {
                $iter = 0;
                do {
                    // Could check iteration count here.
                    $iter += 1;
                    // Compute implicit shift
                    $g = $this->d[$l];
                    $p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
                    $r = hypot($p, 1.0);
                    if ($p < 0) {
                        $r *= -1;
                    }
                    $this->d[$l] = $this->e[$l] / ($p + $r);
                    $this->d[$l + 1] = $this->e[$l] * ($p + $r);
                    $dl1 = $this->d[$l + 1];
                    $h = $g - $this->d[$l];
                    for ($i = $l + 2; $i < $this->n; ++$i) {
                        $this->d[$i] -= $h;
                    }
                    $f += $h;
                    // Implicit QL transformation.
                    $p = $this->d[$m];
                    $c = 1.0;
                    $c2 = $c3 = $c;
                    $el1 = $this->e[$l + 1];
                    $s = $s2 = 0.0;
                    for ($i = $m - 1; $i >= $l; --$i) {
                        $c3 = $c2;
                        $c2 = $c;
                        $s2 = $s;
                        $g  = $c * $this->e[$i];
                        $h  = $c * $p;
                        $r  = hypot($p, $this->e[$i]);
                        $this->e[$i + 1] = $s * $r;
                        $s = $this->e[$i] / $r;
                        $c = $p / $r;
                        $p = $c * $this->d[$i] - $s * $g;
                        $this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
                        // Accumulate transformation.
                        for ($k = 0; $k < $this->n; ++$k) {
                            $h = $this->V[$k][$i + 1];
                            $this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
                            $this->V[$k][$i]     = $c * $this->V[$k][$i] - $s * $h;
                        }
                    }
                    $p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
                    $this->e[$l] = $s * $p;
                    $this->d[$l] = $c * $p;
                    // Check for convergence.
                } while (abs($this->e[$l]) > $eps * $tst1);
            }
            $this->d[$l] = $this->d[$l] + $f;
            $this->e[$l] = 0.0;
        }

        // Sort eigenvalues and corresponding vectors.
        for ($i = 0; $i < $this->n - 1; ++$i) {
            $k = $i;
            $p = $this->d[$i];
            for ($j = $i + 1; $j < $this->n; ++$j) {
                if ($this->d[$j] < $p) {
                    $k = $j;
                    $p = $this->d[$j];
                }
            }
            if ($k != $i) {
                $this->d[$k] = $this->d[$i];
                $this->d[$i] = $p;
                for ($j = 0; $j < $this->n; ++$j) {

This code seems to be duplicated across your project.

                    $p = $this->V[$j][$i];
                    $this->V[$j][$i] = $this->V[$j][$k];
                    $this->V[$j][$k] = $p;
                }
            }
        }
    }


    /**
     *	Nonsymmetric reduction to Hessenberg form.
     *
     *	This is derived from the Algol procedures orthes and ortran,
     *	by Martin and Wilkinson, Handbook for Auto. Comp.,
     *	Vol.ii-Linear Algebra, and the corresponding
     *	Fortran subroutines in EISPACK.
     */
    private function orthes()
    {
        $low  = 0;
        $high = $this->n - 1;

        for ($m = $low + 1; $m <= $high - 1; ++$m) {
            // Scale column.
            $scale = 0.0;
            for ($i = $m; $i <= $high; ++$i) {

This code seems to be duplicated across your project.

                $scale = $scale + abs($this->H[$i][$m - 1]);
            }
            if ($scale != 0.0) {
                // Compute Householder transformation.
                $h = 0.0;
                for ($i = $high; $i >= $m; --$i) {
                    $this->ort[$i] = $this->H[$i][$m - 1] / $scale;
                    $h += $this->ort[$i] * $this->ort[$i];
                }
                $g = sqrt($h);
                if ($this->ort[$m] > 0) {
                    $g *= -1;
                }
                $h -= $this->ort[$m] * $g;
                $this->ort[$m] -= $g;
                // Apply Householder similarity transformation
                // H = (I -u * u' / h) * H * (I -u * u') / h)
                for ($j = $m; $j < $this->n; ++$j) {

This code seems to be duplicated across your project.

                    $f = 0.0;
                    for ($i = $high; $i >= $m; --$i) {
                        $f += $this->ort[$i] * $this->H[$i][$j];
                    }
                    $f /= $h;
                    for ($i = $m; $i <= $high; ++$i) {
                        $this->H[$i][$j] -= $f * $this->ort[$i];
                    }
                }
                for ($i = 0; $i <= $high; ++$i) {

This code seems to be duplicated across your project.

                    $f = 0.0;
                    for ($j = $high; $j >= $m; --$j) {
                        $f += $this->ort[$j] * $this->H[$i][$j];
                    }
                    $f = $f / $h;
                    for ($j = $m; $j <= $high; ++$j) {
                        $this->H[$i][$j] -= $f * $this->ort[$j];
                    }
                }
                $this->ort[$m] = $scale * $this->ort[$m];
                $this->H[$m][$m - 1] = $scale * $g;
            }
        }

        // Accumulate transformations (Algol's ortran).
        for ($i = 0; $i < $this->n; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                $this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
            }
        }
        for ($m = $high - 1; $m >= $low + 1; --$m) {
            if ($this->H[$m][$m - 1] != 0.0) {
                for ($i = $m + 1; $i <= $high; ++$i) {

This code seems to be duplicated across your project.

                    $this->ort[$i] = $this->H[$i][$m - 1];
                }
                for ($j = $m; $j <= $high; ++$j) {
                    $g = 0.0;
                    for ($i = $m; $i <= $high; ++$i) {
                        $g += $this->ort[$i] * $this->V[$i][$j];
                    }
                    // Double division avoids possible underflow
                    $g = ($g / $this->ort[$m]) / $this->H[$m][$m - 1];
                    for ($i = $m; $i <= $high; ++$i) {
                        $this->V[$i][$j] += $g * $this->ort[$i];
                    }
                }
            }
        }
    }

    /**
     * Performs complex division.
     *
     * @param int|float $xr
     * @param int|float $xi
     * @param int|float $yr
     * @param int|float $yi
     */
    private function cdiv($xr, $xi, $yr, $yi)
    {
        if (abs($yr) > abs($yi)) {
            $r = $yi / $yr;
            $d = $yr + $r * $yi;
            $this->cdivr = ($xr + $r * $xi) / $d;
            $this->cdivi = ($xi - $r * $xr) / $d;
        } else {
            $r = $yr / $yi;
            $d = $yi + $r * $yr;
            $this->cdivr = ($r * $xr + $xi) / $d;
            $this->cdivi = ($r * $xi - $xr) / $d;
        }
    }

    /**
     *	Nonsymmetric reduction from Hessenberg to real Schur form.
     *
     *	Code is derived from the Algol procedure hqr2,
     *	by Martin and Wilkinson, Handbook for Auto. Comp.,
     *	Vol.ii-Linear Algebra, and the corresponding
     *	Fortran subroutine in EISPACK.
     */
    private function hqr2()
    {
        //  Initialize
        $nn = $this->n;
        $n  = $nn - 1;
        $low = 0;
        $high = $nn - 1;
        $eps = pow(2.0, -52.0);
        $exshift = 0.0;
        $p = $q = $r = $s = $z = 0;
        // Store roots isolated by balanc and compute matrix norm
        $norm = 0.0;

        for ($i = 0; $i < $nn; ++$i) {
            if (($i < $low) or ($i > $high)) {

Using logical operators such as or instead of || is generally not recommended.

                $this->d[$i] = $this->H[$i][$i];
                $this->e[$i] = 0.0;
            }
            for ($j = max($i - 1, 0); $j < $nn; ++$j) {
                $norm = $norm + abs($this->H[$i][$j]);
            }
        }

        // Outer loop over eigenvalue index
        $iter = 0;
        while ($n >= $low) {
            // Look for single small sub-diagonal element
            $l = $n;
            while ($l > $low) {
                $s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
                if ($s == 0.0) {
                    $s = $norm;
                }
                if (abs($this->H[$l][$l - 1]) < $eps * $s) {
                    break;
                }
                --$l;
            }
            // Check for convergence
            // One root found
            if ($l == $n) {
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->d[$n] = $this->H[$n][$n];
                $this->e[$n] = 0.0;
                --$n;
                $iter = 0;
                // Two roots found
            } elseif ($l == $n - 1) {
                $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
                $p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
                $q = $p * $p + $w;
                $z = sqrt(abs($q));
                $this->H[$n][$n] = $this->H[$n][$n] + $exshift;
                $this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift;
                $x = $this->H[$n][$n];
                // Real pair
                if ($q >= 0) {
                    if ($p >= 0) {
                        $z = $p + $z;
                    } else {
                        $z = $p - $z;
                    }
                    $this->d[$n - 1] = $x + $z;
                    $this->d[$n] = $this->d[$n - 1];
                    if ($z != 0.0) {
                        $this->d[$n] = $x - $w / $z;
                    }
                    $this->e[$n - 1] = 0.0;
                    $this->e[$n] = 0.0;
                    $x = $this->H[$n][$n - 1];
                    $s = abs($x) + abs($z);
                    $p = $x / $s;
                    $q = $z / $s;
                    $r = sqrt($p * $p + $q * $q);
                    $p = $p / $r;
                    $q = $q / $r;
                    // Row modification
                    for ($j = $n - 1; $j < $nn; ++$j) {

This code seems to be duplicated across your project.

                        $z = $this->H[$n - 1][$j];
                        $this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
                        $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
                    }
                    // Column modification
                    for ($i = 0; $i <= $n; ++$i) {

This code seems to be duplicated across your project.

                        $z = $this->H[$i][$n - 1];
                        $this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
                        $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
                    }
                    // Accumulate transformations
                    for ($i = $low; $i <= $high; ++$i) {

This code seems to be duplicated across your project.

                        $z = $this->V[$i][$n - 1];
                        $this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
                        $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
                    }
                    // Complex pair
                } else {
                    $this->d[$n - 1] = $x + $p;
                    $this->d[$n]     = $x + $p;
                    $this->e[$n - 1] = $z;
                    $this->e[$n]     = -$z;
                }
                $n = $n - 2;
                $iter = 0;
                // No convergence yet
            } else {
                // Form shift
                $x = $this->H[$n][$n];
                $y = 0.0;
                $w = 0.0;
                if ($l < $n) {
                    $y = $this->H[$n - 1][$n - 1];
                    $w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
                }
                // Wilkinson's original ad hoc shift
                if ($iter == 10) {
                    $exshift += $x;
                    for ($i = $low; $i <= $n; ++$i) {
                        $this->H[$i][$i] -= $x;
                    }
                    $s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
                    $x = $y = 0.75 * $s;
                    $w = -0.4375 * $s * $s;
                }
                // MATLAB's new ad hoc shift
                if ($iter == 30) {
                    $s = ($y - $x) / 2.0;
                    $s = $s * $s + $w;
                    if ($s > 0) {
                        $s = sqrt($s);
                        if ($y < $x) {
                            $s = -$s;
                        }
                        $s = $x - $w / (($y - $x) / 2.0 + $s);
                        for ($i = $low; $i <= $n; ++$i) {
                            $this->H[$i][$i] -= $s;
                        }
                        $exshift += $s;
                        $x = $y = $w = 0.964;
                    }
                }
                // Could check iteration count here.
                $iter = $iter + 1;
                // Look for two consecutive small sub-diagonal elements
                $m = $n - 2;
                while ($m >= $l) {
                    $z = $this->H[$m][$m];
                    $r = $x - $z;
                    $s = $y - $z;
                    $p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
                    $q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
                    $r = $this->H[$m + 2][$m + 1];
                    $s = abs($p) + abs($q) + abs($r);
                    $p = $p / $s;
                    $q = $q / $s;
                    $r = $r / $s;
                    if ($m == $l) {
                        break;
                    }
                    if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) <
                        $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1])))) {
                        break;
                    }
                    --$m;
                }
                for ($i = $m + 2; $i <= $n; ++$i) {
                    $this->H[$i][$i - 2] = 0.0;
                    if ($i > $m + 2) {
                        $this->H[$i][$i - 3] = 0.0;
                    }
                }
                // Double QR step involving rows l:n and columns m:n
                for ($k = $m; $k <= $n - 1; ++$k) {
                    $notlast = ($k != $n - 1);
                    if ($k != $m) {
                        $p = $this->H[$k][$k - 1];
                        $q = $this->H[$k + 1][$k - 1];
                        $r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
                        $x = abs($p) + abs($q) + abs($r);
                        if ($x != 0.0) {
                            $p = $p / $x;
                            $q = $q / $x;
                            $r = $r / $x;
                        }
                    }
                    if ($x == 0.0) {
                        break;
                    }
                    $s = sqrt($p * $p + $q * $q + $r * $r);
                    if ($p < 0) {
                        $s = -$s;
                    }
                    if ($s != 0) {
                        if ($k != $m) {
                            $this->H[$k][$k - 1] = -$s * $x;
                        } elseif ($l != $m) {
                            $this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
                        }
                        $p = $p + $s;
                        $x = $p / $s;
                        $y = $q / $s;
                        $z = $r / $s;
                        $q = $q / $p;
                        $r = $r / $p;
                        // Row modification
                        for ($j = $k; $j < $nn; ++$j) {

This code seems to be duplicated across your project.

                            $p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
                            if ($notlast) {
                                $p = $p + $r * $this->H[$k + 2][$j];
                                $this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
                            }
                            $this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
                            $this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
                        }
                        // Column modification
                        for ($i = 0; $i <= min($n, $k + 3); ++$i) {

This code seems to be duplicated across your project.

                            $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
                            if ($notlast) {
                                $p = $p + $z * $this->H[$i][$k + 2];
                                $this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r;
                            }
                            $this->H[$i][$k] = $this->H[$i][$k] - $p;
                            $this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q;
                        }
                        // Accumulate transformations
                        for ($i = $low; $i <= $high; ++$i) {

This code seems to be duplicated across your project.

                            $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
                            if ($notlast) {
                                $p = $p + $z * $this->V[$i][$k + 2];
                                $this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
                            }
                            $this->V[$i][$k] = $this->V[$i][$k] - $p;
                            $this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q;
                        }
                    }  // ($s != 0)

56% of this comment could be valid code. Did you maybe forget this after debugging?

                }  // k loop
            }  // check convergence
        }  // while ($n >= $low)

46% of this comment could be valid code. Did you maybe forget this after debugging?

        // Backsubstitute to find vectors of upper triangular form
        if ($norm == 0.0) {
            return;
        }

        for ($n = $nn - 1; $n >= 0; --$n) {
            $p = $this->d[$n];
            $q = $this->e[$n];
            // Real vector
            if ($q == 0) {
                $l = $n;
                $this->H[$n][$n] = 1.0;
                for ($i = $n - 1; $i >= 0; --$i) {
                    $w = $this->H[$i][$i] - $p;
                    $r = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $r = $r + $this->H[$i][$j] * $this->H[$j][$n];
                    }

                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $s = $r;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0.0) {
                            if ($w != 0.0) {
                                $this->H[$i][$n] = -$r / $w;
                            } else {
                                $this->H[$i][$n] = -$r / ($eps * $norm);
                            }
                            // Solve real equations
                        } else {
                            $x = $this->H[$i][$i + 1];
                            $y = $this->H[$i + 1][$i];
                            $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
                            $t = ($x * $s - $z * $r) / $q;
                            $this->H[$i][$n] = $t;
                            if (abs($x) > abs($z)) {
                                $this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
                            } else {
                                $this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
                            }
                        }
                        // Overflow control
                        $t = abs($this->H[$i][$n]);
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n] = $this->H[$j][$n] / $t;
                            }
                        }
                    }
                }
                // Complex vector
            } elseif ($q < 0) {
                $l = $n - 1;
                // Last vector component imaginary so matrix is triangular
                if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n])) {
                    $this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
                    $this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
                } else {
                    $this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
                    $this->H[$n - 1][$n - 1] = $this->cdivr;
                    $this->H[$n - 1][$n]     = $this->cdivi;
                }
                $this->H[$n][$n - 1] = 0.0;
                $this->H[$n][$n]     = 1.0;
                for ($i = $n - 2; $i >= 0; --$i) {
                    // double ra,sa,vr,vi;

37% of this comment could be valid code. Did you maybe forget this after debugging?

                    $ra = 0.0;
                    $sa = 0.0;
                    for ($j = $l; $j <= $n; ++$j) {
                        $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n - 1];
                        $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
                    }
                    $w = $this->H[$i][$i] - $p;
                    if ($this->e[$i] < 0.0) {
                        $z = $w;
                        $r = $ra;
                        $s = $sa;
                    } else {
                        $l = $i;
                        if ($this->e[$i] == 0) {
                            $this->cdiv(-$ra, -$sa, $w, $q);
                            $this->H[$i][$n - 1] = $this->cdivr;
                            $this->H[$i][$n]     = $this->cdivi;
                        } else {
                            // Solve complex equations
                            $x = $this->H[$i][$i + 1];
                            $y = $this->H[$i + 1][$i];
                            $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
                            $vi = ($this->d[$i] - $p) * 2.0 * $q;
                            if ($vr == 0.0 & $vi == 0.0) {

Consider adding parentheses for clarity. Current Interpretation: ($vr == 0.0) & $vi == 0.0, Probably Intended Meaning: $vr == (0.0 & $vi == 0.0)

                                $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
                            }
                            $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
                            $this->H[$i][$n - 1] = $this->cdivr;
                            $this->H[$i][$n]     = $this->cdivi;
                            if (abs($x) > (abs($z) + abs($q))) {
                                $this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
                                $this->H[$i + 1][$n]     = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
                            } else {
                                $this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
                                $this->H[$i + 1][$n - 1] = $this->cdivr;
                                $this->H[$i + 1][$n]     = $this->cdivi;
                            }
                        }
                        // Overflow control
                        $t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
                        if (($eps * $t) * $t > 1) {
                            for ($j = $i; $j <= $n; ++$j) {
                                $this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t;
                                $this->H[$j][$n]     = $this->H[$j][$n] / $t;
                            }
                        }
                    } // end else
                } // end for
            } // end else for complex case
        } // end for

        // Vectors of isolated roots
        for ($i = 0; $i < $nn; ++$i) {
            if ($i < $low | $i > $high) {
                for ($j = $i; $j < $nn; ++$j) {
                    $this->V[$i][$j] = $this->H[$i][$j];
                }
            }
        }

        // Back transformation to get eigenvectors of original matrix
        for ($j = $nn - 1; $j >= $low; --$j) {
            for ($i = $low; $i <= $high; ++$i) {
                $z = 0.0;
                for ($k = $low; $k <= min($j, $high); ++$k) {
                    $z = $z + $this->V[$i][$k] * $this->H[$k][$j];
                }
                $this->V[$i][$j] = $z;
            }
        }
    } // end hqr2

    /**
     * Return the eigenvector matrix
     *
     * @access public
     *
     * @return array
     */
    public function getEigenvectors()
    {
        $vectors = $this->V;

        // Always return the eigenvectors of length 1.0
        $vectors = new Matrix($vectors);
        $vectors = array_map(function ($vect) {
            $sum = 0;
            for ($i = 0; $i < count($vect); ++$i) {

It seems like you are calling the size function count() as part of the test condition. You might want to compute the size beforehand, and not on each iteration.

This code seems to be duplicated across your project.

                $sum += $vect[$i] ** 2;
            }

            $sum = sqrt($sum);
            for ($i = 0; $i < count($vect); ++$i) {

It seems like you are calling the size function count() as part of the test condition. You might want to compute the size beforehand, and not on each iteration.

This code seems to be duplicated across your project.

                $vect[$i] /= $sum;
            }

            return $vect;
        }, $vectors->transpose()->toArray());

        return $vectors;
    }

    /**
     *	Return the real parts of the eigenvalues<br>
     *  d = real(diag(D));
     *
     *	@return array
     */
    public function getRealEigenvalues()
    {
        return $this->d;
    }

    /**
     *	Return the imaginary parts of the eigenvalues <br>
     *  d = imag(diag(D))
     *
     *	@return array
     */
    public function getImagEigenvalues()
    {
        return $this->e;
    }

    /**
     *	Return the block diagonal eigenvalue matrix
     *
     *	@return array
     */
    public function getDiagonalEigenvalues()
    {
        $D = [];

        for ($i = 0; $i < $this->n; ++$i) {
            $D[$i] = array_fill(0, $this->n, 0.0);
            $D[$i][$i] = $this->d[$i];
            if ($this->e[$i] == 0) {
                continue;
            }

            $o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
            $D[$i][$o] = $this->e[$i];
        }

        return $D;
    }
}    //	class EigenvalueDecomposition