图形拟合（一）：矩阵运算函数库

背景

基于矩阵运算实现几何图形拟合功能，但调研过JavaScript版本的数学运算库，几乎都不支持矩阵奇异值分解（SVD），故采用math-js重新封装了一个矩阵运算库。

代码

1. 按需引入math-js的基础函数

import {
  create,
  multiplyDependencies,
  sqrtDependencies,
  forEachDependencies,
  matrixDependencies,
  diagDependencies,
  invDependencies,
  transposeDependencies,
  eigsDependencies,
  sizeDependencies,
  matrixFromColumnsDependencies,
  columnDependencies,
  reshapeDependencies,
  countDependencies
} from 'mathjs';

const {
  multiply,
  sqrt,
  forEach,
  matrix,
  diag,
  inv,
  transpose,
  eigs,
  size,
  matrixFromColumns,
  column,
  reshape,
  count,
} = create({
  multiplyDependencies,
  sqrtDependencies,
  forEachDependencies,
  matrixDependencies,
  diagDependencies,
  invDependencies,
  transposeDependencies,
  eigsDependencies,
  sizeDependencies,
  matrixFromColumnsDependencies,
  columnDependencies,
  reshapeDependencies,
  countDependencies
}, {})

2. 封装常用矩阵运算方法

class Mat {
  constructor() {}

  // 初始化matrix 
  static matrix(row: number, col: number, _data: number[]) {
    const _d = _data.slice(0, row * col)
    return reshape(_d, [row, col]);
  }

  // 对角矩阵
  static diag(mat: any) {
    return diag(mat)
  }

  // 矩阵求逆
  static inv(mat: any) {
    return inv(mat);
  }

  // 矩阵转置
  static transpose(mat: any) {
    return transpose(mat);
  }

  // 矩阵乘法
  static multiply(...mats: any) {
    return multiply(...mats);
  }

  // 矩阵遍历
  static forEach(mat: any, cb: Function) {
    forEach(mat, cb);
  }

  // 矩阵SVD分解
  static svd(mat: any) {
    let A = matrix(mat);
    let s : any = size(mat)
    let c = count(mat)
    let A_T = transpose(A);
    let A_TA = multiply(A_T, A);
    let eigsA_TA : any = eigs(A_TA, 1e-8);
    let U = null
    let V = null
    let E = null
    let valA_TA : any = eigsA_TA.values
    let vecA_TA : any = eigsA_TA.vectors
    let lenValA_TA : any = size(valA_TA)
    lenValA_TA = lenValA_TA._data[0]
    let temp_e = new Array(c).fill(0);
    for (let i = lenValA_TA - 1; i >= 0; --i) {
      let r = lenValA_TA - 1 - i;
      temp_e[r * s[1] + r] = sqrt(valA_TA._data[i]);
    }
    E = reshape(matrix(temp_e), s)

    let temp_v = []
    let temp_u = []
    let diagE = diag(E)

    for (let i = lenValA_TA - 1; i >= 0; i-- ) {
      let col = column(vecA_TA, i)
      temp_v.push(col)
      temp_u.push(multiply(1 / diagE._data[lenValA_TA - 1 - i], A, col))
    }

    V = matrixFromColumns(...temp_v)
    U = matrixFromColumns(...temp_u)

    return {
      U,
      V,
      E
    }
  }

  // 矩阵反向回代
  static backSubst(W: any, U: any, Vt: any, b: any) {
    let Vtt = transpose(Vt)
    let Ut = transpose(U)
    let sV = size(Vtt)
    let sU = size(Ut)
    let invDiagW = inv(diag(diag(W)))
    let newDiag = reshape(invDiagW, [sV._data[1], sU._data[0]])
    return multiply(Vtt, newDiag, Ut, b)
  }

  // 矩阵求解
  static solve(A: any, b: any) {
    let {
      E,
      U,
      V
    } = Mat.svd(A)
    V = transpose(V);

    return Mat.backSubst(E, U, V, b);
  }
}