Return the diagonal matrix of singular values

@access public
@return S
347831
public function getS()
    {
        for ($i = 0; $i < $this->n; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                $S[$i][$j] = 0.0;
            }
            $S[$i][$i] = $this->s[$i];
        }
        return new Matrix($S);
    }
320578
public static int rank( DMatrixRMaj A  ) {
        SingularValueDecomposition_F64<DMatrixRMaj> svd = DecompositionFactory_DDRM.svd(A.numRows,A.numCols,false,true,true);

        if( svd.inputModified() ) {
            A = A.copy();
        }
        if( !svd.decompose(A)) {
            throw new RuntimeException("SVD Failed!");
        }

        int N = svd.numberOfSingularValues();
        double sv[] = svd.getSingularValues();

        double threshold = singularThreshold(sv,N);
        int count = 0;
        for (int i = 0; i < sv.length; i++) {
            if( sv[i] >= threshold ) {
                count++;
            }
        }
        return count;
    }
500618
function svdJs() {
                var A = this;
                var V = Matrix$b.I(A.rows());
                var S = A.transpose();
                var U = Matrix$b.I(A.cols());
                var err = Number.MAX_VALUE;
                var i = 0;
                var maxLoop = 100;

                while(err > 2.2737e-13 && i < maxLoop) {
                    var qr = S.transpose().qrJs();
                    S = qr.R;
                    V = V.x(qr.Q);
                    qr = S.transpose().qrJs();
                    U = U.x(qr.Q);
                    S = qr.R;

                    var e = S.triu(1).unroll().norm();
                    var f = S.diagonal().norm();

                    if(f == 0)
                        f = 1;

                    err = e / f;

                    i++;
                }

                var ss = S.diagonal();
                var s = [];

                for(var i = 1; i <= ss.cols(); i++) {
                    var ssn = ss.e(i);
                    s.push(Math.abs(ssn));

                    if(ssn < 0) {
                        for(var j = 0; j < U.rows(); j++) {
                            V.elements[j][i - 1] = -(V.elements[j][i - 1]);
                        }
                    }
                }

                return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
            }
1194756
public function rref() : self
    {
        [$b, $swaps] = $this->ref();

        $b = $b->asArray();

        $row = $col = 0;

        while ($row < $this->m and $col < $this->n) {
            $t = $b[$row];

            if (abs($t[$col]) == 0) {
                $col++;

                continue 1;
            }

            $divisor = $t[$col];

            if ($divisor != 1) {
                for ($i = 0; $i < $this->n; $i++) {
                    $t[$i] /= $divisor;
                }
            }

            for ($i = $row - 1; $i >= 0; $i--) {
                $scale = $b[$i][$col];

                if ($scale != 0) {
                    for ($j = 0; $j < $this->n; $j++) {
                        $b[$i][$j] += -$scale * $t[$j];
                    }
                }
            }

            $b[$row] = $t;

            $row++;
            $col++;
        }

        return self::quick($b);
    }
1092698
public function getL()
    {
        $L = [];
        for ($i = 0; $i < $this->m; ++$i) {
            for ($j = 0; $j < $this->n; ++$j) {
                if ($i > $j) {
                    $L[$i][$j] = $this->LU[$i][$j];
                } elseif ($i == $j) {
                    $L[$i][$j] = 1.0;
                } else {
                    $L[$i][$j] = 0.0;
                }
            }
        }
        return new Matrix($L);
    }
578026
private function getSamplesMatrix(): Matrix
    {
        $samples = [];
        foreach ($this->samples as $sample) {
            array_unshift($sample, 1);
            $samples[] = $sample;
        }

        return new Matrix($samples);
    }
34458
public function getDeterminant()
    {
        if ($this->determinant !== null) {
            return $this->determinant;
        }

        if (!$this->isSquare()) {
            throw new MatrixException('Matrix is not square matrix');
        }

        $lu = new LUDecomposition($this);

        return $this->determinant = $lu->det();
    }
1092191
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)) {
                $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) {
                        $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) {
                        $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) {
                        $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) {
                            $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) {
                            $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) {
                            $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)
                }  // k loop
            }  // check convergence
        }  // while ($n >= $low)

        // 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;
                    $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) {
                                $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;
            }
        }
    }
1219908
public function median() : ColumnVector
    {
        if ($this->n < 1) {
            throw new RuntimeException('Median is not defined for matrices'
                . ' with less than 1 column.');
        }

        $b = [];

        foreach ($this->a as $row) {
            $mid = intdiv($this->n, 2);

            sort($row);

            if ($this->n % 2 === 1) {
                $median = $row[$mid];
            } else {
                $median = ($row[$mid - 1] + $row[$mid]) / 2.;
            }

            $b[] = $median;
        }

        return ColumnVector::quick($b);
    }
1330754
public function det()
    {
        if ($this->m !== $this->n) {
            throw MatrixException::notSquareMatrix();
        }

        $d = $this->pivsign;
        for ($j = 0; $j < $this->n; ++$j) {
            $d *= $this->LU[$j][$j];
        }

        return $d;
    }
622682
private void makeSingularPositive() {
        numSingular = qralg.getNumberOfSingularValues();
        singularValues = qralg.getSingularValues();

        for( int i = 0; i < numSingular; i++ ) {
            double val = singularValues[i];

            if( val < 0 ) {
                singularValues[i] = -val;

                if( computeU ) {
                    // compute the results of multiplying it by an element of -1 at this location in
                    // a diagonal matrix.
                    int start = i* Ut.numCols;
                    int stop = start+ Ut.numCols;

                    for( int j = start; j < stop; j++ ) {
                        Ut.data[j] = -Ut.data[j];
                    }
                }
            }
        }
    }
929254
public static FloatMatrix SVDValues(FloatMatrix A) {
        int m = A.rows;
        int n = A.columns;
        FloatMatrix S = new FloatMatrix(min(m, n));

        int info = NativeBlas.sgesvd('N', 'N', m, n, A.dup().data, 0, m, S.data, 0, null, 0, 1, null, 0, 1);

        if (info > 0) {
          throw new LapackConvergenceException("GESVD", info + " superdiagonals of an intermediate bidiagonal form failed to converge.");
        }

        return S;
    }
226225
public function getU()
    {
        return new Matrix($this->U, $this->m, min($this->m + 1, $this->n));
    }
241424
public function getDiagonalEigenvalues(): array
    {
        $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;
    }
1325660
public function getBruttoSum()
    {
        $dSum = 0;
        foreach ($this->_aList as $oPrice) {
            $dSum += $oPrice->getBruttoPrice();
        }

        return $dSum;
    }
277697
public Matrix getIminus() {
        Matrix Iminus = new Matrix(size, size);
        int i, j;
        for (i = 0; i < size; i++)
            for (j = 0; j < size; j++)
                Iminus.matrix[i][j] = 0d;
        for (i = 1; i < size; i++)
            Iminus.matrix[i][i - 1] = Math.sqrt(J * J + J - (J - i) * (J - i) - (J - i));
        return Iminus;
    }
1185311
@Nonnull
  @ReturnsMutableCopy
  public Matrix getL ()
  {
    final Matrix aNewMatrix = new Matrix (m_nRows, m_nCols);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();
    for (int nRow = 0; nRow < m_nRows; nRow++)
    {
      final double [] aSrcRow = m_aLU[nRow];
      final double [] aDstRow = aNewArray[nRow];
      for (int nCol = 0; nCol < m_nCols; nCol++)
        if (nRow > nCol)
          aDstRow[nCol] = aSrcRow[nCol];
        else
          if (nRow == nCol)
            aDstRow[nCol] = 1.0;
          else
            aDstRow[nCol] = 0.0;
    }
    return aNewMatrix;
  }