Matrix Inversion

[GWS]

Hi folks,
You never know when a nice inverted matrix might come in useful .. 

Here's a method which takes an original NxN matrix, and calculates the inverse of the array without
using any other storage - the result is produced within the array space itself.

Of course, the original array values are lost, but in most applications it is only the inverse that is required.

Matrix inversion is a bit of a sausage machine, and it's best to check the process has not lost precision
during the calculation. 

To do this, keep a copy of the original array and after the inversion, multiply the original and the inverse arrays.
The result should be an Identity array - that is, where the diagonal values = 1.0 and all other values = 0.0

For example:

1 0 0
0 1 0
0 0 1

Here's the test program code for subroutine Invert( ) :

Code Select Expand



openconsole
cls
autodefine "OFF"
setprecision 5

def N:int
' set an example array to hold 3x3 elements (indices starting at 0)..
N = 3

def A[N,N]:double

' load the array in column order ..
' (this is necessary because Creative stores array values in column major order).
' (when initialising, the second number refers to the column).
A[0,0] = 2,4,1
A[0,1] = 1,3,1
A[0,2] = 0,1,1

' so the actual array A[ ] is ..
' 2 1 0
' 4 3 1
' 1 1 1

' print the original array values ..
for i = 0 to N-1
for j = 0 to N-1
print A[i,j],
next j
print
next i

invert(A[])

' print the inverse values ..
print
print "Inverse"
print
for i = 0 to N-1
for j = 0 to N-1
print A[i,j],
next j
print
next i

do:until inkey$<>""
closeconsole
end

sub invert(N,A[]:double)
' routine to invert an NxN array A[ ] inplace (ie. the inverse is formed using the array A[ ] space).
def i,j,k,ik,i1:int
def kk,jj,ii:int
def d:double

for k = 1 to N
kk = k - 1
d = A[kk,kk]
A[kk,kk] = 1.0
for j = 1 to N
jj = j - 1
A[kk,jj] = A[kk,jj] / d
next j
if (k = N) then goto nextpart
ik = k + 1
for i = ik to N
ii = i - 1
d = A[ii,kk]
A[ii,kk] = 0.0
for j = 1 to N
jj = j - 1
A[ii,jj] = A[ii,jj] - D * A[kk,jj]
next j
next i
next k

nextpart:
ik = N - 1
for k = 1 to ik
kk = k - 1
i1 = k + 1
for i = i1 to N
ii = i - 1
d = A[kk,ii]
A[kk,ii] = 0.0
for j = 1 to N
jj = j - 1
A[kk,jj] = A[kk,jj] - d * A[ii,jj]
next j
next i
next k

return




all the best,

Graham

[ZeroDog]

Interesting.  I just have one question... when would you need a function like this?   I can see that it could be useful,  I just dont know when you would use it... 3D transformations or something perhaps?   What were you working on when you put this together?

[GWS]

Hi ZD .. nice to see you around again .. 

Coo! .. some folk want cream with it  .. not content to just admire a brand new Inverse, they want to know what use is it?   ..

Now I've got to delve back umpteen years when I understood this stuff ... Let's see  ..

The simplest application is to solve linear equations.  Take an easy one ..

Solve the simultaneous equations:

x + 2y = 4
3x − 5y = 1

(when this type of problem gets to involve more than 3 variables, it gets a lot more difficult)

We can write this as: A * B = C  (where A is the coefficient array, and B and C are vector arrays)

Array A =  ( 1        2  )      B = ( x )   , and C = ( 4)
                ( 3       -5 )             ( y )                 ( 1)

To solve the equations for x and y, we can write:     B = invert(A) * C

So here's a little test program to do the job ..

Code Select Expand


openconsole
cls

def N:int
' set an example array to hold 2x2 elements (indices starting at 0)..
N = 2

def A[N,N]:double
def B[N],C[N]:double

setprecision 5

' load the coefficient array in column order ..
' (this is necessary because Creative stores array values in column major order).
' (when initialising, the second number refers to the column).
A[0,0] = 1, 3
A[0,1] = 2,-5

' so the actual coefficient array A[ ] is ..
' 1  2
' 3 -5

' load the known constants array C[ ] ..
C[0] = 4.0
C[1] = 1.0

' calculate the inverse of array A[ ] , storing the inverse in the same array ..
invert(A[])

' now multiply the coefficient array C[ ] by the Inverted array A[ ], storing the result in array B[ ] ..
' Note: array indices start at zero ..
B[0] = A[0,0] * C[0] + A[0,1] * C[1]
B[1] = A[1,0] * C[0] + A[1,1] * C[1]

' print the variable values stored in B[ ] ..

print "The value of x is: ", B[0]
print
print "The value of y is: ", B[1]
print


do:until inkey$<>""
closeconsole
end

sub invert(A[]:double)
' routine to invert an NxN array A[ ] inplace (ie. the inverse is formed using the array A[ ] space).
def i,j,k,ik,i1:int
def kk,jj,ii:int
def d:double

for k = 1 to N
kk = k - 1
d = A[kk,kk]
A[kk,kk] = 1.0
for j = 1 to N
jj = j - 1
A[kk,jj] = A[kk,jj] / d
next j
if (k = N) then goto nextpart
ik = k + 1
for i = ik to N
ii = i - 1
d = A[ii,kk]
A[ii,kk] = 0.0
for j = 1 to N
jj = j - 1
A[ii,jj] = A[ii,jj] - D * A[kk,jj]
next j
next i
next k

nextpart:
ik = N - 1
for k = 1 to ik
kk = k - 1
i1 = k + 1
for i = i1 to N
ii = i - 1
d = A[kk,ii]
A[kk,ii] = 0.0
for j = 1 to N
jj = j - 1
A[kk,jj] = A[kk,jj] - d * A[ii,jj]
next j
next i
next k

return


This gives the results x = 2, and  y = 1

You can check this is correct by back-substitution in the original equations.

That example was easy, but when you get say:

2x - 3y + 2z = 14
4x + 4y - 3z = 6
3x + 2y - 3z = -2

things get a bit more difficult.

Form the coefficient matrix ..

A =  2  -3  2              and the C[ ] vector = 14
        4   4  -3                                             6
        3   2  -3                                            -2

invert A[ ] and multiply C[ ] by the inverse  to get B[ ] .

The answer comes out to B = 4
                                           2
                                           6

so this system of equations is satisfied by:  x = 4, y = 2, and z = 6.

What other practical uses might there be?

Take an electric circuit:  (this is hard to do without a picture ..  )

      __ 24volts ___________________     
      |                         |                     |
      2 ohms           3 ohms             6 ohms
      |                         |                     |
      i1                       i2                    i3
      |                         |                     |
      ----------------------------------------

it's supposed to be a battery supplying 3 circuit branches, each having a resistance.
The currents I1, I2 and I3 are assumed flowing in a downwards direction on the diagram.

Using Kirchhoff's circuit laws, we get:
I1     + I2   + I3   = 0
-2I1  + 3I2          = 24
        - 3I2 + 6I3  = 0

From which we can form the coefficient matrix A[ ] =   1    1    1
                                                                              -2    3    0
                                                                               0   -3    6
and the constants array C[ ] =     0
                                               24
                                                 0

Again we invert A[ ] and multiply C[ ] by it to get the solution array B[ ] for the currents in each branch I1, I2, and I3.

There are lots of practical problems where arrays can be used in this way.

I have another example of static forces in a mechanical system - and another of input/output analysis for a multi-product firm to calculate the optimum product mix.

Lots of fun for everyone ..

Hope that helps, 

Graham



               

[LarryMc]

Graham
It's been many years for me but for the life of me I don't understand.
This is the way I'm seeing it.

I1     + I2   + I3   = 0
-2I1  + 3I2          = 24
        - 3I2 + 6I3  = 0

in your example
With the values given in your diagram it is obvious(to me) that (using I=E/R)
i1=12 amps
i2= 8 amps
i3= 4 amps

so that I1     + I2   + I3   =24 amps not 0
that -2I1  + 3I2 is -2*12 +3*4 which is -24 +12 which equals -12 and not 24
and that  - 3I2 + 6I3   is -3*8+6*4 which equals 0 and is correct.

What am I missing or doing wrong?

Larry

[GWS]

Hi Larry, it's probably my poor attempt at drawing the circuit with just the keyboard ..   

The battery is not across all the resistors in parallel - it's just in series in the left-hand loop.

The 6 and 3 ohm resistors are in parallel, and this group is in series with the 2 ohm resistor.

The battery supplies this series / parallel arrangement.  So I make it that I1 = -6 amps, I2 = 4 amps, and I3 = 2 amps .. 
I have to apologise for the rather ragged code - it's a resurrection of a 30 year old Fortran routine
that I used often in my work.  All I can say in it's defence - is that it seems to work .. 

Fortran used to have constructs like:

       DO 180 I = 1 , N
        ...
          DO 180 J = 0 , 3
           ....

              DO 180 K = 1 , J
                ......

180 CONTINUE

and that takes a bit of thought to convert .. 

all the best, 

Graham

[LarryMc]

With your explanation of the diagram it makes perfectly good sense to me now.

I'm glad it was something that simple as opposed to me "going off the deep end".

Seems I'm having more and more "senior moments"

Thanks

Larry