complexEigs.js (24658B)
1 "use strict"; 2 3 var _interopRequireDefault = require("@babel/runtime/helpers/interopRequireDefault"); 4 5 Object.defineProperty(exports, "__esModule", { 6 value: true 7 }); 8 exports.createComplexEigs = createComplexEigs; 9 10 var _toConsumableArray2 = _interopRequireDefault(require("@babel/runtime/helpers/toConsumableArray")); 11 12 var _object = require("../../../utils/object.js"); 13 14 function _createForOfIteratorHelper(o, allowArrayLike) { var it = typeof Symbol !== "undefined" && o[Symbol.iterator] || o["@@iterator"]; if (!it) { if (Array.isArray(o) || (it = _unsupportedIterableToArray(o)) || allowArrayLike && o && typeof o.length === "number") { if (it) o = it; var i = 0; var F = function F() {}; return { s: F, n: function n() { if (i >= o.length) return { done: true }; return { done: false, value: o[i++] }; }, e: function e(_e) { throw _e; }, f: F }; } throw new TypeError("Invalid attempt to iterate non-iterable instance.\nIn order to be iterable, non-array objects must have a [Symbol.iterator]() method."); } var normalCompletion = true, didErr = false, err; return { s: function s() { it = it.call(o); }, n: function n() { var step = it.next(); normalCompletion = step.done; return step; }, e: function e(_e2) { didErr = true; err = _e2; }, f: function f() { try { if (!normalCompletion && it.return != null) it.return(); } finally { if (didErr) throw err; } } }; } 15 16 function _unsupportedIterableToArray(o, minLen) { if (!o) return; if (typeof o === "string") return _arrayLikeToArray(o, minLen); var n = Object.prototype.toString.call(o).slice(8, -1); if (n === "Object" && o.constructor) n = o.constructor.name; if (n === "Map" || n === "Set") return Array.from(o); if (n === "Arguments" || /^(?:Ui|I)nt(?:8|16|32)(?:Clamped)?Array$/.test(n)) return _arrayLikeToArray(o, minLen); } 17 18 function _arrayLikeToArray(arr, len) { if (len == null || len > arr.length) len = arr.length; for (var i = 0, arr2 = new Array(len); i < len; i++) { arr2[i] = arr[i]; } return arr2; } 19 20 function createComplexEigs(_ref) { 21 var addScalar = _ref.addScalar, 22 subtract = _ref.subtract, 23 flatten = _ref.flatten, 24 multiply = _ref.multiply, 25 multiplyScalar = _ref.multiplyScalar, 26 divideScalar = _ref.divideScalar, 27 sqrt = _ref.sqrt, 28 abs = _ref.abs, 29 bignumber = _ref.bignumber, 30 diag = _ref.diag, 31 inv = _ref.inv, 32 qr = _ref.qr, 33 usolve = _ref.usolve, 34 usolveAll = _ref.usolveAll, 35 equal = _ref.equal, 36 complex = _ref.complex, 37 larger = _ref.larger, 38 smaller = _ref.smaller, 39 matrixFromColumns = _ref.matrixFromColumns, 40 dot = _ref.dot; 41 42 /** 43 * @param {number[][]} arr the matrix to find eigenvalues of 44 * @param {number} N size of the matrix 45 * @param {number|BigNumber} prec precision, anything lower will be considered zero 46 * @param {'number'|'BigNumber'|'Complex'} type 47 * @param {boolean} findVectors should we find eigenvectors? 48 * 49 * @returns {{ values: number[], vectors: number[][] }} 50 */ 51 function complexEigs(arr, N, prec, type, findVectors) { 52 if (findVectors === undefined) { 53 findVectors = true; 54 } // TODO check if any row/col are zero except the diagonal 55 // make sure corresponding rows and columns have similar magnitude 56 // important because of numerical stability 57 // MODIFIES arr by side effect! 58 59 60 var R = balance(arr, N, prec, type, findVectors); // R is the row transformation matrix 61 // arr = A' = R A R⁻¹, A is the original matrix 62 // (if findVectors is false, R is undefined) 63 // (And so to return to original matrix: A = R⁻¹ arr R) 64 // TODO if magnitudes of elements vary over many orders, 65 // move greatest elements to the top left corner 66 // using similarity transformations, reduce the matrix 67 // to Hessenberg form (upper triangular plus one subdiagonal row) 68 // updates the transformation matrix R with new row operationsq 69 // MODIFIES arr by side effect! 70 71 reduceToHessenberg(arr, N, prec, type, findVectors, R); // still true that original A = R⁻¹ arr R) 72 // find eigenvalues 73 74 var _iterateUntilTriangul = iterateUntilTriangular(arr, N, prec, type, findVectors), 75 values = _iterateUntilTriangul.values, 76 C = _iterateUntilTriangul.C; // values is the list of eigenvalues, C is the column 77 // transformation matrix that transforms arr, the hessenberg 78 // matrix, to upper triangular 79 // (So U = C⁻¹ arr C and the relationship between current arr 80 // and original A is unchanged.) 81 82 83 var vectors; 84 85 if (findVectors) { 86 vectors = findEigenvectors(arr, N, C, R, values, prec, type); 87 vectors = matrixFromColumns.apply(void 0, (0, _toConsumableArray2.default)(vectors)); 88 } 89 90 return { 91 values: values, 92 vectors: vectors 93 }; 94 } 95 /** 96 * @param {number[][]} arr 97 * @param {number} N 98 * @param {number} prec 99 * @param {'number'|'BigNumber'|'Complex'} type 100 * @returns {number[][]} 101 */ 102 103 104 function balance(arr, N, prec, type, findVectors) { 105 var big = type === 'BigNumber'; 106 var cplx = type === 'Complex'; 107 var realzero = big ? bignumber(0) : 0; 108 var one = big ? bignumber(1) : cplx ? complex(1) : 1; 109 var realone = big ? bignumber(1) : 1; // base of the floating-point arithmetic 110 111 var radix = big ? bignumber(10) : 2; 112 var radixSq = multiplyScalar(radix, radix); // the diagonal transformation matrix R 113 114 var Rdiag; 115 116 if (findVectors) { 117 Rdiag = Array(N).fill(one); 118 } // this isn't the only time we loop thru the matrix... 119 120 121 var last = false; 122 123 while (!last) { 124 // ...haha I'm joking! unless... 125 last = true; 126 127 for (var i = 0; i < N; i++) { 128 // compute the taxicab norm of i-th column and row 129 // TODO optimize for complex numbers 130 var colNorm = realzero; 131 var rowNorm = realzero; 132 133 for (var j = 0; j < N; j++) { 134 if (i === j) continue; 135 var c = abs(arr[i][j]); // should be real 136 137 colNorm = addScalar(colNorm, c); 138 rowNorm = addScalar(rowNorm, c); 139 } 140 141 if (!equal(colNorm, 0) && !equal(rowNorm, 0)) { 142 // find integer power closest to balancing the matrix 143 // (we want to scale only by integer powers of radix, 144 // so that we don't lose any precision due to round-off) 145 var f = realone; 146 var _c = colNorm; 147 var rowDivRadix = divideScalar(rowNorm, radix); 148 var rowMulRadix = multiplyScalar(rowNorm, radix); 149 150 while (smaller(_c, rowDivRadix)) { 151 _c = multiplyScalar(_c, radixSq); 152 f = multiplyScalar(f, radix); 153 } 154 155 while (larger(_c, rowMulRadix)) { 156 _c = divideScalar(_c, radixSq); 157 f = divideScalar(f, radix); 158 } // check whether balancing is needed 159 // condition = (c + rowNorm) / f < 0.95 * (colNorm + rowNorm) 160 161 162 var condition = smaller(divideScalar(addScalar(_c, rowNorm), f), multiplyScalar(addScalar(colNorm, rowNorm), 0.95)); // apply balancing similarity transformation 163 164 if (condition) { 165 // we should loop once again to check whether 166 // another rebalancing is needed 167 last = false; 168 var g = divideScalar(1, f); 169 170 for (var _j = 0; _j < N; _j++) { 171 if (i === _j) { 172 continue; 173 } 174 175 arr[i][_j] = multiplyScalar(arr[i][_j], f); 176 arr[_j][i] = multiplyScalar(arr[_j][i], g); 177 } // keep track of transformations 178 179 180 if (findVectors) { 181 Rdiag[i] = multiplyScalar(Rdiag[i], f); 182 } 183 } 184 } 185 } 186 } // return the diagonal row transformation matrix 187 188 189 return diag(Rdiag); 190 } 191 /** 192 * @param {number[][]} arr 193 * @param {number} N 194 * @param {number} prec 195 * @param {'number'|'BigNumber'|'Complex'} type 196 * @param {boolean} findVectors 197 * @param {number[][]} R the row transformation matrix that will be modified 198 */ 199 200 201 function reduceToHessenberg(arr, N, prec, type, findVectors, R) { 202 var big = type === 'BigNumber'; 203 var cplx = type === 'Complex'; 204 var zero = big ? bignumber(0) : cplx ? complex(0) : 0; 205 206 if (big) { 207 prec = bignumber(prec); 208 } 209 210 for (var i = 0; i < N - 2; i++) { 211 // Find the largest subdiag element in the i-th col 212 var maxIndex = 0; 213 var max = zero; 214 215 for (var j = i + 1; j < N; j++) { 216 var el = arr[j][i]; 217 218 if (smaller(abs(max), abs(el))) { 219 max = el; 220 maxIndex = j; 221 } 222 } // This col is pivoted, no need to do anything 223 224 225 if (smaller(abs(max), prec)) { 226 continue; 227 } 228 229 if (maxIndex !== i + 1) { 230 // Interchange maxIndex-th and (i+1)-th row 231 var tmp1 = arr[maxIndex]; 232 arr[maxIndex] = arr[i + 1]; 233 arr[i + 1] = tmp1; // Interchange maxIndex-th and (i+1)-th column 234 235 for (var _j2 = 0; _j2 < N; _j2++) { 236 var tmp2 = arr[_j2][maxIndex]; 237 arr[_j2][maxIndex] = arr[_j2][i + 1]; 238 arr[_j2][i + 1] = tmp2; 239 } // keep track of transformations 240 241 242 if (findVectors) { 243 var tmp3 = R[maxIndex]; 244 R[maxIndex] = R[i + 1]; 245 R[i + 1] = tmp3; 246 } 247 } // Reduce following rows and columns 248 249 250 for (var _j3 = i + 2; _j3 < N; _j3++) { 251 var n = divideScalar(arr[_j3][i], max); 252 253 if (n === 0) { 254 continue; 255 } // from j-th row subtract n-times (i+1)th row 256 257 258 for (var k = 0; k < N; k++) { 259 arr[_j3][k] = subtract(arr[_j3][k], multiplyScalar(n, arr[i + 1][k])); 260 } // to (i+1)th column add n-times j-th column 261 262 263 for (var _k = 0; _k < N; _k++) { 264 arr[_k][i + 1] = addScalar(arr[_k][i + 1], multiplyScalar(n, arr[_k][_j3])); 265 } // keep track of transformations 266 267 268 if (findVectors) { 269 for (var _k2 = 0; _k2 < N; _k2++) { 270 R[_j3][_k2] = subtract(R[_j3][_k2], multiplyScalar(n, R[i + 1][_k2])); 271 } 272 } 273 } 274 } 275 276 return R; 277 } 278 /** 279 * @returns {{values: values, C: Matrix}} 280 * @see Press, Wiliams: Numerical recipes in Fortran 77 281 * @see https://en.wikipedia.org/wiki/QR_algorithm 282 */ 283 284 285 function iterateUntilTriangular(A, N, prec, type, findVectors) { 286 var big = type === 'BigNumber'; 287 var cplx = type === 'Complex'; 288 var one = big ? bignumber(1) : cplx ? complex(1) : 1; 289 290 if (big) { 291 prec = bignumber(prec); 292 } // The Francis Algorithm 293 // The core idea of this algorithm is that doing successive 294 // A' = Q⁺AQ transformations will eventually converge to block- 295 // upper-triangular with diagonal blocks either 1x1 or 2x2. 296 // The Q here is the one from the QR decomposition, A = QR. 297 // Since the eigenvalues of a block-upper-triangular matrix are 298 // the eigenvalues of its diagonal blocks and we know how to find 299 // eigenvalues of a 2x2 matrix, we know the eigenvalues of A. 300 301 302 var arr = (0, _object.clone)(A); // the list of converged eigenvalues 303 304 var lambdas = []; // size of arr, which will get smaller as eigenvalues converge 305 306 var n = N; // the diagonal of the block-diagonal matrix that turns 307 // converged 2x2 matrices into upper triangular matrices 308 309 var Sdiag = []; // N×N matrix describing the overall transformation done during the QR algorithm 310 311 var Qtotal = findVectors ? diag(Array(N).fill(one)) : undefined; // n×n matrix describing the QR transformations done since last convergence 312 313 var Qpartial = findVectors ? diag(Array(n).fill(one)) : undefined; // last eigenvalue converged before this many steps 314 315 var lastConvergenceBefore = 0; 316 317 while (lastConvergenceBefore <= 100) { 318 lastConvergenceBefore += 1; // TODO if the convergence is slow, do something clever 319 // Perform the factorization 320 321 var k = 0; // TODO set close to an eigenvalue 322 323 for (var i = 0; i < n; i++) { 324 arr[i][i] = subtract(arr[i][i], k); 325 } // TODO do an implicit QR transformation 326 327 328 var _qr = qr(arr), 329 Q = _qr.Q, 330 R = _qr.R; 331 332 arr = multiply(R, Q); 333 334 for (var _i = 0; _i < n; _i++) { 335 arr[_i][_i] = addScalar(arr[_i][_i], k); 336 } // keep track of transformations 337 338 339 if (findVectors) { 340 Qpartial = multiply(Qpartial, Q); 341 } // The rightmost diagonal element converged to an eigenvalue 342 343 344 if (n === 1 || smaller(abs(arr[n - 1][n - 2]), prec)) { 345 lastConvergenceBefore = 0; 346 lambdas.push(arr[n - 1][n - 1]); // keep track of transformations 347 348 if (findVectors) { 349 Sdiag.unshift([[1]]); 350 inflateMatrix(Qpartial, N); 351 Qtotal = multiply(Qtotal, Qpartial); 352 353 if (n > 1) { 354 Qpartial = diag(Array(n - 1).fill(one)); 355 } 356 } // reduce the matrix size 357 358 359 n -= 1; 360 arr.pop(); 361 362 for (var _i2 = 0; _i2 < n; _i2++) { 363 arr[_i2].pop(); 364 } // The rightmost diagonal 2x2 block converged 365 366 } else if (n === 2 || smaller(abs(arr[n - 2][n - 3]), prec)) { 367 lastConvergenceBefore = 0; 368 var ll = eigenvalues2x2(arr[n - 2][n - 2], arr[n - 2][n - 1], arr[n - 1][n - 2], arr[n - 1][n - 1]); 369 lambdas.push.apply(lambdas, (0, _toConsumableArray2.default)(ll)); // keep track of transformations 370 371 if (findVectors) { 372 Sdiag.unshift(jordanBase2x2(arr[n - 2][n - 2], arr[n - 2][n - 1], arr[n - 1][n - 2], arr[n - 1][n - 1], ll[0], ll[1], prec, type)); 373 inflateMatrix(Qpartial, N); 374 Qtotal = multiply(Qtotal, Qpartial); 375 376 if (n > 2) { 377 Qpartial = diag(Array(n - 2).fill(one)); 378 } 379 } // reduce the matrix size 380 381 382 n -= 2; 383 arr.pop(); 384 arr.pop(); 385 386 for (var _i3 = 0; _i3 < n; _i3++) { 387 arr[_i3].pop(); 388 389 arr[_i3].pop(); 390 } 391 } 392 393 if (n === 0) { 394 break; 395 } 396 } // standard sorting 397 398 399 lambdas.sort(function (a, b) { 400 return +subtract(abs(a), abs(b)); 401 }); // the algorithm didn't converge 402 403 if (lastConvergenceBefore > 100) { 404 var err = Error('The eigenvalues failed to converge. Only found these eigenvalues: ' + lambdas.join(', ')); 405 err.values = lambdas; 406 err.vectors = []; 407 throw err; 408 } // combine the overall QR transformation Qtotal with the subsequent 409 // transformation S that turns the diagonal 2x2 blocks to upper triangular 410 411 412 var C = findVectors ? multiply(Qtotal, blockDiag(Sdiag, N)) : undefined; 413 return { 414 values: lambdas, 415 C: C 416 }; 417 } 418 /** 419 * @param {Matrix} A hessenberg-form matrix 420 * @param {number} N size of A 421 * @param {Matrix} C column transformation matrix that turns A into upper triangular 422 * @param {Matrix} R similarity that turns original matrix into A 423 * @param {number[]} values array of eigenvalues of A 424 * @param {'number'|'BigNumber'|'Complex'} type 425 * @returns {number[][]} eigenvalues 426 */ 427 428 429 function findEigenvectors(A, N, C, R, values, prec, type) { 430 var Cinv = inv(C); 431 var U = multiply(Cinv, A, C); 432 var big = type === 'BigNumber'; 433 var cplx = type === 'Complex'; 434 var zero = big ? bignumber(0) : cplx ? complex(0) : 0; 435 var one = big ? bignumber(1) : cplx ? complex(1) : 1; // turn values into a kind of "multiset" 436 // this way it is easier to find eigenvectors 437 438 var uniqueValues = []; 439 var multiplicities = []; 440 441 var _iterator = _createForOfIteratorHelper(values), 442 _step; 443 444 try { 445 for (_iterator.s(); !(_step = _iterator.n()).done;) { 446 var λ = _step.value; 447 448 var _i4 = indexOf(uniqueValues, λ, equal); 449 450 if (_i4 === -1) { 451 uniqueValues.push(λ); 452 multiplicities.push(1); 453 } else { 454 multiplicities[_i4] += 1; 455 } 456 } // find eigenvectors by solving U − λE = 0 457 // TODO replace with an iterative eigenvector algorithm 458 // (this one might fail for imprecise eigenvalues) 459 460 } catch (err) { 461 _iterator.e(err); 462 } finally { 463 _iterator.f(); 464 } 465 466 var vectors = []; 467 var len = uniqueValues.length; 468 var b = Array(N).fill(zero); 469 var E = diag(Array(N).fill(one)); // eigenvalues for which usolve failed (due to numerical error) 470 471 var failedLambdas = []; 472 473 var _loop = function _loop(i) { 474 var λ = uniqueValues[i]; 475 var S = subtract(U, multiply(λ, E)); // the characteristic matrix 476 477 var solutions = usolveAll(S, b); 478 solutions.shift(); // ignore the null vector 479 // looks like we missed something, try inverse iteration 480 481 while (solutions.length < multiplicities[i]) { 482 var approxVec = inverseIterate(S, N, solutions, prec, type); 483 484 if (approxVec == null) { 485 // no more vectors were found 486 failedLambdas.push(λ); 487 break; 488 } 489 490 solutions.push(approxVec); 491 } // Transform back into original array coordinates 492 493 494 var correction = multiply(inv(R), C); 495 solutions = solutions.map(function (v) { 496 return multiply(correction, v); 497 }); 498 vectors.push.apply(vectors, (0, _toConsumableArray2.default)(solutions.map(function (v) { 499 return flatten(v); 500 }))); 501 }; 502 503 for (var i = 0; i < len; i++) { 504 _loop(i); 505 } 506 507 if (failedLambdas.length !== 0) { 508 var err = new Error('Failed to find eigenvectors for the following eigenvalues: ' + failedLambdas.join(', ')); 509 err.values = values; 510 err.vectors = vectors; 511 throw err; 512 } 513 514 return vectors; 515 } 516 /** 517 * Compute the eigenvalues of an 2x2 matrix 518 * @return {[number,number]} 519 */ 520 521 522 function eigenvalues2x2(a, b, c, d) { 523 // λ± = ½ trA ± ½ √( tr²A - 4 detA ) 524 var trA = addScalar(a, d); 525 var detA = subtract(multiplyScalar(a, d), multiplyScalar(b, c)); 526 var x = multiplyScalar(trA, 0.5); 527 var y = multiplyScalar(sqrt(subtract(multiplyScalar(trA, trA), multiplyScalar(4, detA))), 0.5); 528 return [addScalar(x, y), subtract(x, y)]; 529 } 530 /** 531 * For an 2x2 matrix compute the transformation matrix S, 532 * so that SAS⁻¹ is an upper triangular matrix 533 * @return {[[number,number],[number,number]]} 534 * @see https://math.berkeley.edu/~ogus/old/Math_54-05/webfoils/jordan.pdf 535 * @see http://people.math.harvard.edu/~knill/teaching/math21b2004/exhibits/2dmatrices/index.html 536 */ 537 538 539 function jordanBase2x2(a, b, c, d, l1, l2, prec, type) { 540 var big = type === 'BigNumber'; 541 var cplx = type === 'Complex'; 542 var zero = big ? bignumber(0) : cplx ? complex(0) : 0; 543 var one = big ? bignumber(1) : cplx ? complex(1) : 1; // matrix is already upper triangular 544 // return an identity matrix 545 546 if (smaller(abs(c), prec)) { 547 return [[one, zero], [zero, one]]; 548 } // matrix is diagonalizable 549 // return its eigenvectors as columns 550 551 552 if (larger(abs(subtract(l1, l2)), prec)) { 553 return [[subtract(l1, d), subtract(l2, d)], [c, c]]; 554 } // matrix is not diagonalizable 555 // compute off-diagonal elements of N = A - λI 556 // N₁₂ = 0 ⇒ S = ( N⃗₁, I⃗₁ ) 557 // N₁₂ ≠ 0 ⇒ S = ( N⃗₂, I⃗₂ ) 558 559 560 var na = subtract(a, l1); 561 var nb = subtract(b, l1); 562 var nc = subtract(c, l1); 563 var nd = subtract(d, l1); 564 565 if (smaller(abs(nb), prec)) { 566 return [[na, one], [nc, zero]]; 567 } else { 568 return [[nb, zero], [nd, one]]; 569 } 570 } 571 /** 572 * Enlarge the matrix from n×n to N×N, setting the new 573 * elements to 1 on diagonal and 0 elsewhere 574 */ 575 576 577 function inflateMatrix(arr, N) { 578 // add columns 579 for (var i = 0; i < arr.length; i++) { 580 var _arr$i; 581 582 (_arr$i = arr[i]).push.apply(_arr$i, (0, _toConsumableArray2.default)(Array(N - arr[i].length).fill(0))); 583 } // add rows 584 585 586 for (var _i5 = arr.length; _i5 < N; _i5++) { 587 arr.push(Array(N).fill(0)); 588 arr[_i5][_i5] = 1; 589 } 590 591 return arr; 592 } 593 /** 594 * Create a block-diagonal matrix with the given square matrices on the diagonal 595 * @param {Matrix[] | number[][][]} arr array of matrices to be placed on the diagonal 596 * @param {number} N the size of the resulting matrix 597 */ 598 599 600 function blockDiag(arr, N) { 601 var M = []; 602 603 for (var i = 0; i < N; i++) { 604 M[i] = Array(N).fill(0); 605 } 606 607 var I = 0; 608 609 var _iterator2 = _createForOfIteratorHelper(arr), 610 _step2; 611 612 try { 613 for (_iterator2.s(); !(_step2 = _iterator2.n()).done;) { 614 var sub = _step2.value; 615 var n = sub.length; 616 617 for (var _i6 = 0; _i6 < n; _i6++) { 618 for (var j = 0; j < n; j++) { 619 M[I + _i6][I + j] = sub[_i6][j]; 620 } 621 } 622 623 I += n; 624 } 625 } catch (err) { 626 _iterator2.e(err); 627 } finally { 628 _iterator2.f(); 629 } 630 631 return M; 632 } 633 /** 634 * Finds the index of an element in an array using a custom equality function 635 * @template T 636 * @param {Array<T>} arr array in which to search 637 * @param {T} el the element to find 638 * @param {function(T, T): boolean} fn the equality function, first argument is an element of `arr`, the second is always `el` 639 * @returns {number} the index of `el`, or -1 when it's not in `arr` 640 */ 641 642 643 function indexOf(arr, el, fn) { 644 for (var i = 0; i < arr.length; i++) { 645 if (fn(arr[i], el)) { 646 return i; 647 } 648 } 649 650 return -1; 651 } 652 /** 653 * Provided a near-singular upper-triangular matrix A and a list of vectors, 654 * finds an eigenvector of A with the smallest eigenvalue, which is orthogonal 655 * to each vector in the list 656 * @template T 657 * @param {T[][]} A near-singular square matrix 658 * @param {number} N dimension 659 * @param {T[][]} orthog list of vectors 660 * @param {number} prec epsilon 661 * @param {'number'|'BigNumber'|'Complex'} type 662 * @return {T[] | null} eigenvector 663 * 664 * @see Numerical Recipes for Fortran 77 – 11.7 Eigenvalues or Eigenvectors by Inverse Iteration 665 */ 666 667 668 function inverseIterate(A, N, orthog, prec, type) { 669 var largeNum = type === 'BigNumber' ? bignumber(1000) : 1000; 670 var b; // the vector 671 // you better choose a random vector before I count to five 672 673 var i = 0; 674 675 while (true) { 676 b = randomOrthogonalVector(N, orthog, type); 677 b = usolve(A, b); 678 679 if (larger(norm(b), largeNum)) { 680 break; 681 } 682 683 if (++i >= 5) { 684 return null; 685 } 686 } // you better converge before I count to ten 687 688 689 i = 0; 690 691 while (true) { 692 var c = usolve(A, b); 693 694 if (smaller(norm(orthogonalComplement(b, [c])), prec)) { 695 break; 696 } 697 698 if (++i >= 10) { 699 return null; 700 } 701 702 b = normalize(c); 703 } 704 705 return b; 706 } 707 /** 708 * Generates a random unit vector of dimension N, orthogonal to each vector in the list 709 * @template T 710 * @param {number} N dimension 711 * @param {T[][]} orthog list of vectors 712 * @param {'number'|'BigNumber'|'Complex'} type 713 * @returns {T[]} random vector 714 */ 715 716 717 function randomOrthogonalVector(N, orthog, type) { 718 var big = type === 'BigNumber'; 719 var cplx = type === 'Complex'; // generate random vector with the correct type 720 721 var v = Array(N).fill(0).map(function (_) { 722 return 2 * Math.random() - 1; 723 }); 724 725 if (big) { 726 v = v.map(function (n) { 727 return bignumber(n); 728 }); 729 } 730 731 if (cplx) { 732 v = v.map(function (n) { 733 return complex(n); 734 }); 735 } // project to orthogonal complement 736 737 738 v = orthogonalComplement(v, orthog); // normalize 739 740 return normalize(v, type); 741 } 742 /** 743 * Project vector v to the orthogonal complement of an array of vectors 744 */ 745 746 747 function orthogonalComplement(v, orthog) { 748 var _iterator3 = _createForOfIteratorHelper(orthog), 749 _step3; 750 751 try { 752 for (_iterator3.s(); !(_step3 = _iterator3.n()).done;) { 753 var w = _step3.value; 754 // v := v − (w, v)/∥w∥² w 755 v = subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w)); 756 } 757 } catch (err) { 758 _iterator3.e(err); 759 } finally { 760 _iterator3.f(); 761 } 762 763 return v; 764 } 765 /** 766 * Calculate the norm of a vector. 767 * We can't use math.norm because factory can't handle circular dependency. 768 * Seriously, I'm really fed up with factory. 769 */ 770 771 772 function norm(v) { 773 return abs(sqrt(dot(v, v))); 774 } 775 /** 776 * Normalize a vector 777 * @template T 778 * @param {T[]} v 779 * @param {'number'|'BigNumber'|'Complex'} type 780 * @returns {T[]} normalized vec 781 */ 782 783 784 function normalize(v, type) { 785 var big = type === 'BigNumber'; 786 var cplx = type === 'Complex'; 787 var one = big ? bignumber(1) : cplx ? complex(1) : 1; 788 return multiply(divideScalar(one, norm(v)), v); 789 } 790 791 return complexEigs; 792 }