sqrtm.js 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899
  1. import { isMatrix } from '../../utils/is.js';
  2. import { format } from '../../utils/string.js';
  3. import { arraySize } from '../../utils/array.js';
  4. import { factory } from '../../utils/factory.js';
  5. var name = 'sqrtm';
  6. var dependencies = ['typed', 'abs', 'add', 'multiply', 'map', 'sqrt', 'subtract', 'inv', 'size', 'max', 'identity'];
  7. export var createSqrtm = /* #__PURE__ */factory(name, dependencies, _ref => {
  8. var {
  9. typed,
  10. abs,
  11. add,
  12. multiply,
  13. map,
  14. sqrt,
  15. subtract,
  16. inv,
  17. size,
  18. max,
  19. identity
  20. } = _ref;
  21. var _maxIterations = 1e3;
  22. var _tolerance = 1e-6;
  23. /**
  24. * Calculate the principal square root matrix using the Denman–Beavers iterative method
  25. *
  26. * https://en.wikipedia.org/wiki/Square_root_of_a_matrix#By_Denman–Beavers_iteration
  27. *
  28. * @param {Array | Matrix} A The square matrix `A`
  29. * @return {Array | Matrix} The principal square root of matrix `A`
  30. * @private
  31. */
  32. function _denmanBeavers(A) {
  33. var error;
  34. var iterations = 0;
  35. var Y = A;
  36. var Z = identity(size(A));
  37. do {
  38. var Yk = Y;
  39. Y = multiply(0.5, add(Yk, inv(Z)));
  40. Z = multiply(0.5, add(Z, inv(Yk)));
  41. error = max(abs(subtract(Y, Yk)));
  42. if (error > _tolerance && ++iterations > _maxIterations) {
  43. throw new Error('computing square root of matrix: iterative method could not converge');
  44. }
  45. } while (error > _tolerance);
  46. return Y;
  47. }
  48. /**
  49. * Calculate the principal square root of a square matrix.
  50. * The principal square root matrix `X` of another matrix `A` is such that `X * X = A`.
  51. *
  52. * https://en.wikipedia.org/wiki/Square_root_of_a_matrix
  53. *
  54. * Syntax:
  55. *
  56. * X = math.sqrtm(A)
  57. *
  58. * Examples:
  59. *
  60. * math.sqrtm([[33, 24], [48, 57]]) // returns [[5, 2], [4, 7]]
  61. *
  62. * See also:
  63. *
  64. * sqrt, pow
  65. *
  66. * @param {Array | Matrix} A The square matrix `A`
  67. * @return {Array | Matrix} The principal square root of matrix `A`
  68. */
  69. return typed(name, {
  70. 'Array | Matrix': function ArrayMatrix(A) {
  71. var size = isMatrix(A) ? A.size() : arraySize(A);
  72. switch (size.length) {
  73. case 1:
  74. // Single element Array | Matrix
  75. if (size[0] === 1) {
  76. return map(A, sqrt);
  77. } else {
  78. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  79. }
  80. case 2:
  81. {
  82. // Two-dimensional Array | Matrix
  83. var rows = size[0];
  84. var cols = size[1];
  85. if (rows === cols) {
  86. return _denmanBeavers(A);
  87. } else {
  88. throw new RangeError('Matrix must be square ' + '(size: ' + format(size) + ')');
  89. }
  90. }
  91. default:
  92. // Multi dimensional array
  93. throw new RangeError('Matrix must be at most two dimensional ' + '(size: ' + format(size) + ')');
  94. }
  95. }
  96. });
  97. });