simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

sqrtm.js (2919B)

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