Gauss-Jordan matrix inversion

From Rosetta Code

Gauss-Jordan matrix inversion

Gauss-Jordan matrix inversion is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Task

Invert matrix   A   using Gauss-Jordan method.

A   being an   n by n   matrix.

C#[edit]

 
using System;
 
namespace Rosetta
{
internal class Vector
{
private double[] b;
internal readonly int rows;
 
internal Vector(int rows)
{
this.rows = rows;
b = new double[rows];
}
 
internal Vector(double[] initArray)
{
b = (double[])initArray.Clone();
rows = b.Length;
}
 
internal Vector Clone()
{
Vector v = new Vector(b);
return v;
}
 
internal double this[int row]
{
get { return b[row]; }
set { b[row] = value; }
}
 
internal void SwapRows(int r1, int r2)
{
if (r1 == r2) return;
double tmp = b[r1];
b[r1] = b[r2];
b[r2] = tmp;
}
 
internal double norm(double[] weights)
{
double sum = 0;
for (int i = 0; i < rows; i++)
{
double d = b[i] * weights[i];
sum += d*d;
}
return Math.Sqrt(sum);
}
 
internal void print()
{
for (int i = 0; i < rows; i++)
Console.WriteLine(b[i]);
Console.WriteLine();
}
 
public static Vector operator-(Vector lhs, Vector rhs)
{
Vector v = new Vector(lhs.rows);
for (int i = 0; i < lhs.rows; i++)
v[i] = lhs[i] - rhs[i];
return v;
}
}
 
class Matrix
{
private double[] b;
internal readonly int rows, cols;
 
internal Matrix(int rows, int cols)
{
this.rows = rows;
this.cols = cols;
b = new double[rows * cols];
}
 
internal Matrix(int size)
{
this.rows = size;
this.cols = size;
b = new double[rows * cols];
for (int i = 0; i < size; i++)
this[i, i] = 1;
}
 
internal Matrix(int rows, int cols, double[] initArray)
{
this.rows = rows;
this.cols = cols;
b = (double[])initArray.Clone();
if (b.Length != rows * cols) throw new Exception("bad init array");
}
 
internal double this[int row, int col]
{
get { return b[row * cols + col]; }
set { b[row * cols + col] = value; }
}
 
public static Vector operator*(Matrix lhs, Vector rhs)
{
if (lhs.cols != rhs.rows) throw new Exception("I can't multiply matrix by vector");
Vector v = new Vector(lhs.rows);
for (int i = 0; i < lhs.rows; i++)
{
double sum = 0;
for (int j = 0; j < rhs.rows; j++)
sum += lhs[i,j]*rhs[j];
v[i] = sum;
}
return v;
}
 
internal void SwapRows(int r1, int r2)
{
if (r1 == r2) return;
int firstR1 = r1 * cols;
int firstR2 = r2 * cols;
for (int i = 0; i < cols; i++)
{
double tmp = b[firstR1 + i];
b[firstR1 + i] = b[firstR2 + i];
b[firstR2 + i] = tmp;
}
}
 
//with partial pivot
internal bool InvPartial()
{
const double Eps = 1e-12;
if (rows != cols) throw new Exception("rows != cols for Inv");
Matrix M = new Matrix(rows); //unitary
for (int diag = 0; diag < rows; diag++)
{
int max_row = diag;
double max_val = Math.Abs(this[diag, diag]);
double d;
for (int row = diag + 1; row < rows; row++)
if ((d = Math.Abs(this[row, diag])) > max_val)
{
max_row = row;
max_val = d;
}
if (max_val <= Eps) return false;
SwapRows(diag, max_row);
M.SwapRows(diag, max_row);
double invd = 1 / this[diag, diag];
for (int col = diag; col < cols; col++)
{
this[diag, col] *= invd;
}
for (int col = 0; col < cols; col++)
{
M[diag, col] *= invd;
}
for (int row = 0; row < rows; row++)
{
d = this[row, diag];
if (row != diag)
{
for (int col = diag; col < this.cols; col++)
{
this[row, col] -= d * this[diag, col];
}
for (int col = 0; col < this.cols; col++)
{
M[row, col] -= d * M[diag, col];
}
}
}
}
b = M.b;
return true;
}
 
internal void print()
{
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
Console.Write(this[i,j].ToString()+" ");
Console.WriteLine();
}
}
}
}
 
 
using System;
 
namespace Rosetta
{
class Program
{
static void Main(string[] args)
{
Matrix M = new Matrix(4, 4, new double[] { -1, -2, 3, 2, -4, -1, 6, 2, 7, -8, 9, 1, 1, -2, 1, 3 });
M.InvPartial();
M.print();
}
}
}
 
Output:

-0.913043478260869 0.246376811594203 0.0942028985507246 0.413043478260869 -1.65217391304348 0.652173913043478 0.0434782608695652 0.652173913043478 -0.695652173913043 0.36231884057971 0.0797101449275362 0.195652173913043 -0.565217391304348 0.231884057971014 -0.0289855072463768 0.565217391304348

J[edit]

Solution:

Uses Gauss-Jordan implementation (as described in Reduced_row_echelon_form#J) to find reduced row echelon form of the matrix after augmenting with an identity matrix.

require 'math/misc/linear'
augmentR_I1=: ,. [email protected]@# NB. augment matrix on the right with its Identity matrix
matrix_invGJ=: # }."1 [: [email protected]_I1

Usage:

   ]A =: 1 2 3, 4 1 6,: 7 8 9
1 2 3
4 1 6
7 8 9
matrix_invGJ A
_0.8125 0.125 0.1875
0.125 _0.25 0.125
0.520833 0.125 _0.145833

Julia[edit]

Works with: Julia version 0.6

Built-in LAPACK-based linear solver uses partial-pivoted Gauss elimination):

A = [1 2 3; 4 1 6; 7 8 9]
@show I / A
@show inv(A)

Native implementation:

function gaussjordan(A::Matrix)
size(A, 1) == size(A, 2) || throw(ArgumentError("A must be squared"))
n = size(A, 1)
M = [convert(Matrix{float(eltype(A))}, A) I]
i = 1
local tmp = Vector{eltype(M)}(2n)
# forward
while i ≤ n
if M[i, i] ≈ 0.0
local j = i + 1
while j ≤ n && M[j, i] ≈ 0.0
j += 1
end
if j ≤ n
tmp .= M[i, :]
M[i, :] .= M[j, :]
M[j, :] .= tmp
else
throw(ArgumentError("matrix is singular, cannot compute the inverse"))
end
end
for j in (i + 1):n
M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
end
i += 1
end
i = n
# backward
while i ≥ 1
if M[i, i] ≈ 0.0
local j = i - 1
while j ≥ 1 && M[j, i] ≈ 0.0
j -= 1
end
if j ≥ 1
tmp .= M[i, :]
M[i, :] .= M[j, :]
M[j, :] .= tmp
else
throw(ArgumentError("matrix is singular, cannot compute the inverse"))
end
end
for j in (i - 1):-1:1
M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
end
i -= 1
end
M ./= diag(M) # normalize
return M[:, n+1:2n]
end
 
@show gaussjordan(A)
@assert gaussjordan(A) ≈ inv(A)
 
A = rand(10, 10)
@assert gaussjordan(A) ≈ inv(A)
Output:
I / A = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]
inv(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]
gaussjordan(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]

Kotlin[edit]

This follows the description of Gauss-Jordan elimination in Wikipedia whereby the original square matrix is first augmented to the right by its identity matrix, its reduced row echelon form is then found and finally the identity matrix to the left is removed to leave the inverse of the original square matrix.

// version 1.2.21
 
typealias Matrix = Array<DoubleArray>
 
fun Matrix.inverse(): Matrix {
val len = this.size
require(this.all { it.size == len }) { "Not a square matrix" }
val aug = Array(len) { DoubleArray(2 * len) }
for (i in 0 until len) {
for (j in 0 until len) aug[i][j] = this[i][j]
// augment by identity matrix to right
aug[i][i + len] = 1.0
}
aug.toReducedRowEchelonForm()
val inv = Array(len) { DoubleArray(len) }
// remove identity matrix to left
for (i in 0 until len) {
for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
}
return inv
}
 
fun Matrix.toReducedRowEchelonForm() {
var lead = 0
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
if (colCount <= lead) return
var i = r
 
while (this[i][lead] == 0.0) {
i++
if (rowCount == i) {
i = r
lead++
if (colCount == lead) return
}
}
 
val temp = this[i]
this[i] = this[r]
this[r] = temp
 
if (this[r][lead] != 0.0) {
val div = this[r][lead]
for (j in 0 until colCount) this[r][j] /= div
}
 
for (k in 0 until rowCount) {
if (k != r) {
val mult = this[k][lead]
for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
}
}
 
lead++
}
}
 
fun Matrix.printf(title: String) {
println(title)
val rowCount = this.size
val colCount = this[0].size
 
for (r in 0 until rowCount) {
for (c in 0 until colCount) {
if (this[r][c] == -0.0) this[r][c] = 0.0 // get rid of negative zeros
print("${"% 10.6f".format(this[r][c])} ")
}
println()
}
 
println()
}
 
fun main(args: Array<String>) {
val a = arrayOf(
doubleArrayOf(1.0, 2.0, 3.0),
doubleArrayOf(4.0, 1.0, 6.0),
doubleArrayOf(7.0, 8.0, 9.0)
)
a.inverse().printf("Inverse of A is :\n")
 
val b = arrayOf(
doubleArrayOf( 2.0, -1.0, 0.0),
doubleArrayOf(-1.0, 2.0, -1.0),
doubleArrayOf( 0.0, -1.0, 2.0)
)
b.inverse().printf("Inverse of B is :\n")
}
Output:
Inverse of A is :

 -0.812500    0.125000    0.187500  
  0.125000   -0.250000    0.125000  
  0.520833    0.125000   -0.145833  

Inverse of B is :

  0.750000    0.500000    0.250000  
  0.500000    1.000000    0.500000  
  0.250000    0.500000    0.750000  

Perl 6[edit]

Works with: Rakudo version 2018.03

Uses bits and pieces from other tasks, Reduced row echelon form primarily.

sub gauss-jordan-invert (@m where *.&is-square) {
^@m .map: { @m[$_].append: identity(+@m)[$_] };
@m.&rref[*]»[+@m .. *];
}
 
sub is-square (@m) { so @m == all @m[*] }
 
sub identity ($n) { [ 1, |(0 xx $n-1) ], *.rotate(-1) ... *.tail }
 
# reduced row echelon form (Gauss-Jordan elimination)
sub rref (@m) {
return unless @m;
my ($lead, $rows, $cols) = 0, +@m, +@m[0];
 
for ^$rows -> $r {
$lead < $cols or return @m;
my $i = $r;
until @m[$i;$lead] {
++$i == $rows or next;
$i = $r;
++$lead == $cols and return @m;
}
@m[$i, $r] = @m[$r, $i] if $r != $i;
my $lv = @m[$r;$lead];
@m[$r] »/=» $lv;
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »*» (@m[$n;$lead] // 0);
}
++$lead;
}
@m
}
 
sub rat-or-int ($num) {
return $num unless $num ~~ Rat;
return $num.narrow if $num.narrow.WHAT ~~ Int;
$num.nude.join: '/';
}
 
sub say_it ($message, @array) {
say "\n$message";
$_».&rat-or-int.fmt(" %6s").put for @array;
}
 
my @tests = (
[
[ 1, 2, 3 ],
[ 4, 1, 6 ],
[ 7, 8, 9 ]
],
[
[ 2, -1, 0 ],
[-1, 2, -1 ],
[ 0, -1, 2 ]
],
[
[ -1, -2, 3, 2 ],
[ -4, -1, 6, 2 ],
[ 7, -8, 9, 1 ],
[ 1, -2, 1, 3 ]
],
);
 
 
for @tests -> @matrix {
say_it( 'Original Matrix:', @matrix );
say_it( 'Gauss-Jordan Inverted Matrix:', gauss-jordan-invert @matrix );
say "\n";
}
 
Output:
Original Matrix:
      1       2       3
      4       1       6
      7       8       9

Gauss-Jordan Inverted Matrix:
 -13/16     1/8    3/16
    1/8    -1/4     1/8
  25/48     1/8   -7/48



Original Matrix:
      2      -1       0
     -1       2      -1
      0      -1       2

Gauss-Jordan Inverted Matrix:
    3/4     1/2     1/4
    1/2       1     1/2
    1/4     1/2     3/4



Original Matrix:
     -1      -2       3       2
     -4      -1       6       2
      7      -8       9       1
      1      -2       1       3

Gauss-Jordan Inverted Matrix:
 -21/23   17/69  13/138   19/46
 -38/23   15/23    1/23   15/23
 -16/23   25/69  11/138    9/46
 -13/23   16/69   -2/69   13/23


Sidef[edit]

Uses the rref(M) function from Reduced row echelon form.

Translation of: Perl 6
func gauss_jordan_invert (M) {
 
var I = M.len.of {|i|
M.len.of {|j|
i == j ? 1 : 0
}
}
 
var A = gather {
^M -> each {|i| take(M[i] + I[i]) }
}
 
rref(A).map { .last(M.len) }
}
 
var A = [
[-1, -2, 3, 2],
[-4, -1, 6, 2],
[ 7, -8, 9, 1],
[ 1, -2, 1, 3],
]
 
say gauss_jordan_invert(A).map {
.map { "%6s" % .as_rat }.join(" ")
}.join("\n")
Output:
-21/23   17/69  13/138   19/46
-38/23   15/23    1/23   15/23
-16/23   25/69  11/138    9/46
-13/23   16/69   -2/69   13/23

zkl[edit]

This uses GSL to invert a matrix via LU decomposition, not Gauss-Jordan.

var [const] GSL=Import.lib("zklGSL");    // libGSL (GNU Scientific Library)
m:=GSL.Matrix(3,3).set(1,2,3, 4,1,6, 7,8,9);
i:=m.invert();
i.format(10,4).println("\n");
(m*i).format(10,4).println();
Output:
   -0.8125,    0.1250,    0.1875
    0.1250,   -0.2500,    0.1250
    0.5208,    0.1250,   -0.1458

    1.0000,    0.0000,    0.0000
   -0.0000,    1.0000,    0.0000
   -0.0000,    0.0000,    1.0000
m:=GSL.Matrix(3,3).set(2,-1,0, -1,2,-1, 0,-1,2);
m.invert().format(10,4).println("\n");
Output:
    0.7500,    0.5000,    0.2500
    0.5000,    1.0000,    0.5000
    0.2500,    0.5000,    0.7500