pinv.js 5.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191
  1. "use strict";
  2. Object.defineProperty(exports, "__esModule", {
  3. value: true
  4. });
  5. exports.createPinv = void 0;
  6. var _is = require("../../utils/is.js");
  7. var _array = require("../../utils/array.js");
  8. var _factory = require("../../utils/factory.js");
  9. var _string = require("../../utils/string.js");
  10. var _object = require("../../utils/object.js");
  11. var name = 'pinv';
  12. var dependencies = ['typed', 'matrix', 'inv', 'deepEqual', 'equal', 'dotDivide', 'dot', 'ctranspose', 'divideScalar', 'multiply', 'add', 'Complex'];
  13. var createPinv = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
  14. var typed = _ref.typed,
  15. matrix = _ref.matrix,
  16. inv = _ref.inv,
  17. deepEqual = _ref.deepEqual,
  18. equal = _ref.equal,
  19. dotDivide = _ref.dotDivide,
  20. dot = _ref.dot,
  21. ctranspose = _ref.ctranspose,
  22. divideScalar = _ref.divideScalar,
  23. multiply = _ref.multiply,
  24. add = _ref.add,
  25. Complex = _ref.Complex;
  26. /**
  27. * Calculate the Moore–Penrose inverse of a matrix.
  28. *
  29. * Syntax:
  30. *
  31. * math.pinv(x)
  32. *
  33. * Examples:
  34. *
  35. * math.pinv([[1, 2], [3, 4]]) // returns [[-2, 1], [1.5, -0.5]]
  36. * math.pinv([[1, 0], [0, 1], [0, 1]]) // returns [[1, 0, 0], [0, 0.5, 0.5]]
  37. * math.pinv(4) // returns 0.25
  38. *
  39. * See also:
  40. *
  41. * inv
  42. *
  43. * @param {number | Complex | Array | Matrix} x Matrix to be inversed
  44. * @return {number | Complex | Array | Matrix} The inverse of `x`.
  45. */
  46. return typed(name, {
  47. 'Array | Matrix': function ArrayMatrix(x) {
  48. var size = (0, _is.isMatrix)(x) ? x.size() : (0, _array.arraySize)(x);
  49. switch (size.length) {
  50. case 1:
  51. // vector
  52. if (_isZeros(x)) return ctranspose(x); // null vector
  53. if (size[0] === 1) {
  54. return inv(x); // invertible matrix
  55. } else {
  56. return dotDivide(ctranspose(x), dot(x, x));
  57. }
  58. case 2:
  59. // two dimensional array
  60. {
  61. if (_isZeros(x)) return ctranspose(x); // zero matrixx
  62. var rows = size[0];
  63. var cols = size[1];
  64. if (rows === cols) {
  65. try {
  66. return inv(x); // invertible matrix
  67. } catch (err) {
  68. if (err instanceof Error && err.message.match(/Cannot calculate inverse, determinant is zero/)) {
  69. // Expected
  70. } else {
  71. throw err;
  72. }
  73. }
  74. }
  75. if ((0, _is.isMatrix)(x)) {
  76. return matrix(_pinv(x.valueOf(), rows, cols), x.storage());
  77. } else {
  78. // return an Array
  79. return _pinv(x, rows, cols);
  80. }
  81. }
  82. default:
  83. // multi dimensional array
  84. throw new RangeError('Matrix must be two dimensional ' + '(size: ' + (0, _string.format)(size) + ')');
  85. }
  86. },
  87. any: function any(x) {
  88. // scalar
  89. if (equal(x, 0)) return (0, _object.clone)(x); // zero
  90. return divideScalar(1, x);
  91. }
  92. });
  93. /**
  94. * Calculate the Moore–Penrose inverse of a matrix
  95. * @param {Array[]} mat A matrix
  96. * @param {number} rows Number of rows
  97. * @param {number} cols Number of columns
  98. * @return {Array[]} pinv Pseudoinverse matrix
  99. * @private
  100. */
  101. function _pinv(mat, rows, cols) {
  102. var _rankFact2 = _rankFact(mat, rows, cols),
  103. C = _rankFact2.C,
  104. F = _rankFact2.F; // TODO: Use SVD instead (may improve precision)
  105. var Cpinv = multiply(inv(multiply(ctranspose(C), C)), ctranspose(C));
  106. var Fpinv = multiply(ctranspose(F), inv(multiply(F, ctranspose(F))));
  107. return multiply(Fpinv, Cpinv);
  108. }
  109. /**
  110. * Calculate the reduced row echelon form of a matrix
  111. *
  112. * Modified from https://rosettacode.org/wiki/Reduced_row_echelon_form
  113. *
  114. * @param {Array[]} mat A matrix
  115. * @param {number} rows Number of rows
  116. * @param {number} cols Number of columns
  117. * @return {Array[]} Reduced row echelon form
  118. * @private
  119. */
  120. function _rref(mat, rows, cols) {
  121. var M = (0, _object.clone)(mat);
  122. var lead = 0;
  123. for (var r = 0; r < rows; r++) {
  124. if (cols <= lead) {
  125. return M;
  126. }
  127. var i = r;
  128. while (_isZero(M[i][lead])) {
  129. i++;
  130. if (rows === i) {
  131. i = r;
  132. lead++;
  133. if (cols === lead) {
  134. return M;
  135. }
  136. }
  137. }
  138. var _ref2 = [M[r], M[i]];
  139. M[i] = _ref2[0];
  140. M[r] = _ref2[1];
  141. var val = M[r][lead];
  142. for (var j = 0; j < cols; j++) {
  143. M[r][j] = dotDivide(M[r][j], val);
  144. }
  145. for (var _i = 0; _i < rows; _i++) {
  146. if (_i === r) continue;
  147. val = M[_i][lead];
  148. for (var _j = 0; _j < cols; _j++) {
  149. M[_i][_j] = add(M[_i][_j], multiply(-1, multiply(val, M[r][_j])));
  150. }
  151. }
  152. lead++;
  153. }
  154. return M;
  155. }
  156. /**
  157. * Calculate the rank factorization of a matrix
  158. *
  159. * @param {Array[]} mat A matrix (M)
  160. * @param {number} rows Number of rows
  161. * @param {number} cols Number of columns
  162. * @return {{C: Array, F: Array}} rank factorization where M = C F
  163. * @private
  164. */
  165. function _rankFact(mat, rows, cols) {
  166. var rref = _rref(mat, rows, cols);
  167. var C = mat.map(function (_, i) {
  168. return _.filter(function (_, j) {
  169. return j < rows && !_isZero(dot(rref[j], rref[j]));
  170. });
  171. });
  172. var F = rref.filter(function (_, i) {
  173. return !_isZero(dot(rref[i], rref[i]));
  174. });
  175. return {
  176. C: C,
  177. F: F
  178. };
  179. }
  180. function _isZero(x) {
  181. return equal(add(x, Complex(1, 1)), add(0, Complex(1, 1)));
  182. }
  183. function _isZeros(arr) {
  184. return deepEqual(add(arr, Complex(1, 1)), add(multiply(arr, 0), Complex(1, 1)));
  185. }
  186. });
  187. exports.createPinv = createPinv;