Important edit
Problem is solved. Please, take a look at my own answer in this StackOverflow question to know how.
However, here is the new (and correctly working) code :
Displayer
The displayer is the same as below.
My correct and working implementation
/**
* 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))
}
I'm trying to implement a Scala version of the Gauss-Jordan elimination to invert a matrix (NB : mutable collections and imperative paradigm are used to simplify the implementation - I tried to write the algorithm without, but it's almost impossible, consequently to the fact the algorithm contains nested steps).
My problem
The identity matrix is not well transformed into the result of the inversion. In other words : the transformation of the identity matrix into the inverted matrix (which is the result of the Gauss-Jordan's elimination) is uncorrect.
Example
Consider this matrix (A) :
(2.0, -1.0, 0.0)
(-1.0, 2.0, -1.0)
(0.0, -1.0, 2.0)
And this one (B) :
(1.0, 0.0, 0.0)
(0.0, 1.0, 0.0)
(0.0, 0.0, 1.0)
If we apply Gauss-Jordan elimination, A becomes :
(1.0, 0.0, 0.0)
(0.0, 1.0, 0.0)
(0.0, 0.0, 1.0)
If we apply Gauss-Jordan elimination, B becomes :
(0.75 0.5 0.25)
(0.5 1 0.5)
(0.25 0.5 0.75)
If we apply my implementation, there is no problem for A, since I obtain the following matrix :
(1.0, 0.0, 0.0)
(0.0, 1.0, 0.0)
(0.0, 0.0, 1.0)
If we apply my implementation, however, B is not well transformed, since I obtain the following matrix :
(1.0, 0.5, 0.0)
(1.0, 0.5, 0.6666666666666666)
(0.0, 1.0, 0.33333333333333337)
Gauss-Jordan's elimination : explanations about this algorithm
It proceeds in 3 steps, column by column. Those steps are :
- We find the max^1 coefficient in the current column^2. If it equals 0, it means the matrix can't be inverted and the algorithm returns this error. Otherwise, we swap the line containing the max coefficient with the line containing the pivot : in other words, we change the pivot with the max coefficient of the column (NB : the whole lines are swapped). ^1 : max is the used function only for divisions accuracy reasons (divisions done in STEP 2). Another function would be the random one.
^2 : the max coefficient in the current column is found from the (z+1)th line, where z is the ID of the last pivot we used (i.e. : the ID of the last worked column)
We divide the whole line containing the pivot we got at STEP 1 by the pivot, to set the pivot to 1 (in later sentences, the expression "the pivot" systematicaly refers to this pivot we got at STEP 1). By the way, note the less important fact that the other coefficients of this same line are also divided (cf. "we divide the whole line").
We subtract each whole line of the current column by itself times the pivot's line, to set all the current column's coefficient to 0. By the way, note the less important fact that the other coefficients of these same lines are also subtracted (cf. "we subtract each whole line").
STEP 3 and STEP 2 are realized within STEP 1 (i.e. : those are nested STEPS). STEP 3 must be realized after STEP 2 (to make use of the pivot's value = 1 in the {subtractions and multiplications} realized in STEP 3.
Gauss-Jordan's elimination : my unworking implementation
Input
val m : Matrix = new Matrix(Seq(Seq(2, -1, 0), Seq(-1, 2, -1), Seq(0, -1, 2)))
Unworking implementation of this algorithm
Displayer
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)
My unworking implementation
/**
* 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))
}
Execution and Output
You can find the execution of my implementation with this input here : https://jsfiddle.net/wwhdu32x/
Troubleshooting
You can find a troubleshooting here : https://jsfiddle.net/wwhdu32x/1/ (comments beginning by "ERROR" are written - NB : this troubleshooting only concerns the first iteration, i.e. the first column).
My question
Why my identity matrix isn't well transformed ? How could I deal with it ?
