Matrix Inversion

Question

I'm developing a program for class that reads a matrix and then prints it with an identity matrix. Then step by step I have to reduce it into its identity matrix while applying the same calculations onto the identity matrix and therefore getting the inverse of that matrix.

I'm stuck at the row reducing section of the code. I'm trying to create a loop that row reduces an nxn matrix into an identity matrix all in one shot instead of creating one for each column. I need some help with the loops that row reduces the matrix column by column.

Anyone have any ideas. (I'm a beginner at programming, I'm using C language, posted solutions are okay, but I'm looking for an algorithm to take down all the row reducing in one shot).

Answer 1

Wikipedia has nice pseudocode for Gaussian Elimination (which is what you've described).

Gaussian Elimination Pseudocode

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k
    end if
  end for
  if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]
    Now A[i,j] will have the value 1.
    for u := i+1 to m do
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1
  end if
  j := j + 1
end while

Answer 2

You might look at the Wikipedia article on Gauss-Jordan elimination http://en.wikipedia.org/wiki/Gauss%E2%80%93Jordan%5Felimination for ideas on what to do. This is a classic method and if you understand it you will be able to solve the algorithmic part of your problem. There are several example codes there as well.

Answer 3

Why don't you try using: http://www.codeproject.com/KB/cs/CSML.aspx

Answer 4

What you're describing is called Gauss Jordan elimination.

Write the (n x n) matrix A next to an (n x n) identity matrix.

You can perform three row operations:

1. Multiply row i by a constant c != 0
2. Interchange rows i and j
3. Add c times row i to row j

Just be aware that this is NOT a good way to find an inverse, because it suffers from roundoff errors. Without pivoting, you could end up with a bad answer.

1. Start with the first row.
2. Divide every element in the ith row by its diagonal element A[i,i], making it equal to 1.
3. Loop over all the rows below the ith row
4. Multiply the ith row by the A[j,i] and subtract the result from the jth row, making the values below the diagonal equal to 0.
5. Repeat until you have a matrix with all ones on the diagonal and zeroes below the diagonal.
6. Do similar row operations to make the values above the diagonal equal to zero, the same way you did for the rows below the diagonal.

When you see the identity matrix where your original A was, the right matrix will equal the inverse. Do a multiplication to prove to yourself that this is true.

Get a solid example and do it out by hand before you start coding. Write unit tests to prove that your code works once it's done.

Answer 5

drug this up from an old college project... Maybe it'll help

typedef struct matrix
  {
     flag type;        /* type of matrix: {RECTANGULAR, SQUARE, BANDED } */
     ushort num_rows;  /* number of rows in th matrix  */
     union cols 
     {
     ushort num_cols;  /* Number of cols in Rectanular matrix */
     ushort band;      /* Bandwidth in a square, symmetric, banded matrix */
     }
     double *m_val;   /* ptr to start of storage of actual matrix values */
   } MAT;

/* row_swap ---- swap two rows of an nxn UNBANDED  matrix   */
void row_swap(mat, r1, r2)
    MAT *mat;
    ushort r1, r2;
    {
    double *m = mat->m_val;
    ushort n = m->num_row;
    INDEX i;
    double temp;
    for (i=0; i < n; i++)
      SWAP(*(m + n*r1 +i), *(m + n*r2 +i), temp);
    }

/*  reduce --- function to do Gauss- Jordon Reduction on matrix */
static void reduce(m, in, ref)
    MAT *m *in;
    double *ref;
    {
    ushort n = mat->num_row;
    INDEX r,c;
    int sr,sc; 
    double key, mat = m->m_val, inv = in->m_val;

    sr =(ref-mat)/n;
    sc=(ref-mat) % n;
    scalprod(1.0/(*ref), inv+n*sr, n);
    scalprod(1.0/(*ref), ref-sc, n);

    for ( r=0; r LT n; r++)
    {
        if ( r != sr && *(mat+n*r+sc) != 0.0 )
           {
           key = *(mat+n*r+sc);
             for (c=0; c < n; c++)
         {
             *(inv+n*r+c) -= key * *(inv+n*sr+c);
            *(mat+n*r+c) -= key * *(mat+n*sr+c);
         }
           }
        }
    }

/* find_non_zero -- this finds and returns a ptr to next non-zero entry 
 *   below and to the right of the last one */
static double *find_non_zero(start, last, n)
   ushort n;
   double *start, *last;
   {
   double *i = last;              /* i is ptr to matrix entries */
   while (*i == 0.0 && (i-start) < n*n)
      {
      if (n+i > start+n*n)  /*  last row without non-zero entry,  */
         return (-1);        /*  Matrix must be singular           */
      else 
         i += n;            /*  this is not the last row, goto next row */
      }
   if (i >= start + n*n)   /* we have departed the matrix */
     return (0);
   else                    /* we found a non-zero entry, return it */ 
     return(i);
   }


/*   invert    --  function to invert nxn matrix  */
flag invert(m, i)
    MAT *m, *i;
    {
    double *ref, *new, *mat = m->m_val, *inv = i->m_val;
    ushort n = mat->num_row;
    INDEX i, j, row;
    ushort new_row;
    if (det(mat,n) == 0.0) return 0;
    for (i=0; i < n; i++)
     for (j=0; j < n; j++)
      *(inv+n*i+j) = (i == j)? 1.0: 0.0;
    ref = mat;
    for (row=0; row < n; row++)
      {
      new = find_non_zero(mat,ref,n);
      if ( new == -1) 
         {
         scr_print(" This matrix is singular, not invertible ");
         break;
         } 
      new_row = (new-mat)/n;
      if (new_row != row)
         {
         row_swap(mat, new_row, row, n);
         row_swap(inv, new_row, row, n);
         }
      reduce(mat,inv,ref,n);
      ref += n+1;
      }
    }

Answer 6

Pardon my earlier joke about the patent. ;-)
I was merely trying to impress on Wert, the OP, that shy of a few [typically local] tricks aimed at either improving performance or numerical accuracy, the multi-step process for converting a matrix to its reduced row echelon, is pretty much the way to go.

Indeed it is particularly important to use this simple and mechanical process (rather than variations there upon, with various local optimizations), because of the OP's lack of familiarity with the C language. First you walk, then you run...

As explained by others, the assignment seems to be the implementation of the Gauss-Jordan algorithm, whereby we don't only perform all the necessary substitutions on the original matrix, but we also carry these on what starts as the identity matrix, affixed to the right of the original matrix. The result of this process is that when the original matrix (the left half of the "augmented matrix") is converted to the identity matrix (assuming it is non singular), the right side, where we had started with the identity, contains the inverse of the main matrix. Because this inverse matrix contains a summary of the individual operations done to reduce the original matrix, it then can allow to do thins "in one shot" (at least at linear algebra level, for this still involves multiple cell-level arithmetic operations), this could be the source of the confusion and the origin of Wert's quest for the "one shot" grail.