Householder.php 1.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657
  1. <?php
  2. namespace MathPHP\LinearAlgebra;
  3. use MathPHP\Exception;
  4. /**
  5. * A Householder transformation (also known as a Householder reflection or elementary reflector) is a linear
  6. * transformation that describes a reflection about a plane or hyperplane containing the origin.
  7. */
  8. class Householder
  9. {
  10. /**
  11. * Householder matrix transformation
  12. *
  13. * u = x ± αe where α = ‖x‖ and sgn(α) = sgn(x)
  14. *
  15. * uuᵀ
  16. * Q = I - 2 * -----
  17. * uᵀu
  18. *
  19. * https://en.wikipedia.org/wiki/Householder_transformation
  20. *
  21. * @param NumericMatrix $A source Matrix
  22. *
  23. * @return NumericMatrix
  24. *
  25. * @throws Exception\MathException
  26. */
  27. public static function transform(NumericMatrix $A): NumericMatrix
  28. {
  29. $m = $A->getM();
  30. $I = MatrixFactory::identity($m);
  31. // x is the leftmost column of A
  32. $x = $A->submatrix(0, 0, $m - 1, 0);
  33. // α is the square root of the sum of squares of x with the correct sign
  34. $sign = $x[0][0] >= 0 ? 1 : -1;
  35. $α = $sign * $x->frobeniusNorm();
  36. // e is the first column of I
  37. $e = $I->submatrix(0, 0, $m - 1, 0);
  38. // u = x ± αe
  39. $u = $e->scalarMultiply($α)->add($x);
  40. $uᵀ = $u->transpose();
  41. $uᵀu = $uᵀ->multiply($u)->get(0, 0);
  42. $uuᵀ = $u->multiply($uᵀ);
  43. if ($uᵀu == 0) {
  44. return $I;
  45. }
  46. // We scale $uuᵀ and subtract it from the identity matrix
  47. return $I->subtract($uuᵀ->scalarMultiply(2 / $uᵀu));
  48. }
  49. }