C Program To Find Inverse Of A Matrix

Finding the inverse of a matrix is a fundamental operation in linear algebra, with wide-ranging applications in various fields. In this article, you will learn how to implement a C program to calculate the inverse of a square matrix using the determinant and adjoint method.

Problem Statement

A matrix inverse exists only for square matrices (n x n) that are non-singular, meaning their determinant is not zero. Given a square matrix A, its inverse, denoted A⁻¹, is another square matrix such that A * A⁻¹ = I, where I is the identity matrix. The challenge is to compute A⁻¹ given A, especially considering floating-point arithmetic and the conditions for existence.

Example

Consider a simple 2x2 matrix:

A = | 4  7 |
    | 2  6 |

Its determinant is (4*6) - (7*2) = 24 - 14 = 10. The adjoint matrix is:

Adj(A) = | 6  -7 |
         | -2  4 |

The inverse A⁻¹ is (1/Determinant) * Adj(A):

A⁻¹ = (1/10) * | 6  -7 |   = | 0.6  -0.7 |
               | -2  4 |     | -0.2  0.4 |

Background & Knowledge Prerequisites

To understand and implement a matrix inverse program, readers should be familiar with:

• C Programming Basics: Variables, data types, loops (for, while), conditional statements (if-else), functions, arrays (especially 2D arrays).
• Linear Algebra Concepts:
• Matrices: Understanding rows, columns, and elements.
• Determinant: How to calculate the determinant of 2x2 and 3x3 matrices, and the concept of cofactor expansion for larger matrices.
• Cofactor: The signed minor of a matrix element.
• Adjoint Matrix: The transpose of the cofactor matrix.
• Identity Matrix: A square matrix with ones on the main diagonal and zeros elsewhere.
• Singular Matrix: A matrix with a determinant of zero, which does not have an inverse.

Use Cases or Case Studies

Matrix inversion is crucial in several practical domains:

• Solving Systems of Linear Equations: If you have a system Ax = B, where A is a matrix of coefficients, x is a vector of unknowns, and B is a vector of constants, then x = A⁻¹B.
• Computer Graphics: Used in transformations (rotation, scaling, translation) to reverse operations or convert between different coordinate systems.
• Statistics and Econometrics: Essential for calculating regression coefficients in multivariate analysis, especially in the least squares method.
• Cryptography: Some encryption algorithms rely on matrix operations, where decryption involves inverse matrices.
• Robotics: Used in kinematics to determine joint angles from end-effector positions (inverse kinematics).

Solution Approaches

For calculating the inverse of a matrix, common approaches include:

1. Determinant and Adjoint Method: Suitable for smaller matrices (2x2, 3x3). It involves calculating the determinant, the cofactor matrix, the adjoint matrix (transpose of cofactor), and then scaling by 1/determinant.
2. Gaussian Elimination (Gauss-Jordan Elimination): More general and computationally efficient for larger matrices. This method transforms the augmented matrix [A|I] into [I|A⁻¹] through a series of elementary row operations.
3. LU Decomposition: Decomposes a matrix A into a lower triangular matrix L and an upper triangular matrix U (A = LU). Solving Ly = B and Ux = y is more efficient than direct inversion for multiple right-hand sides.

For this article, we will focus on the Determinant and Adjoint Method as it clearly demonstrates the underlying mathematical principles for smaller matrices.

Finding Inverse Using Determinant and Adjoint

The inverse of a matrix A is given by the formula: A⁻¹ = (1/det(A)) * adj(A), where det(A) is the determinant of A and adj(A) is the adjoint of A.

One-line summary: Calculate the determinant and adjoint of the matrix, then divide each element of the adjoint by the determinant to get the inverse.

Code Example

This C program implements the determinant and adjoint method for a 3x3 matrix.

// Matrix Inverse using Determinant and Adjoint
#include <stdio.h>

#define N 3 // Define the size of the matrix (N x N)

// Function to get cofactor of matrix[p][q] in temp[][]
// n is current dimension of matrix
void getCofactor(double mat[N][N], double temp[N][N], int p, int q, int n) {
    int i = 0, j = 0;
    // Looping for each element of the matrix
    for (int row = 0; row < n; row++) {
        for (int col = 0; col < n; col++) {
            // Copying into temporary matrix only those elements which are not in given row and column
            if (row != p && col != q) {
                temp[i][j++] = mat[row][col];
                // Row is full, so increment row index and reset col index
                if (j == n - 1) {
                    j = 0;
                    i++;
                }
            }
        }
    }
}

// Recursive function for finding determinant of matrix.
// n is current dimension of mat[][].
double determinant(double mat[N][N], int n) {
    double D = 0; // Initialize result

    // Base case : if matrix contains single element
    if (n == 1)
        return mat[0][0];

    double temp[N][N]; // To store cofactors
    int sign = 1;      // To store sign multiplier

    // Iterate for each element of the first row
    for (int f = 0; f < n; f++) {
        // Getting Cofactor of mat[0][f]
        getCofactor(mat, temp, 0, f, n);
        D += sign * mat[0][f] * determinant(temp, n - 1);

        // Terms are to be added with alternate signs
        sign = -sign;
    }
    return D;
}

// Function to calculate and store adjoint of matrix[N][N] in adj[N][N]
void adjoint(double mat[N][N], double adj[N][N]) {
    if (N == 1) {
        adj[0][0] = 1;
        return;
    }

    int sign = 1;
    double temp[N][N];

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            // Get cofactor of mat[i][j]
            getCofactor(mat, temp, i, j, N);

            // sign is (-1)^(i+j)
            sign = ((i + j) % 2 == 0) ? 1 : -1;

            // Transpose of cofactor matrix is adjoint
            adj[j][i] = (sign) * (determinant(temp, N - 1));
        }
    }
}

// Function to calculate and store inverse, returns false if
// matrix is singular
int inverse(double mat[N][N], double inv[N][N]) {
    // Find determinant of matrix
    double det = determinant(mat, N);
    if (det == 0) {
        printf("Singular matrix, cannot find its inverse.\\n");
        return 0;
    }

    // Find adjoint
    double adj[N][N];
    adjoint(mat, adj);

    // Find inverse using formula "inverse(A) = adj(A)/det(A)"
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            inv[i][j] = adj[i][j] / det;
        }
    }
    return 1;
}

// Generic function to display a matrix
void display(double mat[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%8.4f ", mat[i][j]);
        }
        printf("\\n");
    }
}

int main() {
    // Step 1: Define the matrix for which to find the inverse
    double mat[N][N] = {{1, 2, 3},
                        {4, 5, 6},
                        {7, 8, 9}}; // This matrix is singular (determinant is 0)

    double mat2[N][N] = {{1, 2, 3},
                         {0, 1, 4},
                         {5, 6, 0}}; // This matrix is non-singular

    double inverseMat[N][N];

    printf("Original Matrix 1:\\n");
    display(mat);

    // Step 2: Attempt to find the inverse of the first matrix
    if (inverse(mat, inverseMat)) {
        printf("\\nInverse of Matrix 1:\\n");
        display(inverseMat);
    } else {
        printf("\\nMatrix 1 is singular, inverse does not exist.\\n");
    }

    printf("\\n-------------------------\\n");

    printf("\\nOriginal Matrix 2:\\n");
    display(mat2);

    // Step 3: Attempt to find the inverse of the second matrix
    if (inverse(mat2, inverseMat)) {
        printf("\\nInverse of Matrix 2:\\n");
        display(inverseMat);
    } else {
        printf("\\nMatrix 2 is singular, inverse does not exist.\\n");
    }

    return 0;
}

Run

Sample Output

Original Matrix 1:
  1.0000   2.0000   3.0000 
  4.0000   5.0000   6.0000 
  7.0000   8.0000   9.0000 
Singular matrix, cannot find its inverse.

Matrix 1 is singular, inverse does not exist.

-------------------------

Original Matrix 2:
  1.0000   2.0000   3.0000 
  0.0000   1.0000   4.0000 
  5.0000   6.0000   0.0000 

Inverse of Matrix 2:
-24.0000  18.0000   5.0000 
 20.0000 -15.0000  -4.0000 
 -5.0000   4.0000   1.0000

Stepwise Explanation for Clarity

1. getCofactor(mat, temp, p, q, n) Function:
   • This helper function generates a submatrix (cofactor matrix) by excluding the element at mat[p][q].

It iterates through the original matrix mat and copies elements into temp only if they are not in row p or column q. This temp matrix will have a dimension of (n-1) x (n-1).

1. determinant(mat, n) Function:
   • This is a recursive function to calculate the determinant of an n x n matrix.

Base Case: If n is 1 (a 1x1 matrix), its determinant is simply the single element.

Recursive Step: For larger matrices, it iterates through the first row. For each element mat[0][f]:

It finds the cofactor matrix using getCofactor.

It recursively calls determinant on this cofactor matrix.

The contribution to the total determinant is sign * mat[0][f] * determinant(cofactor_matrix).

The sign alternates between +1 and -1 for each term.

1. adjoint(mat, adj) Function:
   • Calculates the adjoint of the input matrix mat and stores it in adj.

For a 1x1 matrix, the adjoint is [[1]].

For larger matrices, it iterates through each element mat[i][j]:

It finds the cofactor matrix temp for mat[i][j].

It calculates the determinant of this temp matrix.

The sign for the cofactor is determined by (-1)^(i+j).

The adjoint element adj[j][i] is set to (signed cofactor determinant). Notice adj[j][i] instead of adj[i][j] – this performs the transpose automatically, as the adjoint is the transpose of the cofactor matrix.

1. inverse(mat, inv) Function:
   • This is the main function to compute the inverse.

It first calculates the determinant of the input matrix mat using the determinant function.

Singular Matrix Check: If the determinant is 0, the matrix is singular, and an inverse does not exist. The function prints a message and returns 0 (false).

If the determinant is non-zero, it calculates the adjoint matrix using the adjoint function.

Finally, it computes the inverse inv by dividing each element of the adjoint matrix by the determinant. It returns 1 (true) indicating success.

1. display(mat) Function:
   • A utility function to print the matrix elements in a formatted way.

1. main() Function:
   • Initializes two 3x3 matrices: one singular (mat) and one non-singular (mat2).

Calls inverse() for both matrices and displays the result (either the inverse or a singular matrix message).

Conclusion

Finding the inverse of a matrix is a vital operation in many computational fields. The determinant and adjoint method provides a clear, mathematically grounded approach suitable for smaller matrices, allowing for a direct understanding of the underlying concepts. While computationally intensive for very large matrices, it serves as an excellent illustrative example for learning linear algebra in a programming context.

Summary

• A matrix inverse A⁻¹ exists only for non-singular square matrices (determinant ≠ 0).
• The formula for the inverse is A⁻¹ = (1/det(A)) * adj(A).
• Calculating the inverse involves:
• Computing the matrix's determinant.
• Finding the cofactor for each element.
• Constructing the adjoint matrix (transpose of the cofactor matrix).
• Dividing each element of the adjoint by the determinant.
• Functions like getCofactor, determinant, and adjoint break down the problem into manageable parts.
• The provided C code demonstrates this method for 3x3 matrices, handling the case of singular matrices.