Solving Ax=b using the Inverse

For invertible matrices, the matrix inverse $A^{-1}$ can be calculated. This matrix inverse offers a direct method for solving the fundamental linear system $Ax = b$.

Remember, the equation $Ax = b$ represents a linear transformation $A$ applied to an unknown vector $x$, resulting in a known vector $b$. Our goal is to find the original vector $x$. If the matrix $A$ is square and invertible, we can think of its inverse $A^{-1}$ as the transformation that "undoes" the effect of $A$.

Consider the system: $Ax = b$

Assuming $A$ is invertible, we can pre-multiply both sides of the equation by $A^{-1}$. It's important to multiply from the left on both sides to maintain equality, as matrix multiplication is generally not commutative.

$A^{-1}(Ax) = A^{-1}b$

Using the associative property of matrix multiplication, we can regroup the terms on the left side:

$(A^{-1}A)x = A^{-1}b$

By the definition of the matrix inverse, we know that $A^{-1}A = I$, where $I$ is the identity matrix:

$Ix = A^{-1}b$

Finally, since multiplying any vector by the identity matrix leaves the vector unchanged ($Ix = x$), we arrive at the solution for $x$:

$x = A^{-1}b$

This elegant result tells us that if we can find the inverse of the matrix $A$, we can find the solution vector $x$ simply by multiplying $A^{-1}$ by the vector $b$.

Conditions for Using the Inverse Method

This method relies entirely on the existence of the matrix inverse $A^{-1}$. Therefore, it's applicable only when:

1. A is a square matrix: The number of equations must equal the number of unknowns (i.e., $A$ must be $n \times n$). An inverse is only defined for square matrices.
2. A is invertible (non-singular): The determinant of $A$ must be non-zero ($\det(A) \neq 0$). If the determinant is zero, the matrix does not have an inverse. In this situation, the system $Ax=b$ either has no solutions or infinitely many solutions, and this method cannot find a unique solution.

If these conditions are met, the system $Ax=b$ has a unique solution given by $x = A^{-1}b$.

Example: Solving a Simple System

Let's apply this to a simple 2x2 system:

$$\begin{align*} 2x_1 + 3x_2 &= 8 \\ 1x_1 + 4x_2 &= 9 \end{align*}$$

In matrix form, this is $Ax = b$, where:

$$A = \begin{bmatrix} 2 & 3 \\ 1 & 4 \end{bmatrix}, \quad x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad b = \begin{bmatrix} 8 \\ 9 \end{bmatrix}$$

First, we check if $A$ is invertible by calculating its determinant: $\det(A) = (2)(4) - (3)(1) = 8 - 3 = 5$. Since the determinant is 5 (which is non-zero), the inverse exists.

In the previous section ("Calculating the Inverse of a Matrix"), we might have found (or can calculate now using the formula for a 2x2 inverse) that:

$$A^{-1} = \frac{1}{\det(A)} \begin{bmatrix} d & -b \\ -c & a \end{bmatrix} = \frac{1}{5} \begin{bmatrix} 4 & -3 \\ -1 & 2 \end{bmatrix} = \begin{bmatrix} 4/5 & -3/5 \\ -1/5 & 2/5 \end{bmatrix}$$

Now, we can find $x$ using the formula $x = A^{-1}b$:

$$x = \begin{bmatrix} 4/5 & -3/5 \\ -1/5 & 2/5 \end{bmatrix} \begin{bmatrix} 8 \\ 9 \end{bmatrix}$$

Performing the matrix-vector multiplication:

$$x = \begin{bmatrix} (4/5)(8) + (-3/5)(9) \\ (-1/5)(8) + (2/5)(9) \end{bmatrix} = \begin{bmatrix} (32/5) - (27/5) \\ (-8/5) + (18/5) \end{bmatrix} = \begin{bmatrix} 5/5 \\ 10/5 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix}$$

The solution is therefore $x_1 = 1$ and $x_2 = 2$. We can quickly verify this: $2(1) + 3(2) = 2+6 = 8$ and $1(1) + 4(2) = 1+8 = 9$. The solution is correct.

Implementation with NumPy

While manual calculation is instructive for small matrices, we typically rely on numerical libraries like NumPy for systems encountered in machine learning. NumPy provides functions to compute the inverse and perform the necessary matrix multiplication efficiently.

import numpy as np

# Define the matrix A and vector b from the example
A = np.array([[2, 3],
              [1, 4]])

b = np.array([[8],   # Define b as a column vector (2x1)
              [9]])

# Check if A is square
rows, cols = A.shape
if rows != cols:
    print("Matrix A must be square to find an inverse.")
else:
    # Calculate the determinant
    det_A = np.linalg.det(A)
    print(f"Determinant of A: {det_A:.2f}")

    # Check if the determinant is close to zero (within machine precision)
    if np.isclose(det_A, 0):
        print("Matrix A is singular (or nearly singular).")
        print("Cannot solve using the inverse method.")
    else:
        # Calculate the inverse of A using np.linalg.inv()
        try:
            A_inv = np.linalg.inv(A)
            print("\nInverse of A (A^{-1}):")
            print(A_inv)

            # Calculate the solution x = A_inv * b
            # Use the @ operator for matrix multiplication
            x = A_inv @ b
            # Alternatively: x = np.dot(A_inv, b)

            print("\nSolution vector x = A^{-1}b:")
            print(x)

            # Verification: Check if A @ x is close to the original b
            print("\nVerification (A @ x):")
            print(A @ x)
            print(f"\nIs A @ x close to b? {np.allclose(A @ x, b)}")

        except np.linalg.LinAlgError:
             # Handle cases where inverse computation fails numerically
             print("NumPy encountered an error computing the inverse.")

Output of the Python code:

Determinant of A: 5.00

Inverse of A (A^{-1}):
[[ 0.8 -0.6]
 [-0.2  0.4]]

Solution vector x = A^{-1}b:
[[1.]
 [2.]]

Verification (A @ x):
[[8.]
 [9.]]

Is A @ x close to b? True

This code replicates our manual calculation steps using NumPy. It first checks if A is square, then calculates the determinant using np.linalg.det(). If the determinant is non-zero, it computes the inverse using np.linalg.inv() and then finds the solution x by multiplying A_inv with b using the @ operator. The np.allclose() function is useful for verifying the result, accounting for potential small floating-point inaccuracies.

A Note on Practical Computation

While solving $Ax=b$ using the explicit inverse $x = A^{-1}b$ is clear and mathematically correct, it's often not the most efficient or numerically stable approach in practice, especially for large matrices common in machine learning applications.

1. Computational Cost: Calculating the inverse $A^{-1}$ typically requires more computational operations (flops) than methods like Gaussian elimination or LU decomposition, which are designed specifically to solve the system $Ax=b$ directly.
2. Numerical Stability: Computer arithmetic uses finite precision (floating-point numbers). The process of computing an inverse can sometimes amplify small precision errors present in the matrix $A$. Direct solving methods are often engineered to minimize the accumulation of such errors, leading to more accurate solutions.

For these reasons, libraries like NumPy and SciPy offer functions such as np.linalg.solve(A, b). This function solves $Ax=b$ for $x$ directly, typically using an efficient and stable underlying algorithm (like LU decomposition). It avoids computing the full inverse $A^{-1}$.

# Example using np.linalg.solve (generally preferred)
try:
    x_solve = np.linalg.solve(A, b)
    print("\nSolution using np.linalg.solve(A, b):")
    print(x_solve)
    print(f"Is solution from solve() close to b? {np.allclose(A @ x_solve, b)}")
except np.linalg.LinAlgError:
    print("\nMatrix A is singular. Cannot solve using np.linalg.solve.")

Output of the np.linalg.solve example:

Solution using np.linalg.solve(A, b):
[[1.]
 [2.]]
Is solution from solve() close to b? True

You should generally prefer np.linalg.solve(A, b) over np.linalg.inv(A) @ b when your primary goal is just to find the solution vector $x$.

However, understanding the inverse method remains essential. It provides fundamental theoretical insight into linear systems, appears in many analytical derivations (for instance, in the derivation of the Normal Equations for linear regression), and is perfectly suitable for smaller systems or situations where the inverse matrix $A^{-1}$ itself is needed for further analysis.