eigs.js 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256
  1. import { factory } from '../../utils/factory.js';
  2. import { format } from '../../utils/string.js';
  3. import { createComplexEigs } from './eigs/complexEigs.js';
  4. import { createRealSymmetric } from './eigs/realSymetric.js';
  5. import { typeOf, isNumber, isBigNumber, isComplex, isFraction } from '../../utils/is.js';
  6. var name = 'eigs';
  7. // The absolute state of math.js's dependency system:
  8. var dependencies = ['config', 'typed', 'matrix', 'addScalar', 'equal', 'subtract', 'abs', 'atan', 'cos', 'sin', 'multiplyScalar', 'divideScalar', 'inv', 'bignumber', 'multiply', 'add', 'larger', 'column', 'flatten', 'number', 'complex', 'sqrt', 'diag', 'qr', 'usolve', 'usolveAll', 'im', 're', 'smaller', 'matrixFromColumns', 'dot'];
  9. export var createEigs = /* #__PURE__ */factory(name, dependencies, _ref => {
  10. var {
  11. config,
  12. typed,
  13. matrix,
  14. addScalar,
  15. subtract,
  16. equal,
  17. abs,
  18. atan,
  19. cos,
  20. sin,
  21. multiplyScalar,
  22. divideScalar,
  23. inv,
  24. bignumber,
  25. multiply,
  26. add,
  27. larger,
  28. column,
  29. flatten,
  30. number,
  31. complex,
  32. sqrt,
  33. diag,
  34. qr,
  35. usolve,
  36. usolveAll,
  37. im,
  38. re,
  39. smaller,
  40. matrixFromColumns,
  41. dot
  42. } = _ref;
  43. var doRealSymetric = createRealSymmetric({
  44. config,
  45. addScalar,
  46. subtract,
  47. column,
  48. flatten,
  49. equal,
  50. abs,
  51. atan,
  52. cos,
  53. sin,
  54. multiplyScalar,
  55. inv,
  56. bignumber,
  57. complex,
  58. multiply,
  59. add
  60. });
  61. var doComplexEigs = createComplexEigs({
  62. config,
  63. addScalar,
  64. subtract,
  65. multiply,
  66. multiplyScalar,
  67. flatten,
  68. divideScalar,
  69. sqrt,
  70. abs,
  71. bignumber,
  72. diag,
  73. qr,
  74. inv,
  75. usolve,
  76. usolveAll,
  77. equal,
  78. complex,
  79. larger,
  80. smaller,
  81. matrixFromColumns,
  82. dot
  83. });
  84. /**
  85. * Compute eigenvalues and eigenvectors of a matrix. The eigenvalues are sorted by their absolute value, ascending.
  86. * An eigenvalue with multiplicity k will be listed k times. The eigenvectors are returned as columns of a matrix –
  87. * the eigenvector that belongs to the j-th eigenvalue in the list (eg. `values[j]`) is the j-th column (eg. `column(vectors, j)`).
  88. * If the algorithm fails to converge, it will throw an error – in that case, however, you may still find useful information
  89. * in `err.values` and `err.vectors`.
  90. *
  91. * Syntax:
  92. *
  93. * math.eigs(x, [prec])
  94. *
  95. * Examples:
  96. *
  97. * const { eigs, multiply, column, transpose } = math
  98. * const H = [[5, 2.3], [2.3, 1]]
  99. * const ans = eigs(H) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}
  100. * const E = ans.values
  101. * const U = ans.vectors
  102. * multiply(H, column(U, 0)) // returns multiply(E[0], column(U, 0))
  103. * const UTxHxU = multiply(transpose(U), H, U) // diagonalizes H
  104. * E[0] == UTxHxU[0][0] // returns true
  105. *
  106. * See also:
  107. *
  108. * inv
  109. *
  110. * @param {Array | Matrix} x Matrix to be diagonalized
  111. *
  112. * @param {number | BigNumber} [prec] Precision, default value: 1e-15
  113. * @return {{values: Array|Matrix, vectors: Array|Matrix}} Object containing an array of eigenvalues and a matrix with eigenvectors as columns.
  114. *
  115. */
  116. return typed('eigs', {
  117. Array: function Array(x) {
  118. var mat = matrix(x);
  119. return computeValuesAndVectors(mat);
  120. },
  121. 'Array, number|BigNumber': function ArrayNumberBigNumber(x, prec) {
  122. var mat = matrix(x);
  123. return computeValuesAndVectors(mat, prec);
  124. },
  125. Matrix: function Matrix(mat) {
  126. var {
  127. values,
  128. vectors
  129. } = computeValuesAndVectors(mat);
  130. return {
  131. values: matrix(values),
  132. vectors: matrix(vectors)
  133. };
  134. },
  135. 'Matrix, number|BigNumber': function MatrixNumberBigNumber(mat, prec) {
  136. var {
  137. values,
  138. vectors
  139. } = computeValuesAndVectors(mat, prec);
  140. return {
  141. values: matrix(values),
  142. vectors: matrix(vectors)
  143. };
  144. }
  145. });
  146. function computeValuesAndVectors(mat, prec) {
  147. if (prec === undefined) {
  148. prec = config.epsilon;
  149. }
  150. var size = mat.size();
  151. if (size.length !== 2 || size[0] !== size[1]) {
  152. throw new RangeError('Matrix must be square (size: ' + format(size) + ')');
  153. }
  154. var arr = mat.toArray();
  155. var N = size[0];
  156. if (isReal(arr, N, prec)) {
  157. coerceReal(arr, N);
  158. if (isSymmetric(arr, N, prec)) {
  159. var _type = coerceTypes(mat, arr, N);
  160. return doRealSymetric(arr, N, prec, _type);
  161. }
  162. }
  163. var type = coerceTypes(mat, arr, N);
  164. return doComplexEigs(arr, N, prec, type);
  165. }
  166. /** @return {boolean} */
  167. function isSymmetric(arr, N, prec) {
  168. for (var i = 0; i < N; i++) {
  169. for (var j = i; j < N; j++) {
  170. // TODO proper comparison of bignum and frac
  171. if (larger(bignumber(abs(subtract(arr[i][j], arr[j][i]))), prec)) {
  172. return false;
  173. }
  174. }
  175. }
  176. return true;
  177. }
  178. /** @return {boolean} */
  179. function isReal(arr, N, prec) {
  180. for (var i = 0; i < N; i++) {
  181. for (var j = 0; j < N; j++) {
  182. // TODO proper comparison of bignum and frac
  183. if (larger(bignumber(abs(im(arr[i][j]))), prec)) {
  184. return false;
  185. }
  186. }
  187. }
  188. return true;
  189. }
  190. function coerceReal(arr, N) {
  191. for (var i = 0; i < N; i++) {
  192. for (var j = 0; j < N; j++) {
  193. arr[i][j] = re(arr[i][j]);
  194. }
  195. }
  196. }
  197. /** @return {'number' | 'BigNumber' | 'Complex'} */
  198. function coerceTypes(mat, arr, N) {
  199. /** @type {string} */
  200. var type = mat.datatype();
  201. if (type === 'number' || type === 'BigNumber' || type === 'Complex') {
  202. return type;
  203. }
  204. var hasNumber = false;
  205. var hasBig = false;
  206. var hasComplex = false;
  207. for (var i = 0; i < N; i++) {
  208. for (var j = 0; j < N; j++) {
  209. var el = arr[i][j];
  210. if (isNumber(el) || isFraction(el)) {
  211. hasNumber = true;
  212. } else if (isBigNumber(el)) {
  213. hasBig = true;
  214. } else if (isComplex(el)) {
  215. hasComplex = true;
  216. } else {
  217. throw TypeError('Unsupported type in Matrix: ' + typeOf(el));
  218. }
  219. }
  220. }
  221. if (hasBig && hasComplex) {
  222. console.warn('Complex BigNumbers not supported, this operation will lose precission.');
  223. }
  224. if (hasComplex) {
  225. for (var _i = 0; _i < N; _i++) {
  226. for (var _j = 0; _j < N; _j++) {
  227. arr[_i][_j] = complex(arr[_i][_j]);
  228. }
  229. }
  230. return 'Complex';
  231. }
  232. if (hasBig) {
  233. for (var _i2 = 0; _i2 < N; _i2++) {
  234. for (var _j2 = 0; _j2 < N; _j2++) {
  235. arr[_i2][_j2] = bignumber(arr[_i2][_j2]);
  236. }
  237. }
  238. return 'BigNumber';
  239. }
  240. if (hasNumber) {
  241. for (var _i3 = 0; _i3 < N; _i3++) {
  242. for (var _j3 = 0; _j3 < N; _j3++) {
  243. arr[_i3][_j3] = number(arr[_i3][_j3]);
  244. }
  245. }
  246. return 'number';
  247. } else {
  248. throw TypeError('Matrix contains unsupported types only.');
  249. }
  250. }
  251. });