From db86b2ecd02ca28f94c3e9569891813bad46bf26 Mon Sep 17 00:00:00 2001 From: Glen Whitney Date: Tue, 6 May 2025 16:53:11 -0700 Subject: [PATCH] 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 | 22 ++- src/vector/__test__/arithmetic.spec.js | 28 ++++ src/vector/__test__/type.spec.js | 20 +++ src/vector/arithmetic.js | 197 ++++++++++++++++++++++++- src/vector/type.js | 125 +++++++++++++++- 5 files changed, 379 insertions(+), 13 deletions(-) diff --git a/src/number/type.js b/src/number/type.js index d8e6806..0d48348 100644 --- a/src/number/type.js +++ b/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)) + }) + }) ] diff --git a/src/vector/__test__/arithmetic.spec.js b/src/vector/__test__/arithmetic.spec.js index 5dcc3ca..ce888d0 100644 --- a/src/vector/__test__/arithmetic.spec.js +++ b/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]) diff --git a/src/vector/__test__/type.spec.js b/src/vector/__test__/type.spec.js index 369033a..1208ffc 100644 --- a/src/vector/__test__/type.spec.js +++ b/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) + }) }) diff --git a/src/vector/arithmetic.js b/src/vector/arithmetic.js index 1687138..71332c2 100644 --- a/src/vector/arithmetic.js +++ b/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'), diff --git a/src/vector/type.js b/src/vector/type.js index 88cbe3f..ad549d5 100644 --- a/src/vector/type.js +++ b/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)) })