∑MathIntermediate

Matrix Inverse

Key Points

  • •
    A matrix inverse undoes the effect of a linear transformation, just like dividing by a number undoes multiplication.
  • •
    A square matrix A is invertible if and only if det(A) ≠ 0, which is equivalent to having full rank and linearly independent rows/columns.
  • •
    Computing an inverse numerically is usually done by elimination-based methods (Gauss–Jordan or LU with pivoting), not by cofactors/adjugate.
  • •
    The cost to invert an n × n matrix is Θ() time and Θ() space with classical algorithms.
  • •
    In floating point, partial pivoting is essential to reduce numerical error; avoid explicitly forming unless truly needed.
  • •
    To solve , prefer LU (or QR) solves over explicitly computing ; it is faster and more stable.
  • •
    Over finite fields (e.g., modulo a prime p), the same elimination steps work by replacing division with modular inverses.
  • •
    Key identities include (AB)^{-1} = , ()^{-1} = ()^T, and det() = 1/det(A).

Prerequisites

  • →Systems of linear equations — Matrix inversion is closely tied to solving Ax = b; understanding linear systems clarifies why and when inverses exist.
  • →Matrix multiplication and identity — The definition of an inverse relies on matrix multiplication and the identity matrix I.
  • →Determinant and rank — Invertibility conditions det(A) ≠ 0 and rank(A) = n are fundamental tests.
  • →Row operations and Gaussian elimination — Practical inversion uses elimination; familiarity with row operations is essential.
  • →Floating-point arithmetic and numerical stability — Pivoting and conditioning matter when computing inverses with real numbers.
  • →Modular arithmetic (finite fields) — Over 𝔽_p, division becomes modular inversion; elimination steps adapt directly.
  • →Algorithmic complexity (Big-O) — Helps compare inversion with solving and understand Θ(n^3) costs.
  • →Vector and matrix norms — Needed to understand condition numbers and error bounds.

Detailed Explanation

Tap terms for definitions

01Overview

The inverse of a matrix is the matrix analogue of a reciprocal. For a square matrix A, a matrix A^{-1} is called its inverse if multiplying by it from either side returns the identity matrix I: AA^{-1} = A^{-1}A = I. Conceptually, if A represents a linear transformation that rotates, scales, or shears vectors, then A^{-1} reverses that transformation. Not every square matrix is invertible: invertibility requires that A be nonsingular, which has many equivalent characterizations, such as nonzero determinant, full rank, and linearly independent rows and columns. In practical computing, inverses are rarely found via the cofactor/adjugate formula because it is numerically unstable and inefficient for large matrices. Modern algorithms use elimination-based techniques such as Gauss–Jordan elimination or LU factorization with partial pivoting, which share the same asymptotic complexity as solving a linear system. These methods run in Θ(n^3) time and Θ(n^2) space for an n × n dense matrix. In floating-point arithmetic, pivoting and scaling are vital to control round-off error. In modular arithmetic (like programming contests), the same steps work by replacing division with modular inverses modulo a prime p. Understanding matrix inverses underpins solutions to systems of equations, coordinate changes, control and robotics, optimization, statistics, cryptography, and more.

02Intuition & Analogies

Imagine a reversible action: putting on a jacket can be undone by taking it off. If a matrix A is the “put on the jacket” transformation, then A^{-1} is the “take it off” transformation. Apply A, then A^{-1}, and you are back to where you started: that’s AA^{-1} = I. But some actions are not reversible: if you blend two paint colors, you can’t unmix them to get the originals. Similarly, a matrix that squashes 3D space onto a plane loses information and cannot be inverted.

Think of a matrix as a machine that maps input vectors to output vectors. If two different inputs produce the same output, the machine can’t be reversed—you wouldn’t know which input to choose when going backwards. Invertible matrices avoid this problem: they are one-to-one and onto, so every output corresponds to exactly one input. Geometrically, invertible matrices stretch, rotate, reflect, or shear space without collapsing any dimension to zero thickness. The determinant measures the “volume scaling” of the transformation. If det(A) = 0, the transformation flattens space along some direction, like pressing dough with a rolling pin—information is lost, and no perfect undo exists. If det(A) ≠ 0, there’s no flattening, and an undo is possible.

In computation, finding the inverse is like solving the equation Ax = b for many different right-hand sides b. Instead of inverting outright, we often solve those systems directly—which is like using the machine in reverse without physically building a second machine. Over integers modulo a prime, undoing is still possible, but “division” becomes multiplying by a modular inverse instead of a real-number reciprocal.

03Formal Definition

Let A ∈ be a square matrix over a field F (e.g., ℝ, ℂ, or a finite field 𝔽_p). A is invertible (nonsingular) if there exists B ∈ such that A, where is the identity matrix. The inverse, denoted , is unique when it exists. Equivalent conditions for invertibility include: det(A) ≠ 0; rank(A) = n; the columns (and rows) of A are linearly independent; the linear map x ↦ Ax is bijective; zero is not an eigenvalue of A; all singular values are strictly positive. Key algebraic identities include: (AB)^{-1} = for invertible A, B; ()^{-1} = ()^T; (A^*)^{-1} = ()^* where * denotes conjugate transpose; det() = 1/det(A); and ()^{-1} = A. The classical (but impractical) adjugate formula is = (1/ det(A)) adj(A), where adj(A) is the transpose of the cofactor matrix. In numerical linear algebra, one usually computes an LU factorization P with a permutation matrix P from partial pivoting. Then one can solve by forward/back substitution or compute by solving n systems for the columns of the identity. Over finite fields 𝔽_p (p prime), the same elimination steps define inverses by replacing reciprocals with modular inverses.

04When to Use

  • Solving linear systems: Conceptually, x = A^{-1}b. Practically, do not form A^{-1}; factor A (e.g., LU or QR) and solve. This is faster, more accurate, and avoids extra rounding error.
  • Coordinate transforms: Changing bases or frames in graphics, robotics, and computer vision often requires applying the inverse transformation to map coordinates back.
  • Control and estimation: State estimation (e.g., Kalman filters) and control laws frequently need solving linear systems and sometimes explicit inverses of small matrices (2×2, 3×3) where closed forms are convenient.
  • Symbolic or exact arithmetic: Over rationals or modulo a prime p (cryptography, coding theory, competitive programming), forming A^{-1} with Gauss–Jordan over exact arithmetic can be appropriate.
  • Precomputation: If the same matrix A is reused to transform many vectors and memory allows, storing an explicit A^{-1} can amortize costs and reduce per-query latency.
  • Block methods and low-rank updates: Identities like the Sherman–Morrison–Woodbury formula efficiently update inverses when A changes by a low-rank term.
    Avoid inverting when you only need Ax = b a few times numerically; instead, factor once and solve. Inverse is also inappropriate when A is ill-conditioned; prefer regularization or more stable factorizations (QR/SVD).

⚠️Common Mistakes

  • Computing A^{-1} when you only need to solve Ax = b. This wastes time and increases rounding error. Prefer LU/QR solves.
  • Skipping pivoting in floating-point elimination. Without partial pivoting, tiny pivots can explode rounding errors; always select a suitable pivot (max magnitude in column).
  • Using the adjugate/cofactor formula for n > 3. It is slow (superlinear factorial-like growth) and numerically unstable due to catastrophic cancellation.
  • Treating near-singular matrices as invertible. A tiny determinant (or huge condition number Îş(A)) means solutions are highly sensitive to errors. Check conditioning or use SVD and regularization.
  • In modular arithmetic, forgetting that division is not defined. Replace division by multiplication with a modular inverse, and ensure pivots are nonzero modulo p; otherwise, swap rows.
  • Confusing (AB)^{-1} with A^{-1}B^{-1}. The correct identity is (AB)^{-1} = B^{-1}A^{-1}; the order reverses.
  • Not propagating row swaps to the right-hand side when using LU with pivoting, leading to incorrect solutions.
  • Using a single global epsilon for all scales. Scale-aware checks (relative to row/column norms) are more robust than a fixed threshold.

Key Formulas

Inverse definition

Explanation: Multiplying a matrix by its inverse (on either side) yields the identity matrix, meaning the inverse undoes A.

Determinant test

Explanation: A nonzero determinant indicates no dimension is collapsed, so the matrix can be reversed.

Rank test

Explanation: Full rank means columns (and rows) are linearly independent, ensuring a unique inverse exists.

Adjugate formula

Explanation: Theoretical formula for inverses using cofactors; rarely used in practice for n > 3 due to instability and cost.

Product inverse

Explanation: The order reverses when inverting a product; used in algebraic manipulations and block factorizations.

Transpose/adjoint inverse

Explanation: Inverses commute with transpose and conjugate transpose; useful in symmetric/Hermitian contexts.

Determinant of inverse

Explanation: The volume scaling of the inverse is the reciprocal of the original scaling.

Condition number

Explanation: Measures sensitivity: a larger Îş means solutions x to Ax = b change more with small changes in A or b.

LU factorization cost

Explanation: Classical LU with partial pivoting requires about floating-point operations for dense matrices.

System solution via inverse

Explanation: Conceptual solution to in practice, we compute x by solving with factorizations rather than forming .

Woodbury identity

Explanation: Updates the inverse efficiently when A receives a low-rank modification UCV.

Classical inversion complexity

Explanation: Dense matrix inversion with elimination-based methods runs in time proportional to n cubed.

Involution of inverse

Explanation: The inverse of the inverse returns the original matrix.

Complexity Analysis

For dense n × n matrices, all classical elimination-based methods to compute or apply an inverse are Θ() time and Θ() space. Gauss–Jordan elimination augments A with I and performs forward and backward elimination to reduce [A ]. Each pivot step touches O(n) rows of length O(n), and there are n pivots, yielding Θ() arithmetic. Memory is dominated by storing the augmented matrix, Θ(). LU factorization with partial pivoting computes P U in about floating-point operations (flops). If you only need to solve , the subsequent forward/back substitutions cost O() per right-hand side. To explicitly form via LU, you solve n right-hand sides (the columns of I), adding about n · O() = O() work. The total remains Θ() (more precisely around flops). In practice, the constant factors, memory access patterns, and BLAS-3 (cache-friendly) implementations matter greatly. Over finite fields modulo a prime p, the same asymptotic costs apply; modular multiplication/addition are O(1) on typical integer hardware, while modular inverse uses fast exponentiation O(log p), but only for n pivots, so the total is still Θ( + n log p). For sparse matrices, specialized sparse factorizations can be much faster than dense Θ() by exploiting zero structure. Asymptotically faster algorithms exist via fast matrix multiplication (exponent ω < 2.373), giving O(n^ inversion, but they are impractical for most real-world sizes. Space remains Θ() for dense methods, dominated by matrix storage; additional overhead for pivot indices is O(n).

Code Examples

Gauss–Jordan matrix inverse with partial pivoting (double)
1#include <bits/stdc++.h>
2using namespace std;
3
4// Compute the inverse of a square matrix A using Gauss–Jordan elimination
5// with partial pivoting. Returns pair(success, inverse).
6// If singular (or numerically singular), success = false.
7pair<bool, vector<vector<double>>> invertMatrix(const vector<vector<double>>& A) {
8    int n = (int)A.size();
9    const double EPS = 1e-12; // tolerance for detecting near-zero pivots
10
11    // Build augmented matrix [A | I]
12    vector<vector<double>> aug(n, vector<double>(2*n, 0.0));
13    for (int i = 0; i < n; ++i) {
14        for (int j = 0; j < n; ++j) aug[i][j] = A[i][j];
15        aug[i][n + i] = 1.0; // identity on the right
16    }
17
18    // Gauss–Jordan elimination to reduce [A|I] to [I|A^{-1}]
19    for (int col = 0; col < n; ++col) {
20        // Partial pivoting: find row with max |aug[row][col]|
21        int pivot = col;
22        for (int r = col + 1; r < n; ++r) {
23            if (fabs(aug[r][col]) > fabs(aug[pivot][col])) pivot = r;
24        }
25        if (fabs(aug[pivot][col]) < EPS) {
26            return {false, {}}; // singular or ill-conditioned w.r.t. EPS
27        }
28        if (pivot != col) swap(aug[pivot], aug[col]);
29
30        // Normalize pivot row so pivot becomes 1
31        double div = aug[col][col];
32        for (int j = 0; j < 2*n; ++j) aug[col][j] /= div;
33
34        // Eliminate this column in all other rows
35        for (int r = 0; r < n; ++r) if (r != col) {
36            double factor = aug[r][col];
37            if (factor != 0.0) {
38                for (int j = 0; j < 2*n; ++j) aug[r][j] -= factor * aug[col][j];
39            }
40        }
41    }
42
43    // Extract the inverse from the right half
44    vector<vector<double>> inv(n, vector<double>(n, 0.0));
45    for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) inv[i][j] = aug[i][n + j];
46    return {true, inv};
47}
48
49int main() {
50    // Example: invert a 3x3 matrix
51    vector<vector<double>> A = {
52        {1, 2, 3},
53        {0, 1, 4},
54        {5, 6, 0}
55    };
56
57    auto [ok, inv] = invertMatrix(A);
58    if (!ok) {
59        cout << "Matrix is singular or ill-conditioned.\n";
60        return 0;
61    }
62
63    cout.setf(std::ios::fixed); cout << setprecision(6);
64    cout << "Inverse:\n";
65    for (auto &row : inv) {
66        for (double x : row) cout << setw(12) << x << ' ';
67        cout << '\n';
68    }
69
70    // Verify by computing A * A^{-1} ≈ I
71    int n = (int)A.size();
72    vector<vector<double>> I(n, vector<double>(n, 0.0));
73    for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j)
74        for (int k = 0; k < n; ++k) I[i][j] += A[i][k] * inv[k][j];
75
76    cout << "\nProduct A*inv(A) (should be I):\n";
77    for (int i = 0; i < n; ++i) {
78        for (int j = 0; j < n; ++j) cout << setw(12) << I[i][j] << ' ';
79        cout << '\n';
80    }
81}
82

Build the augmented matrix [A | I] and perform Gauss–Jordan elimination with partial pivoting. Each pivot step scales the pivot row and eliminates the pivot column in all other rows. The right block becomes A^{-1} if A is nonsingular. We verify by multiplying A with the computed inverse to approximate the identity.

Time: O(n^3)Space: O(n^2)
Solve Ax = b via LU with partial pivoting (don’t form the inverse)
1#include <bits/stdc++.h>
2using namespace std;
3
4// LU factorization with partial pivoting: produces LU in-place and pivot vector 'piv'.
5bool luDecompose(const vector<vector<double>>& A, vector<vector<double>>& LU, vector<int>& piv) {
6    int n = (int)A.size();
7    const double EPS = 1e-12;
8    LU = A;
9    piv.assign(n, 0);
10
11    for (int k = 0; k < n; ++k) {
12        // Find pivot row
13        int p = k;
14        for (int r = k + 1; r < n; ++r) if (fabs(LU[r][k]) > fabs(LU[p][k])) p = r;
15        if (fabs(LU[p][k]) < EPS) return false; // singular/ill-conditioned
16        if (p != k) swap(LU[p], LU[k]);
17        piv[k] = p;
18        // Eliminate below
19        for (int i = k + 1; i < n; ++i) {
20            LU[i][k] /= LU[k][k];
21            for (int j = k + 1; j < n; ++j) LU[i][j] -= LU[i][k] * LU[k][j];
22        }
23    }
24    return true;
25}
26
27// Solve Ax = b using LU (with pivot vector piv recording row swaps).
28vector<double> luSolve(const vector<vector<double>>& LU, const vector<int>& piv, const vector<double>& b) {
29    int n = (int)LU.size();
30    vector<double> x = b;
31    // Apply the recorded row swaps to b: Pb
32    for (int k = 0; k < n; ++k) if (piv[k] != k) swap(x[k], x[piv[k]]);
33    // Forward substitution for Ly = Pb (L has unit diagonal, L entries are LU[i][j] for i>j)
34    for (int i = 0; i < n; ++i) {
35        for (int j = 0; j < i; ++j) x[i] -= LU[i][j] * x[j];
36        // since diag(L)=1, no division here
37    }
38    // Back substitution for Ux = y (U on/above diagonal in LU)
39    for (int i = n - 1; i >= 0; --i) {
40        for (int j = i + 1; j < n; ++j) x[i] -= LU[i][j] * x[j];
41        x[i] /= LU[i][i];
42    }
43    return x;
44}
45
46int main() {
47    // Example system Ax = b
48    vector<vector<double>> A = {
49        {3, 1, -1},
50        {2, 4, 1},
51        {-1, 2, 5}
52    };
53    vector<double> b = {4, 1, 1};
54
55    vector<vector<double>> LU; vector<int> piv;
56    if (!luDecompose(A, LU, piv)) {
57        cout << "Matrix is singular or ill-conditioned.\n";
58        return 0;
59    }
60
61    vector<double> x = luSolve(LU, piv, b);
62    cout.setf(std::ios::fixed); cout << setprecision(6);
63    cout << "Solution x to Ax=b (without forming A^{-1}):\n";
64    for (double xi : x) cout << xi << ' ';
65    cout << '\n';
66
67    // Check residual ||Ax - b||
68    vector<double> r(b.size(), 0.0);
69    for (size_t i = 0; i < A.size(); ++i) for (size_t j = 0; j < A.size(); ++j) r[i] += A[i][j] * x[j];
70    for (size_t i = 0; i < r.size(); ++i) r[i] -= b[i];
71    double norm = 0.0; for (double v : r) norm = max(norm, fabs(v));
72    cout << "Max-abs residual: " << norm << '\n';
73}
74

Factor A as PA = LU with partial pivoting, then solve Ly = Pb and Ux = y. This avoids explicitly computing the inverse, which is recommended for performance and numerical stability when solving systems.

Time: O(n^3) to factor, plus O(n^2) per right-hand sideSpace: O(n^2)
Matrix inverse modulo a prime p (Gauss–Jordan over finite field)
1#include <bits/stdc++.h>
2using namespace std;
3using int64 = long long;
4
5const int64 MOD = 1000000007LL; // must be prime for Fermat-based inverses
6
7int64 normmod(int64 x) { x %= MOD; if (x < 0) x += MOD; return x; }
8
9int64 modmul(int64 a, int64 b) { return ( (__int128)normmod(a) * normmod(b) ) % MOD; }
10
11int64 modpow(int64 a, int64 e) {
12    int64 r = 1 % MOD; a = normmod(a);
13    while (e > 0) {
14        if (e & 1) r = modmul(r, a);
15        a = modmul(a, a);
16        e >>= 1;
17    }
18    return r;
19}
20
21int64 modinv(int64 a) { return modpow(a, MOD - 2); } // Fermat's little theorem (MOD prime)
22
23pair<bool, vector<vector<int64>>> invertMatrixMod(vector<vector<int64>> A) {
24    int n = (int)A.size();
25    vector<vector<int64>> aug(n, vector<int64>(2*n, 0));
26    for (int i = 0; i < n; ++i) {
27        for (int j = 0; j < n; ++j) aug[i][j] = normmod(A[i][j]);
28        aug[i][n + i] = 1;
29    }
30
31    for (int col = 0; col < n; ++col) {
32        // Find a nonzero pivot in this column (any will do over a field)
33        int pivot = -1;
34        for (int r = col; r < n; ++r) if (aug[r][col] != 0) { pivot = r; break; }
35        if (pivot == -1) return {false, {}}; // singular modulo MOD
36        if (pivot != col) swap(aug[pivot], aug[col]);
37
38        // Scale pivot row to make pivot == 1
39        int64 inv = modinv(aug[col][col]);
40        for (int j = 0; j < 2*n; ++j) aug[col][j] = modmul(aug[col][j], inv);
41
42        // Eliminate this column in all other rows
43        for (int r = 0; r < n; ++r) if (r != col) {
44            int64 factor = aug[r][col];
45            if (factor != 0) {
46                for (int j = 0; j < 2*n; ++j) {
47                    aug[r][j] = normmod(aug[r][j] - modmul(factor, aug[col][j]));
48                }
49            }
50        }
51    }
52
53    vector<vector<int64>> inv(n, vector<int64>(n));
54    for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) inv[i][j] = aug[i][n + j];
55    return {true, inv};
56}
57
58int main() {
59    // Example over MOD prime: invert a 2x2 matrix
60    vector<vector<int64>> A = {{4, 7}, {2, 6}}; // det = 4*6 - 7*2 = 10 (invertible mod 1e9+7)
61    auto [ok, inv] = invertMatrixMod(A);
62    if (!ok) { cout << "Singular modulo MOD.\n"; return 0; }
63    cout << "Inverse modulo " << MOD << ":\n";
64    for (auto &row : inv) {
65        for (auto v : row) cout << v << ' ';
66        cout << '\n';
67    }
68
69    // Verify A * inv(A) ≡ I (mod MOD)
70    int n = (int)A.size();
71    vector<vector<int64>> I(n, vector<int64>(n, 0));
72    for (int i = 0; i < n; ++i)
73        for (int j = 0; j < n; ++j)
74            for (int k = 0; k < n; ++k)
75                I[i][j] = normmod(I[i][j] + (__int128)A[i][k] * inv[k][j]);
76
77    cout << "Product A*inv(A) mod MOD:\n";
78    for (int i = 0; i < n; ++i) {
79        for (int j = 0; j < n; ++j) cout << I[i][j] << ' ';
80        cout << '\n';
81    }
82}
83

Performs Gauss–Jordan elimination over the finite field 𝔽_p. Division is implemented by multiplying with modular inverses (Fermat’s little theorem since MOD is prime). Row swaps handle zero pivots. The result satisfies A·A^{-1} ≡ I (mod p).

Time: O(n^3) + O(n log p) for n modular inversesSpace: O(n^2)
Closed-form 2×2 inverse (safe check with determinant)
1#include <bits/stdc++.h>
2using namespace std;
3
4bool invert2x2(const array<array<double,2>,2>& A, array<array<double,2>,2>& inv) {
5    double a = A[0][0], b = A[0][1], c = A[1][0], d = A[1][1];
6    double det = a*d - b*c;
7    const double EPS = 1e-12;
8    if (fabs(det) < EPS) return false; // singular or nearly singular
9    double id = 1.0 / det;
10    inv = {array<double,2>{ d*id, -b*id }, array<double,2>{ -c*id, a*id }};
11    return true;
12}
13
14int main() {
15    array<array<double,2>,2> A = {{{4, 7}, {2, 6}}};
16    array<array<double,2>,2> inv;
17    if (!invert2x2(A, inv)) { cout << "Singular.\n"; return 0; }
18    cout.setf(std::ios::fixed); cout << setprecision(6);
19    cout << "Inverse of 2x2:\n";
20    cout << inv[0][0] << ' ' << inv[0][1] << '\n';
21    cout << inv[1][0] << ' ' << inv[1][1] << '\n';
22}
23

Demonstrates the classical adjugate formula for a 2×2 matrix: A^{-1} = (1/(ad−bc)) [[d, −b], [−c, a]]. For very small matrices, the closed form is convenient and efficient; always check det ≠ 0.

Time: O(1)Space: O(1)