algorithm - 我的 Gauss-Jordan 消元法有什么问题？

Question

重要修改

问题解决了。请看一下我自己在这个 StackOverflow 问题中的回答，以了解如何做。

但是，这是新的(并且可以正常工作的)代码:

显示

显示器同下。

我正确且有效的实现

/**
  * Returns the identity matrix of the specified dimension
  * @param size the number of columns (i.e. the number of rows) of the desired identity matrix
  * @return the identity matrix of the specified dimension
  */
def getIdentityMatrix(size : Int): scala.collection.mutable.Seq[scala.collection.mutable.Seq[Double]] = {
  scala.collection.mutable.Seq.tabulate(size)(r => scala.collection.mutable.Seq.tabulate(size)(c => if(r == c) 1.0 else 0.0))
}

  /**
    * This algorithm processes column by column.
    * STEP 1. It finds the greatest coefficient for the current column (called 'a') and, if it equals 0, returns NULL (since the matrix
    * can't be inverted) ; otherwise (STEP 2.), it swaps the pivot's line with this new line and the pivot becomes the adequate coefficient
    * of this new line
    * STEP 3. It divides the pivot's line by the pivot
    * STEP 4. It sets each of the current column's coefficient to 0 by subtracting the corresponding lines by the pivot's line
    * @return
    */
  def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {

    // We get first the matrix to be inverted, second the identity one
    val mutable_being_inversed_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = scala.collection.mutable.Seq(content.map(ms => scala.collection.mutable.Seq(ms:_*)):_*)
    val identity_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = getIdentityMatrix(content.length)  // We get the identity matrix. It will be modified
                                                                      // as the original matrix will.

    var id_last_pivot : Int = 0  // ID of the last pivot, i.e. ID of the current column
    content.indices.foreach(general_id_column => {
      println("Current column : " + general_id_column)

      //  STEP 1.
      val id_line_with_max_coefficient_in_this_column = (id_last_pivot until content.length).maxBy(id_line_in_this_column => Math.abs(mutable_being_inversed_matrix(id_line_in_this_column)(general_id_column)))

      if(mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)(general_id_column) == 0) {
        println("The Gauss-Jordan elimination's algorithm returns an error : indeed, the matrix can't be inverted")

      } else {

        //  STEP 2.
        val tmp_line : scala.collection.mutable.Seq[Double] = mutable_being_inversed_matrix(id_last_pivot)
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)
        mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column) = tmp_line

        val identity_tmp_line : scala.collection.mutable.Seq[Double] = identity_matrix(id_last_pivot)
        identity_matrix(id_last_pivot) = identity_matrix(id_line_with_max_coefficient_in_this_column)
        identity_matrix(id_line_with_max_coefficient_in_this_column) = identity_tmp_line
        println("\nSWAP DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

        //  STEP 3.
        val tmp = mutable_being_inversed_matrix(id_last_pivot)(general_id_column)
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / tmp)
        identity_matrix(id_last_pivot) = identity_matrix(id_last_pivot).map(coefficient => coefficient / tmp)

        println("\nDIVISION DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

        //  STEP 4.
        content.indices.foreach(id_line => {
          val tmp = mutable_being_inversed_matrix(id_line)(general_id_column)

          if(id_line != id_last_pivot) {
            content.indices.foreach(id_column => {
              mutable_being_inversed_matrix(id_line)(id_column) -= mutable_being_inversed_matrix(id_last_pivot)(id_column) * tmp
              identity_matrix(id_line)(id_column) -= identity_matrix(id_last_pivot)(id_column) * tmp
            })
          }

        })

        println("\nSUBTRACTION & MULTIPLICATION DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))
        println()

        id_last_pivot += 1

      }

    })

    (new Matrix(identity_matrix), new Matrix(mutable_being_inversed_matrix))
  }

---

我正在尝试实现 Gauss-Jordan 消去法的 Scala 版本来反转矩阵(注意:可变集合和命令式范例用于简化实现 - 我试图在没有的情况下编写算法，但这几乎是不可能的，因此，该算法包含嵌套步骤)。

我的问题

单位矩阵没有很好地转化为求逆的结果。换句话说:将单位矩阵变换为倒矩阵(这是Gauss-Jordan消去法的结果)是不正确的。

例子

考虑这个矩阵 (A) :

(2.0, -1.0, 0.0)

(-1.0, 2.0, -1.0)

(0.0, -1.0, 2.0)

还有这个(B):

(1.0, 0.0, 0.0)

(0.0, 1.0, 0.0)

(0.0, 0.0, 1.0)

如果我们应用 Gauss-Jordan 消去法，A 变为:

(1.0, 0.0, 0.0)

(0.0, 1.0, 0.0)

(0.0, 0.0, 1.0)

如果我们应用 Gauss-Jordan 消去法，B 变为:

(0.75 0.5 0.25)

(0.5 1 0.5)

(0.25 0.5 0.75)

如果我们应用我的实现，A 就没有问题，因为我获得了以下矩阵:

(1.0, 0.0, 0.0)

(0.0, 1.0, 0.0)

(0.0, 0.0, 1.0)

但是，如果我们应用我的实现，B 并没有很好地转换，因为我获得了以下矩阵:

(1.0, 0.5, 0.0)

(1.0, 0.5, 0.6666666666666666)

(0.0, 1.0, 0.33333333333333337)

Gauss-Jordan消元法:算法解释

它分 3 个步骤逐列进行。这些步骤是:

1. 我们在当前列 ^2 中找到最大 ^1 系数。如果它等于 0，则表示矩阵无法反转，算法会返回此错误。否则，我们将包含最大系数的行与包含主元的行交换:换句话说，我们用列的最大系数更改主元(注意:交换整行)。 ^1 : max 是仅出于除法精度原因而使用的函数(在步骤 2 中完成的除法)。另一个函数是随机函数。

^2 : 当前列中的最大系数是从第 (z+1) 行找到的，其中 z 是我们使用的最后一个基准的 ID(即:最后一个工作列的 ID)

1. 我们将包含我们在第 1 步获得的枢轴的整条线除以枢轴，将枢轴设置为 1(在后面的句子中，表达“枢轴”系统地指代我们在第 1 步获得的这个枢轴).顺便说一下，注意一个不太重要的事实，即同一行的其他系数也被划分(参见“我们划分整条线”)。

2. 我们减去当前列的每一整行乘以主元的行，将所有当前列的系数设置为 0。顺便说一下，请注意一个不太重要的事实，即这些相同行的其他系数是也减去(参见“我们减去每一整行”)。

STEP 3 和 STEP 2 在 STEP 1 中实现(即:它们是嵌套的 STEPS)。 STEP 3 必须在 STEP 2 之后实现(以便在 STEP 3 实现的{减法和乘法}中使用 pivot 的值 = 1。

Gauss-Jordan 消元法:我的无效实现

输入

val m : Matrix = new Matrix(Seq(Seq(2, -1, 0), Seq(-1, 2, -1), Seq(0, -1, 2)))

该算法的无效实现

显示

val m : Matrix = new Matrix(Seq(Seq(2, -1, 0), Seq(-1, 2, -1), Seq(0, -1, 2)))
println("ORIGINAL MATRIX =\n" + m)
println
val result : (Matrix, Matrix) = m.getGaussJordanInvertedMatrix
println()
println("RESULT =\n" + Console.BLUE + "Original matrix :\n" + Console.RESET + result._2 + Console.RED + "\nIdentity matrix :\n" + Console.RESET + result._1)

我的无效实现

/**
  * Returns the identity matrix of the specified dimension
  * @param size the number of columns (i.e. the number of rows) of the desired identity matrix
  * @return the identity matrix of the specified dimension
  */
def getIdentityMatrix(size : Int): scala.collection.mutable.Seq[scala.collection.mutable.Seq[Double]] = {
  scala.collection.mutable.Seq.tabulate(size)(r => scala.collection.mutable.Seq.tabulate(size)(c => if(r == c) 1.0 else 0.0))
}

  /**
    * This algorithm processes column by column.
    * STEP 1. It finds the greatest coefficient for the current column (called 'a') and, if it equals 0, returns NULL (since the matrix
    * can't be inverted) ; otherwise (STEP 2.), it swaps the pivot's line with this new line and the pivot becomes the adequate coefficient
    * of this new line
    * STEP 3. It divides the pivot's line by the pivot
    * STEP 4. It sets each of the current column's coefficient to 0 by subtracting the corresponding lines by the pivot's line
    * @return
    */
  def getGaussJordanInvertedMatrix: (Matrix, Matrix) = {

    // We get first the matrix to be inverted, second the identity one
    val mutable_being_inversed_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = scala.collection.mutable.Seq(content.map(ms => scala.collection.mutable.Seq(ms:_*)):_*)
    val identity_matrix : collection.mutable.Seq[collection.mutable.Seq[Double]] = getIdentityMatrix(content.length)  // We get the identity matrix. It will be modified
                                                                      // as the original matrix will.

    var id_last_pivot : Int = 0  // ID of the last pivot, i.e. ID of the current column
    content.indices.foreach(general_id_column => {
      println("Current column : " + general_id_column)

      //  STEP 1.
      val id_line_with_max_coefficient_in_this_column = (id_last_pivot until content.length).maxBy(id_line_in_this_column => Math.abs(mutable_being_inversed_matrix(id_line_in_this_column)(general_id_column)))

      if(mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)(general_id_column) == 0) {
        println("The Gauss-Jordan elimination's algorithm returns an error : indeed, the matrix can't be inverted")

      } else {

        //  STEP 2.
        val tmp_line : scala.collection.mutable.Seq[Double] = mutable_being_inversed_matrix(id_last_pivot)
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column)
        mutable_being_inversed_matrix(id_line_with_max_coefficient_in_this_column) = tmp_line

        val identity_tmp_line : scala.collection.mutable.Seq[Double] = identity_matrix(id_last_pivot)
        identity_matrix(id_last_pivot) = identity_matrix(id_line_with_max_coefficient_in_this_column)
        identity_matrix(id_line_with_max_coefficient_in_this_column) = identity_tmp_line
        println("\nSWAP DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

        //  STEP 3.
        mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / mutable_being_inversed_matrix(id_last_pivot)(general_id_column))
        identity_matrix(id_last_pivot) = identity_matrix(id_last_pivot).map(coefficient => coefficient / mutable_being_inversed_matrix(id_last_pivot)(general_id_column))

        println("\nDIVISION DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))

        //  STEP 4.
        content.indices.foreach(id_line => {
          val tmp = mutable_being_inversed_matrix(id_line)(general_id_column)

          if(id_line != id_last_pivot) {
            content.indices.foreach(id_column => {
              mutable_being_inversed_matrix(id_line)(id_column) -= mutable_being_inversed_matrix(id_last_pivot)(id_column) * tmp
              identity_matrix(id_line)(id_column) -= mutable_being_inversed_matrix(id_last_pivot)(id_column) * tmp
            })
          }

        })

        println("\nSUBTRACTION & MULTIPLICATION DONE")
        println(Console.BLUE + "Original matrix :\n" + Console.RESET + mutable_being_inversed_matrix.mkString("\n"))
        println(Console.RED + "Identity matrix :\n" + Console.RESET + identity_matrix.mkString("\n"))
        println()

        id_last_pivot += 1

      }

    })

    (new Matrix(identity_matrix), new Matrix(mutable_being_inversed_matrix))
  }

执行与输出

您可以在此处找到使用此输入执行我的实现:https://jsfiddle.net/wwhdu32x/

疑难解答

您可以在此处找到故障排除:https://jsfiddle.net/wwhdu32x/1/ (以“ERROR”开头的注释被写入 - 注意:此故障排除仅涉及第一次迭代，即第一列)。

我的问题

为什么我的单位矩阵没有很好地转换？我该如何处理？

Answer 1

问题解决了。该问题已更新，其中包括新代码(旧代码仍然可用，以便进行比较)。有两个错误(下面的“STEP XYZ”引用了相应源代码的 STEP，而不是这个 StackOverflow 问题中提到的步骤，呈现方式略有不同):

1. 关于单位矩阵的减法没有使用单位矩阵的系数(第 4 步)。错误修复:identity_matrix(id_line)(id_column) -= identity_matrix(id_last_pivot)(id_column) * tmp

2. 其次，在第 3 步中，我忘记将主元存储在一个临时变量中，以便用它除以两个矩阵(原始矩阵和恒等矩阵)。在不存储的情况下，主元的值在对原始矩阵进行除法后发生了变化。错误修复:

       val tmp = mutable_being_inversed_matrix(id_last_pivot)(general_id_column)
       mutable_being_inversed_matrix(id_last_pivot) = mutable_being_inversed_matrix(id_last_pivot).map(coefficient => coefficient / tmp)
       identity_matrix(id_last_pivot) = identity_matrix(id_last_pivot).map(coefficient => coefficient / tmp)