Repository: josdejong/mathjs
Pull request #1086 (2018-04-18 22:33:01)

Adds matrix function math.sqrtm(A)
## Summary

Adds the `math.sqrtm(A)` matrix function to compute the principal square root of a square matrix. This is similar to MATLAB's [`sqrtm`](https://uk.mathworks.com/help/matlab/ref/sqrtm.html).

This PR is somewhat related to #263.

#### Denman–Beavers method

The implementation being current used is based on the [Denman–Beavers iterative method](https://en.wikipedia.org/wiki/Square_root_of_a_matrix#By_Denman%E2%80%93Beavers_iteration):

![Denman–Beavers Y_{k+1}](https://latex.codecogs.com/svg.latex?Y_%7Bk+1%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%20%28%20Y_k%20+%20Z%5E%7B-1%7D_k%20%29)

![Denman–Beavers Z_{k+1}](https://latex.codecogs.com/svg.latex?Z_%7Bk+1%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%20%28%20Z_k%20+%20Y%5E%7B-1%7D_k%20%29)

#### Babylonian method

The [Babylonian iterative method](https://en.wikipedia.org/wiki/Square_root_of_a_matrix#By_the_Babylonian_method) is also implemented:

![Babylonian method](https://latex.codecogs.com/svg.latex?X_%7Bk+1%7D%20%3D%20%5Cfrac%7B1%7D%7B2%7D%20%28%20X_k%20+%20A%20X%5E%7B-1%7D_k%20%29)

> *"However, unlike Denman–Beavers iteration, the Babylonian method is numerically unstable and more likely to fail to converge."* [2]

---

#### References
1. https://en.wikipedia.org/wiki/Square_root_of_a_matrix
2. Higham, N.J., 1986. [Newton’s method for the matrix square root](http://www.ams.org/journals/mcom/1986-46-174/S0025-5718-1986-0829624-5/S0025-5718-1986-0829624-5.pdf). *Mathematics of Computation*, 46(174), pp.537-549.

ferrolho: @harrysarson I think all of your suggestions have now been addressed. I do agree they were quite relevant! Thanks! :+1: 
ferrolho: That is correct, `_denmanBeavers` is preferable to `_babylonianMethod` for numerical stability.

I started to implement `_babylonianMethod` before `_denmanBeavers` so I thought *"Why not include both implementations?"*
ferrolho: Up to you guys. :stuck_out_tongue: 
josdejong: Thanks Henrique, your implementation looks really neat, nice job :+1: 

Thanks @harrysarson for reviewing the PR.

@ferrolho a few more feedback points:

1. > If we are not using _babylonianMethod does it make sense to remove it from the source code? It would make maintaining sqrtm simpler going forward.
    
    Agree, if `_denmanBeavers` is simply the best then let's keep simple and just keep this implementation.

2. Right now you can input a Matrix or an Array, and the output is always a Matrix. I think it would be good to change this behavior such that for Matrix input, the output is a Matrix, and for Array inputs, the output will be an Array, so the behavior is consistent with other functions in math.js.

3. The function is well tested. It would be neat to add one more test: testing a matrix containing BigNumbers instead of regular numbers. Looking at the implementation I think it just works, but it would be good to validate this.
ChristopherChudzicki: Two pretty minor questions:
1. Right now the tolerance is hard-coded as `1e-6`, right? Is it worth having a configurable tolerance, with `1e-6` as the default?
2. In the `toTex()` method, do you really want the outer parentheses? I'm not sure what's standard typographically, but `math.parse("[[1,2],[3,4]]").toTex()` yields `{\begin{bmatrix}1&2\\3&4\\\end{bmatrix}}^{2}` not `\left( \begin{bmatrix}1&2\\3&4\\\end{bmatrix} \right)^{2}`. (That is, squaring has no outer parentheses.)
ferrolho: @josdejong I believe the latest commit addressed both points `1` and `3`. Regarding point `2`, I am not so sure if I understood your feedback, as this method is directed towards matrices. For element-wise square root we already have `math.sqrtm(x)`.

I have also edited the LaTeX conversion as pointed out by @ChristopherChudzicki.
josdejong: Thanks for your updates @ferrolho . I indeed do not mean element-wise computation, sorry for the confusion. The behavior I would expect is basically this:

```js
// Array in -> Array out
var a_in = [[33, 24], [48, 57]]
var a_out = math.sqrtm(a_in)   // outputs a 2D Array [[5,  2], [ 4,  7]]
math.typeof(a_in)              // outputs String 'Array'
math.typeof(a_out)             // outputs String 'Array'

// Matrix in -> Matrix out
var m_in = math.matrix([[33, 24], [48, 57]])
var m_out = math.sqrtm(m_in)   // outputs a 2D Matrix [[5,  2], [ 4,  7]]
math.typeof(m_in)              // outputs String 'Matrix'
math.typeof(m_out)             // outputs String 'Matrix'
```
Does that make sense?
ferrolho: @josdejong it does make sense. My bad for not understanding it at first...

The last commit should ensure data type persistency. Thanks for your detailed previous comment.
josdejong: @ferrolho thanks, looks good! that is even less changes in the code needed than I expected :D 

I see one more issue which is important: in the unit tests, at a number of tests `math.deepEqual` is used, which should be `assert.deepEqual`. Currently, these tests do not throw an error at all if the actual result does not match the expected result.
ferrolho: @josdejong I believe I have fixed the issue.

However, I do not seem to be able to get the BigNumber unit test to pass — *cf.* the comment above. Am I doing something stupid again? :upside_down_face: 
josdejong: Thanks @ferrolho, nice job :+1:

Patch:
diff --git a/lib/function/matrix/index.js b/lib/function/matrix/index.js
index e08d885c8..0a30c8a93 100644
--- a/lib/function/matrix/index.js
+++ b/lib/function/matrix/index.js
@@ -18,6 +18,7 @@ module.exports = [
   require('./resize'),
   require('./size'),
   require('./sort'),
+  require('./sqrtm'),
   require('./squeeze'),
   require('./subset'),
   require('./trace'),
diff --git a/lib/function/matrix/sqrtm.js b/lib/function/matrix/sqrtm.js
new file mode 100644
index 000000000..acd0a03bf
--- /dev/null
+++ b/lib/function/matrix/sqrtm.js
@@ -0,0 +1,110 @@
+'use strict';
+
+var array  = require('../../utils/array');
+var latex  = require('../../utils/latex');
+var string = require('../../utils/string');
+
+function factory(type, config, load, typed) {
+  var matrix   = load(require('../../type/matrix/function/matrix'));
+  var abs      = load(require('../arithmetic/abs'));
+  var add      = load(require('../arithmetic/add'));
+  var divide   = load(require('../arithmetic/divide'));
+  var multiply = load(require('../arithmetic/multiply'));
+  var sqrt     = load(require('../arithmetic/sqrt'));
+  var subtract = load(require('../arithmetic/subtract'));
+  var inv      = load(require('../matrix/inv'));
+  var size     = load(require('../matrix/size'));
+  var max      = load(require('../statistics/max'));
+  var eye      = load(require('./eye'));
+
+  /**
+   * Calculate the principal square root of a square matrix.
+   * The principal square root matrix `X` of another matrix `A` is such that `X * X = A`.
+   *
+   * https://en.wikipedia.org/wiki/Square_root_of_a_matrix
+   * 
+   * Syntax:
+   *
+   *     X = math.sqrtm(A)
+   *
+   * Examples:
+   *
+   *     math.sqrtm([[1, 2], [3, 4]]); // returns [[-2, 1], [1.5, -0.5]]
+   *
+   * See also:
+   *
+   *     sqrt, pow, abs, max
+   *
+   * @param  {Array | Matrix} A   The square matrix `A`
+   * @return {Array | Matrix}     The principal square root of matrix `A`
+   */
+  var sqrtm = typed('sqrtm', {
+    'Array | Matrix': function (A) {
+      var size = type.isMatrix(A) ? A.size() : array.size(A);
+      switch (size.length) {
+        case 1:
+          // Single element Array | Matrix
+          if (size[0] == 1) {
+            return sqrt(A);
+          }
+          else {
+            throw new RangeError('Matrix must be square ' +
+            '(size: ' + string.format(size) + ')');
+          }
+
+        case 2:
+          // Two-dimensional Array | Matrix
+          var rows = size[0];
+          var cols = size[1];
+          if (rows == cols) {
+            return _denmanBeavers(A);
+          }
+          else {
+            throw new RangeError('Matrix must be square ' +
+            '(size: ' + string.format(size) + ')');
+          }
+      }
+    }
+  });
+
+  var _maxIterations = 1e3;
+  var _tolerance = 1e-6;
+
+  /**
+   * Calculate the principal square root matrix using the Denman–Beavers iterative method
+   * 
+   * https://en.wikipedia.org/wiki/Square_root_of_a_matrix#By_Denman–Beavers_iteration
+   * 
+   * @param  {Array | Matrix} A   The square matrix `A`
+   * @return {Array | Matrix}     The principal square root of matrix `A`
+   * @private
+   */
+  function _denmanBeavers(A) {
+    var error;
+    var iterations = 0;
+
+    var Y = A;
+    var Z = eye(size(A));
+
+    do {
+      var Y_k = Y;
+      Y = multiply(0.5, add(Y_k, inv(Z)));
+      Z = multiply(0.5, add(Z, inv(Y_k)));
+
+      error = max(abs(subtract(Y, Y_k)));
+
+      if (error > _tolerance && ++iterations > _maxIterations) {
+        throw new Error('computing square root of matrix: iterative method could not converge');
+      }
+    } while (error > _tolerance);
+
+    return Y;
+  }
+
+  sqrtm.toTex = {1: '{${args[0]}}' + latex.operators['pow'] + '{\\frac{1}{2}}'};
+
+  return sqrtm;
+}
+
+exports.name = 'sqrtm';
+exports.factory = factory;


Test patch:
diff --git a/test/function/matrix/sqrtm.test.js b/test/function/matrix/sqrtm.test.js
new file mode 100644
index 000000000..a686d99a4
--- /dev/null
+++ b/test/function/matrix/sqrtm.test.js
@@ -0,0 +1,72 @@
+// test sqrtm
+var assert = require('assert'),
+    approx = require('../../../tools/approx'),
+    math   = require('../../../index');
+
+describe('sqrtm', function () {
+
+  var A  = [[ 5,  2], [ 4,  7]];
+  var AA = [[33, 24], [48, 57]];
+
+  var B  = [[1.5667, 1.7408], [2.6112, 4.1779]];
+  var BB = [[ 7, 10], [15, 22]];
+
+  it('should return the principal square root of a matrix', function() {
+    approx.deepEqual(math.sqrtm(AA), A);
+    approx.deepEqual(math.sqrtm(BB), B);
+
+    approx.deepEqual(
+      math.sqrtm(
+        [[ 5, -4,  1,  0,  0],
+         [-4,  6, -4,  1,  0],
+         [ 1, -4,  6, -4,  1],
+         [ 0,  1, -4,  6, -4],
+         [ 0,  0,  1, -4,  6]]),
+        [[ 2.001450032980806, -0.997069854672086,  0.004226841220338,  0.004648098208404,  0.003162179891248],
+         [-0.997069854672086,  2.006191486385745, -0.990363307725271,  0.011838782789231,  0.009369460553432],
+         [ 0.004226841220338, -0.990363307725271,  2.017072849046023, -0.974622709753106,  0.026274470491696],
+         [ 0.004648098208404,  0.011838782789231, -0.974622709753106,  2.050268428894304, -0.919971837350421],
+         [ 0.003162179891248,  0.009369460553432,  0.026274470491696, -0.919971837350421,  2.269992000979243]]);
+  });
+
+  it('should return the principal square root of a matrix with just one value', function() {
+    assert.deepEqual(math.sqrtm([ 4]), [2]);
+    assert.deepEqual(math.sqrtm([16]), [4]);
+    assert.deepEqual(math.sqrtm([20.25]), [4.5]);
+  });
+
+  it('should return the principal square root of a matrix of big numbers', function() {
+    assert.deepEqual(math.round(math.sqrtm(math.bignumber(AA)), 20), math.bignumber(A));
+  });
+
+  it('math.pow(math.sqrtm(A), 2) should equal A', function() {
+    approx.deepEqual(math.pow(math.sqrtm(A), 2), A);
+    approx.deepEqual(math.pow(math.sqrtm(B), 2), B);
+    approx.deepEqual(math.pow(math.sqrtm(AA), 2), AA);
+    approx.deepEqual(math.pow(math.sqrtm(BB), 2), BB);
+  });
+
+  it('should throw an error in case of non-square matrices', function() {
+    assert.throws(function () { math.sqrtm([1, 2, 3]) }, /Matrix must be square/);
+    assert.throws(function () { math.sqrtm([[1, 2, 3]]) }, /Matrix must be square/);
+    assert.throws(function () { math.sqrtm([[1, 2, 3], [4, 5, 6]]) }, /Matrix must be square/);
+  });
+
+  it('should LaTeX sqrtm', function () {
+    var expression = math.parse('sqrtm([[33, 24], [48, 57]])');
+    assert.equal(expression.toTex(), '{\\begin{bmatrix}33&24\\\\48&57\\\\\\end{bmatrix}}^{\\frac{1}{2}}');
+  });
+
+  it('should return the result in the same format as the input', function() {
+    assert.equal(math.typeof(math.sqrtm(A)),  'Array');
+    assert.equal(math.typeof(math.sqrtm(B)),  'Array');
+    assert.equal(math.typeof(math.sqrtm(AA)), 'Array');
+    assert.equal(math.typeof(math.sqrtm(BB)), 'Array');
+
+    assert.equal(math.typeof(math.sqrtm(math.matrix(A))),  'Matrix');
+    assert.equal(math.typeof(math.sqrtm(math.matrix(B))),  'Matrix');
+    assert.equal(math.typeof(math.sqrtm(math.matrix(AA))), 'Matrix');
+    assert.equal(math.typeof(math.sqrtm(math.matrix(BB))), 'Matrix');
+  });
+
+});