feat: implement invert for vectors and matrices

Here invert returns the pseudoinverse when the input is not invertible.
  To this end, also implements (rudimentary) conversion of Complex to NumberT,
  an identity method that produces identity matrices, determinant, and
  the adjoint operation on matrices (conjugate transpose).

src/number/type.js

@@ -1,8 +1,10 @@
 import {plain} from './helpers.js'
-import {BooleanT} from '#boolean/BooleanT.js'
-import {Returns} from '#core/Type.js'
+import {NumberT} from './NumberT.js'
+
+import {Returns, ReturnType} from '#core/Type.js'
 import {match} from '#core/TypePatterns.js'
-import {NumberT} from '#number/NumberT.js'
+import {BooleanT} from '#boolean/BooleanT.js'
+import {Complex} from '#complex/Complex.js'
 
 const num = f => Returns(NumberT, f)
 
@@ -10,5 +12,17 @@ export const number = [
    plain(a => a),
    // conversions from Boolean should be consistent with one and zero:
    match(BooleanT, num(p => p ? NumberT.one : NumberT.zero)),
-   match([], num(() => 0))
+   match([], num(() => 0)),
+   match(Complex, (math, C) => {
+      const im = math.im.resolve(C)
+      const re = math.re.resolve(C)
+      const compNum = math.number.resolve(ReturnType(re))
+      const isZ = math.isZero.resolve(ReturnType(im))
+      return num(z => {
+         if (!isZ(im(z))) {
+            throw new RangeError(`can't convert Complex ${z} to number`)
+         }
+         return compNum(re(z))
+      })
+   })
 ]

src/vector/__test__/arithmetic.spec.js

@@ -24,6 +24,30 @@ describe('Vector arithmetic functions', () => {
       assert.deepStrictEqual(
          add([[1, 2], [4, 2]], [0, -1]), [[1, 1], [4, 1]])
    })
+   it('(pseudo)inverts matrices', () => {
+      const inv = math.invert
+      // inverses
+      assert.deepStrictEqual(inv([3, 4, 5]), [3/50, 2/25, 1/10])
+      assert.deepStrictEqual(inv([[4]]), [[0.25]])
+      assert.deepStrictEqual(inv([[5, 2], [-7, -3]]), [[3, 2], [-7, -5]])
+      assert(math.equal(
+         inv([[3, 0, 2], [2, 1, 0], [1, 4, 2]]),
+         [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]]))
+      // pseudoinverses
+      assert.deepStrictEqual(inv([[1, 0], [1, 0]]), [[1/2, 1/2], [0, 0]])
+      assert.deepStrictEqual(
+         inv([[1, 0], [0, 1], [0, 1]]), [[1, 0, 0], [0, 1/2, 1/2]])
+      assert.deepStrictEqual(
+         inv([[1, 0, 0, 0, 2],
+              [0, 0, 3, 0, 0],
+              [0, 0, 0, 0, 0],
+              [0, 4, 0, 0, 0]]),
+         [[1/5,   0, 0,   0],
+          [  0,   0, 0, 1/4],
+          [  0, 1/3, 0,   0],
+          [  0,   0, 0,   0],
+          [2/5,   0, 0,   0]])
+   })
    it('multiplies vectors and matrices', () => {
       const mult = math.multiply
       const pyth = [3, 4, 5]
@@ -38,6 +62,10 @@ describe('Vector arithmetic functions', () => {
       assert.deepStrictEqual(
          mult(mat32, [[1, 2], [3, 4]]),
          [[-8, -10], [-4, -4], [0, 2]])
+      assert(math.equal(
+         mult([[3, 0, 2], [2, 1, 0], [1, 4, 2]],
+              [[1/10, 2/5, -1/10], [-1/5, 1/5, 1/5], [7/20, -3/5, 3/20]]),
+         [[1, 0, 0], [0, 1, 0], [0, 0, 1]]))
    })
    it('negates a vector', () => {
       assert.deepStrictEqual(math.negate([-3, 4, -5]), [3, -4, 5])

src/vector/__test__/type.spec.js

@@ -22,4 +22,24 @@ describe('Vector type functions', () => {
       assert.deepStrictEqual(
          tsp([[1, 2, 3], [4, 5, 6]]), [[1, 4], [2, 5], [3, 6]])
    })
+   it('can take adjoint (conjugate transpose) of a matrix', () => {
+      const cx = math.complex
+      assert.deepStrictEqual(
+         math.adjoint([[cx(1, 1), cx(2, 2)], [cx(3, 3), cx(4, 4)]]),
+         [[cx(1, -1), cx(3, -3)], [cx(2, -2), cx(4, -4)]])
+   })
+   it('generates identity from an example matrix or a number of rows', () => {
+      const id = math.identity
+      const cx = math.complex
+      assert.deepStrictEqual(id(2), [[1, 0], [0, 1]])
+      assert.deepStrictEqual(id(cx(2)), [[cx(1), cx(0)], [cx(0), cx(1)]])
+      assert.deepStrictEqual(
+         id([[1, 2, 3]]),
+         [[1, 0 , 0], [0, 1, 0], [0, 0, 1]])
+   })
+   it('takes the determinant of a matrix', () => {
+      assert.strictEqual(
+         math.determinant([[6, 1, 1], [4, -2, 5], [2, 8, 7]]),
+         -306)
+   })
 })

src/vector/arithmetic.js

@@ -3,7 +3,7 @@ import {
 } from './helpers.js'
 import {Vector} from './Vector.js'
 
-import {ReturnType} from '#core/Type.js'
+import {Returns, ReturnType} from '#core/Type.js'
 import {match} from '#core/TypePatterns.js'
 import {ReturnsAs} from '#generic/helpers.js'
 
@@ -19,6 +19,201 @@ export const abs = promoteUnary('abs')
 export const add = promoteBinary('add')
 
 export const dotMultiply = promoteBinary('multiply')
+export const invert = match(Vector, (math, V, strategy) => {
+   if (V.vectorDepth > 2) {
+      throw new TypeError(
+         'invert not implemented for arrays of dimension > 2')
+   }
+   const normsq = math.normsq.resolve(V, strategy)
+   const NormT = ReturnType(normsq)
+   const zNorm = math.zero(NormT)
+   const isRealZ = math.isZero.resolve(NormT, strategy)
+   if (V.vectorDepth === 1) {
+      const invNorm = math.invert.resolve(NormT, strategy)
+      const scalarMult = math.multiply.resolve(
+         [V, ReturnType(invNorm)], strategy)
+      return ReturnsAs(scalarMult, v => {
+         const nsq = normsq(v)
+         if (isRealZ(nsq)) return Array(v.length).fill(zNorm)
+         return scalarMult(v, invNorm(nsq))
+      })
+   }
+   // usual matrix situation, want to find a matrix whose product with v
+   // is the identity, or as close as we can get to that if the rank is
+   // deficient. We use the Moore-Penrose pseudoinverse.
+   const clone = math.clone.resolve(V, strategy)
+   const VComp = V.Component
+   const Elt = VComp.Component
+   const invElt = math.invert.resolve(Elt, strategy)
+   const det = math.determinant.resolve(V, strategy)
+   const neg = math.negate.resolve(Elt, strategy)
+   const multMM = math.multiply.resolve([V, V], strategy)
+   const multMS = math.multiply.resolve([V, ReturnType(invElt)], strategy)
+   const multVS = math.multiply.resolve([VComp, ReturnType(invElt)], strategy)
+   const multSS = math.multiply.resolve([Elt, Elt], strategy)
+   const sub = math.subtract.resolve([Elt, Elt], strategy)
+   const id = math.identity.resolve(V, strategy)
+   const abs = math.abs.resolve(Elt, strategy)
+   const gt = math.larger.resolve([NormT, NormT], strategy)
+   const isEltZ = math.isZero.resolve(Elt, strategy)
+   const zElt = math.zero(Elt)
+   const adj = math.adjoint.resolve(V, strategy)
+   const methods = {
+      invElt, det, isRealZ, neg, multMM, multMS, multVS, multSS,
+      id, abs, isEltZ, zElt, gt, sub, adj, clone
+   }
+   if (ReturnType(abs) !== NormT) {
+      throw new TypeError('type inconsistency in matrix invert')
+   }
+   return Returns(V, m => {
+      const rows = m.length
+      const cols = m[0].length
+      const nsq = normsq(m)
+      if (isRealZ(nsq)) { // all-zero matrix
+         const retval = []
+         for (let ix = 0; ix < cols; ++ix) {
+            retval.push(Array(rows).fill(zNorm))
+         }
+         return retval
+      }
+      if (rows == cols) {
+         // the inv helper will return falsy if not invertible
+         const retval = inv(m, rows, methods)
+         if (retval) return retval
+      }
+
+      return pinv(m, rows, cols, methods)
+   })
+})
+
+// Returns the inverse of a rows×rows matrix or false if not invertible
+// Note: destroys m in the inversion process.
+function inv(
+   origm,
+   rows,
+   {
+      invElt, det, isRealZ, neg, multMS, multVS, multSS,
+      id, abs, isEltZ, zElt, gt, sub, clone
+   }
+) {
+   switch (rows) {
+   case 1: return [[invElt(origm[0][0])]]
+   case 2: {
+      const dt = det(origm)
+      if (isRealZ(dt)) return false
+      const divisor = invElt(dt)
+      const [[a, b], [c, d]] = origm
+      return multMS([[d, neg(b)], [neg(c), a]], divisor)
+   }
+   default: { // Gauss-Jordan elimination
+      const m = clone(origm)
+      const B = id(m)
+      const cols = rows
+      // Loop over columns, performing row reductions:
+      for (let c = 0; c < cols; ++c) {
+         // Pivot: Find row r that has the largest entry in column c, and
+         // swap row c and r:
+         let colMax = abs(m[c][c])
+         let rMax = c
+         for (let r = c + 1; r < rows; ++r) {
+            const mag = abs(m[r][c])
+            if (gt(mag, colMax)) {
+               colMax = mag
+               rMax = r
+            }
+         }
+         if (isRealZ(colMax)) return false
+         if (rMax !== c) {
+            [m[c], m[rMax], B[c], B[rMax]]
+               = [m[rMax], m[c], B[rMax], B[c]]
+         }
+         // Normalize the cth row:
+         const normalizer = invElt(m[c][c])
+         const mc = multVS(m[c], normalizer)
+         m[c] = mc
+         const Bc = multVS(B[c], normalizer)
+         B[c] = Bc
+
+         // Eliminate nonzero values on other rows at column c
+         for (let r = 0; r < rows; ++r) {
+            if (r === c) continue
+            const mr = m[r]
+            const Br = B[r]
+            const mrc = mr[c]
+            if (!isEltZ(mr[c])) {
+               // Subtract Arc times row c from row r to eliminate A[r][c]
+               mr[c] = zElt
+               for (let s = c + 1; s < cols; ++s) {
+                  mr[s] = sub(mr[s], multSS(mrc, mc[s]))
+               }
+               for (let s = 0; s < cols; ++s) {
+                  Br[s] = sub(Br[s], multSS(mrc, Bc[s]))
+               }
+            }
+         }
+      }
+
+      return B
+   }}
+}
+
+// Calculates Moore-Penrose pseudoinverse
+// uses rank factorization per mathjs; SVD appears to be considered better
+// but not worth the effort to implement for this prototype
+function pinv(m, rows, cols, methods) {
+   const {C, F} = rankFactorization(m, rows, cols, methods)
+   const {multMM, adj} = methods
+   const Cstar = adj(C)
+   const Cpinv = multMM(inv(multMM(Cstar, C), Cstar.length, methods), Cstar)
+   const Fstar = adj(F)
+   const Fpinv = multMM(Fstar, inv(multMM(F, Fstar), F.length, methods))
+   return multMM(Fpinv, Cpinv)
+}
+
+// warning: destroys m in computing the row-reduced echelon form
+// TODO: this code should be merged with inv to the extent possible. It's
+// a very similar process.
+function rref(origm, rows, cols, methods) {
+   const {isEltZ, invElt, multVS, zElt, sub, multSS, clone} = methods
+   const m = clone(origm)
+   let lead = -1
+   for (let r = 0; r < rows && ++lead < cols; ++r) {
+      if (cols <= lead) return m
+      let i = r
+      while (isEltZ(m[i][lead])) {
+         if (++i === rows) {
+            i = r
+            if (++lead === cols) return m
+         }
+      }
+
+      if (i !== r) [m[i], m[r]] = [m[r], m[i]]
+
+      let normalizer = invElt(m[r][lead])
+      const mr = multVS(m[r], normalizer)
+      m[r] = mr
+      for (let i = 0; i < rows; ++i) {
+         if (i === r) continue
+         const mi = m[i]
+         const toRemove = mi[lead]
+         mi[lead] = zElt
+         for (let j = lead + 1; j < cols; ++j) {
+            mi[j] = sub(mi[j], multSS(toRemove, mr[j]))
+         }
+      }
+   }
+   return m
+}
+
+function rankFactorization(m, rows, cols, methods) {
+   const RREF = rref(m, rows, cols, methods)
+   const {isEltZ} = methods
+   const rankRows = RREF.map(row => row.some(elt => !isEltZ(elt)))
+   const C = m.map(row => row.filter((_, j) => j < rows && rankRows[j]))
+   const F = RREF.filter((_, i) => rankRows[i])
+   return {C, F}
+}
+
 export const multiply = [
    distributeFirst('multiply'),
    distributeSecond('multiply'),

src/vector/type.js

@@ -1,5 +1,5 @@
 import {Vector} from './Vector.js'
-import {Returns, Unknown} from '#core/Type.js'
+import {OneOf, Returns, ReturnType, Unknown} from '#core/Type.js'
 import {Any, Multiple, match} from '#core/TypePatterns.js'
 
 export const vector = match(Multiple(Any), (math, [TV]) => {
@@ -9,16 +9,125 @@ export const vector = match(Multiple(Any), (math, [TV]) => {
    return Returns(Vector(CompType), v => v)
 })
 
-export const transpose = match(Vector, (_math, V) => {
-   const wrapV = V.vectorDepth === 1
-   const Mat = wrapV ? Vector(V) : V
-   return Returns(Mat, v => {
-      if (wrapV) v = [v]
+export const determinant = match(Vector(Vector), (math, M, strategy) => {
+   const Elt = M.Component.Component
+   const cloneElt = math.clone.resolve(Elt, strategy)
+   const mult = math.multiply.resolve([Elt, Elt], strategy)
+   const sub = math.subtract.resolve(
+      [ReturnType(mult), ReturnType(mult)], strategy)
+   const isZ = math.isZero.resolve(Elt, strategy)
+   const zElt = math.zero(Elt)
+   const clone = math.clone.resolve(M, strategy)
+   const div = math.divide.resolve([ReturnType(sub), Elt], strategy)
+   const neg = math.negate.resolve(Elt, strategy)
+   return Returns(OneOf(ReturnType(clone), ReturnType(sub), Elt), origm => {
+      const rows = origm.length
+      switch (rows) {
+      case 1: return cloneElt(origm[0][0])
+      case 2: {
+         const [[a, b], [c, d]] = origm
+         return sub(mult(a, d), mult(b, c))
+      }
+      default: { // Bareiss algorithm
+         const m = clone(origm)
+         let negated = false
+         const rowIndices = [...Array(rows).keys()] // track row indices
+         // because the algorithm may swap rows
+         for (let k = 0; k < rows; ++k) {
+            let k_ = rowIndices[k]
+            if (isZ(m[k_][k])) {
+               let _k = k
+               while (++_k < rows) {
+                  if (!isZ(m[rowIndices[_k]][k])) {
+                     k_ = rowIndices[_k]
+                     rowIndices[_k] = rowIndices[k]
+                     rowIndices[k] = k_
+                     negated = !negated
+                     break
+                  }
+               }
+               if (_k === rows) return zElt
+            }
+            const piv = m[k_][k] // we now know nonzero
+            const piv_ = k === 0 ? 1 : m[rowIndices[k-1]][k-1]
+            for (let i = k + 1; i < rows; ++i) {
+               const i_ = rowIndices[i]
+               for (let j = k + 1; j < rows; ++j) {
+                  m[i_][j] = div(
+                     sub(mult(m[i_][j], piv), mult(m[i_][k], m[k_][j])),
+                     piv_)
+               }
+            }
+         }
+         const det = m[rowIndices[rows - 1]][rows - 1]
+         return negated ? neg(det) : det
+      }}
+   })
+})
+
+function identitizer(cols, zero, one) {
+   const retval = []
+   for (let ix = 0; ix < cols; ++ix) {
+      const row = Array(cols).fill(zero)
+      row[ix] = one
+      retval.push(row)
+   }
+   return retval
+}
+
+export const identity = [
+   match(Any, (math, V) => {
+      const toNum = math.number.resolve(V)
+      const zero = math.zero(V)
+      const one = math.one(V)
+      return Returns(Vector(Vector(V)), n => identitizer(toNum(n), zero, one))
+   }),
+   match(Vector, (math, V) => {
+      switch (V.vectorDepth) {
+      case 1: {
+         const Elt = V.Component
+         const one = math.one(Elt)
+         return Returns(math.typeOf(one), () => one)
+      }
+      case 2: {
+         const Elt = V.Component.Component
+         const zero = math.zero(Elt)
+         const one = math.one(Elt)
+         return Returns(V, m => identitizer(
+            m.length ? m[0].length : 0, zero, one))
+      }
+      default:
+         throw new RangeError(
+            `'identity' not implemented on ${V.vectorDepth} dimensional arrays`)
+      }
+   })
+]
+
+// transposes a 2D matrix
+function transposer(wrap, eltFun) {
+   return v => {
+      if (wrap) v = [v]
       const cols = v.length ? v[0].length : 0
       const retval = []
       for (let ix = 0; ix < cols; ++ix) {
-         retval.push(v.map(row => row[ix]))
+         retval.push(v.map(row => eltFun(row[ix])))
       }
       return retval
-   })
+   }
+}
+
+export const transpose = match(Vector, (_math, V) => {
+   const wrapV = V.vectorDepth === 1
+   const Mat = wrapV ? Vector(V) : V
+   return Returns(Mat, transposer(wrapV, elt => elt))
+})
+
+// or with conjugation:
+export const adjoint = match(Vector, (math, V, strategy) => {
+   const wrapV = V.vectorDepth === 1
+   const VComp = V.Component
+   const Elt = wrapV ? VComp : VComp.Component
+   const conj = math.conj.resolve(Elt, strategy)
+   const Mat = Vector(Vector(ReturnType(conj)))
+   return Returns(Mat, transposer(wrapV, conj))
 })