探讨JS编程语言下的高斯消元法

高斯消元法

引言

之前看到一篇博客：数列找规律的问题，这篇博客说的是如何用解五元一次方程组的方式来获取数列（长度为5）的拟合曲线。所以想到如何去解一个$n$元一次方程组，遂有此文。

基本思想

通过一系列的加减消元运算，直到得到类似$kx=b$的式子，然后逐一回代求解$x$向量（用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解）
将方程组的增广矩阵$A^{\prime }=(A|b)$实施初等行变化为行的最简行，此时该最简行作为增广矩阵对应的方程组与原方程组同解，这样通过解简化的阶梯形矩阵所对应的方程组就求出原方程组的解。

编程步骤

1. 消元,将矩阵化简为上三角形式
2. 判断解的情况
3. 回代

解的判断

• 判断一
  无解： 当消元完毕后，发现一行系数都为0，但是常数项不为0，此时无解
  多解： 当消元完毕后，发现有多行系数，常数项均为0，此时多解，有几行全为0，就有几个自由元，即变量的值可以任取，有无数种情况可以满足给出的方程组

• 判断二
  当$n$元一次方程组的系数矩阵$A$可逆时，可求出方程组解$X=A^{-1}b$，实际上这也是方程组的唯一解。如果方程组系数矩阵$A$不可逆或者$A$不是方阵时，
  齐次线性方程组的解只有两种情况：唯一解（零解）和无穷多解（非零解）

算法实现

function elimination(mat) { 
  for (let i = 1,len = mat.length; i < len; i++) {
    for (let j = 0; j < i; j++) {
      // TODO: 被除数不能为0, 每次运算时，必须保证对角线上的元素不为0(即运算中的分母不为0)
      const coe = mat[i][j] / mat[j][j] // 系数
  
      mat[i].forEach((val, idx) => {
        mat[i][idx] -= mat[j][idx] * coe
      })
    }
  }

  console.log("elimination -> mat\n", mat)
}

/**
 * 返回当前方程组的解的状态  
 * 解的状态：无解，多解，唯一解
 * @param {Matrix} mat 处理后的行阶梯阵
 */
function getSolutionStatus(mat) {
  let status = 'unique' // 解的状态：无解，多解，唯一解

  if (mat.some(row => row.slice(0, -1).every(val => val === 0) && row[row.length - 1] !== 0)) status = 'no'

  if (mat.some(row => row.every(val => val === 0))) status = 'multiple'

  return status
}
/**
 * 回代
 * @param {Matrix} mat 处理后的矩阵
 */
function back_substitution(mat) {
  const rowLen = mat.length
  const columnLen = mat[0].length
  const result = Array(rowLen).fill()

  for (let t = rowLen - 1; t >= 0; t--) {
    const divisor = result.reduceRight((pre, cur, idx) => {
      if (cur === undefined) {
        return pre
      } else {
        return pre - mat[t][idx] * cur
      }
    }, mat[t][columnLen - 1])

    result[t] = divisor / mat[t][t]
  }

  console.log("functiongaussian_elimination -> result\n", result)

  return result
}

/**
 * 高斯消元法
 * @param {Array[[]]} mat 方程组的增广矩阵，默认方程组的系数矩阵都是合法数值
 */
function gaussian_elimination(mat) {
  // mat 为增广矩阵
  // 校验参数
  const rowLen = mat.length
  const columnLen = mat[0].length

  if (columnLen - rowLen !== 1) {
    throw new Error('mat required is the augmented matrix')
  }

  // 消元
  elimination(mat)

  // 解的判断
  const solutionStatus = getSolutionStatus(mat) // 解的状态：无解，多解，唯一解


  if (solutionStatus === 'no') {
    return -1
  } else if (solutionStatus === 'multiple') {
    return []
  } else {
    // 回代
    return back_substitution(mat)
  }
}
// 测试数据
const arr = [
  [2, 1, -1, 8],
  [-3, -1, 2, -11],
  [-2, 1, 2, -3]
] // output: [2, 3, -1]

const arr2 = [
  [1, 1, 1, 1, 1, 1],
  [16, 8, 4, 2, 1, 3],
  [81, 27, 9, 3, 1, 5],
  [256, 64, 16, 4, 1, 7],
  [625, 125, 25, 5, 1, 2020]
] // output: [ 83.79166666667061,-837.9166666667011,2932.7083333334544,-4187.5833333335095,2010.0000000000857 ]

const arr3 = [
  [3, 2, 1, 6],
  [2, 2, 2, 4],
  [4, -2, -2, 2]
] // output: [ 0.9999999999999997, 2.0000000000000004, -1.0000000000000002 ]

const arr4 = [
  [47, 28, 19],
  [89, 53, 36]
] // output: [ 1, -1 ]

gaussian_elimination(arr5)

顺序消去法虽然编程操作简单，但是存在以下两方面限制:
(1). 每次运算时，必须保证对角线上的元素不为0(即运算中的分母不为0)，否则算法无法继续进行。
(2). 即使的值不为零，但如果绝对值很小，由于第$k$次运算中在分母位置，因此作除数会引起很大的误差，从而影响算法的稳定性。
为了减少计算过程中舍入误差对方程组求解的影响，因此是否可以选择绝对值尽可能大的主元作为除数，基于这种思想就有了高斯消去法的改进型:列主元消去法和全主元消去法。

列主元消去法

在第$k$步消元前，先找出$k$行下所有第$k$列元素最大的非零元素$a_{r,k}$，将第$r$行与第$k$行进行整行交换，这样既不影响原方程的解，也可以将绝对值最大的$a_{r,k}$作为主元，放在除数的位置上，尽可能减小引入误差。

全主元消去法与列主元消去法类似，只不过列主元消去法是从第$k$列中选取一个最大的元素，与第$k$行进行交换。而全主元消去法是从第k行第k列开始的右下角矩阵中所有元素中选取一个最大的元素作为主元，同时交换$r$行与$c$列，从而保证稳定性。
之前将算法分解为三个步骤：消元，判断解，回代。这样我们只需要写列主元消去法，后面的步骤可以直接调用之前写好的函数，列主元消去法代码如下：

/**
 * 列主元消去法
 * @param {Matrix} mat 方程组矩阵
 */
function columnElimination(mat) {
  const swap = (a, b) => {
    const temp = mat[a]

    mat[a] = mat[b]
    mat[b] = temp
  }
  const getMaxColumnIndex = (ci) => {
    let result = 0
    for (let i = 1, len = mat.length; i < len; i++) {
      if (Math.abs(mat[i][ci]) > Math.abs(mat[result][ci])) result = i
    }

    return result
  }

  for (let i = 1,len = mat.length; i < len; i++) {
    for (let j = 0; j < i; j++) {
      const maxCI = getMaxColumnIndex(j)

	  // 矩阵交换两行（列）不需要变号，这和行列式不同。因为：行列式是值，而矩阵是表，每一行可以看做一个记录
      if (maxCI !== j) swap(maxCI, j)

      const coe = mat[i][j] / mat[j][j] // 系数
  
      mat[i].forEach((val, idx) => {
        mat[i][idx] -= mat[j][idx] * coe
      })
    }
  }

  console.log("elimination -> mat\n", mat)
}

数列找规律要点

$f(x)$ 函数中使用的是多项式，如果你想弄得复杂点，$f(x)$中也可以使用指数函数、三角函数等等。当然，弄得太复杂的话可能会很难解出方程中常数的值来，如果你只是想凑一个公式，让数列的下一项得到任意指定的值，那么多项式就足够了。
另外，任意长度的数列找规律都可以使用这个方法，原因上面说了，这本质上就是一个给定若干散点求拟合函数的问题，使用拉格朗日插值法等方法，总是能找到符合要求的多项式函数。

参考链接

• 线性代数 —— 高斯消元法
• n元线性方程组求解
• 用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解
• 高斯消元法详解
• 数列找规律的问题