Generalization: The Pseudoinverse - Executable Engineering: Interactive Numerical Methods in Python

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

import numpy as np

In previous sections, we solved square matrix systems ( $n$ equations, $n$ unknowns). But what if we have more equations than unknowns?

\mathbf{A} \mathbf{x} = \mathbf{b}

(1)

where $\mathbf{A}$ is an $m \times n$ rectangular matrix (with $m > n$ ). Such systems are called overdetermined (more equations than variables), and we usually cannot find a single solution $\mathbf{x}$ that perfectly satisfies every equation simultaneously. Instead, we can expand our notion of finding a solution to finding a best-fit.

Minimizing the Error¶

Define the residual, $\mathbf{r}$ , of the system as the difference between the left and right side of the equation:

\mathbf{r} = \mathbf{A}\mathbf{x} - \mathbf{b}

(2)

If an exact solution existed, $\mathbf{r}=\mathbf{0}$ (the equation is satisfied for each element).

In lieu of an exact solution we can find one that minimizes all the elements of $\mathbf{r}$ , i.e. minimizes the quantity:

\frac{1}{2} \mathbf{r} \cdot \mathbf{r} = \frac{1}{2} (\mathbf{A}\mathbf{x}-\mathbf{b})^T (\mathbf{A}\mathbf{x}-\mathbf{b})

(3)

where the factor of $\frac{1}{2}$ is arbitrary but useful as seen below.

The Normal Equations¶

To find the minimum of our residul quantity $\frac{1}{2} \mathbf{r} \cdot \mathbf{r}$ , we take the derivative with respect to $\mathbf{x}$ and set it to zero:

\begin{aligned} \frac{1}{2} \frac{\mathrm{d} \mathbf{r} \cdot \mathbf{r}}{\mathrm{d}\mathbf{x}} &= 0 \\ \mathbf{A}^T (\mathbf{A}\mathbf{x} - \mathbf{b}) &= 0 \\ (\mathbf{A}^T \mathbf{A}) \mathbf{x} &= \mathbf{A}^T \mathbf{b} \end{aligned}

(4)

This final equation is called the Normal Equations. Notice that even though $\mathbf{A}$ is rectangular ( $m \times n$ ), the product $(\mathbf{A}^T \mathbf{A})$ is always a square $n \times n$ matrix!

The Normal Equations linear system with a square matrix and can be solved using our usual toolkit for linear systems. However, we can go a step further.

The Moore-Penrose Pseudoinverse¶

Since $(\mathbf{A}^T \mathbf{A})$ is a square matrix, we can invert it to solve for $\mathbf{x}$ :

\mathbf{x} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{b}

(5)

We define the pseudoinverse, denoted $\mathbf{A}^+$ , as:

\mathbf{A}^+ = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T

(6)

So our solution is simply $\mathbf{x} = \mathbf{A}^+ \mathbf{b}$ .

The pseudoinverse generalizes the inverse to rectangular matrices (both overdetermined and underdetermined). If the matrix actually is square and invertible, the pseudoinverse reduces mathematically to the standard exact inverse!

In practice, the $\mathbf{A}^T \mathbf{A}$ matrix, while square, is ill-conditioned. The pseudoinverse is calculated using a matrix decomposition called singular value decomposition, a generalization of eigendecomposition for rectangular matricies, which is more stable.

The number function is ‘np.linalg.pinv(A)’

Terminology¶

Consistent System: A system that has at least one exact solution satisfying all equations.
Inconsistent System: A system with solution that satisfies all equations simultaneously.
- Overdetermined: More equations than unknowns ( $m > n$ ). Usually inconsistent. We use the pseudoinverse to find the best fit solution that minimizes the error.
- Underdetermined: Fewer equations than unknowns ( $m < n$ ). Has an infinite number of solutions. The pseudoinverse finds the solution with the smallest magnitude: $\mathbf{x}\cdot \mathbf{x}$ .

Example: Consistent Overdetermined System¶

Let’s look at a system with 3 equations but only 2 variables.

(1) $20 c + 50 t = 700$
(2) $c + t = 20$
(3) $50 c + 20 t = 700$

\begin{pmatrix} 20 & 50 \\ 1 & 1 \\ 50 & 20 \end{pmatrix} \begin{pmatrix} c \\ t \end{pmatrix} = \begin{pmatrix} 700 \\ 20 \\ 700 \end{pmatrix}

(7)

Since it’s still a linear system, lets inspect this graphically:

A = [[20, 50], [1, 1], [50, 20]]
b = [700, 20, 700]

fig = exe.visual_solve_2d(A, b, method='pseudo')
fig.show()

Example: Inconsistent Overdetermined System¶

Usually, 3 lines will not cross at exactly the same point.

(1) $2 c + t = 15$
(2) $c + -t = 10$
(3) $c + 2 t = 14$

A = [[2, 1], [1,-1], [1, 2]]
b = [15, 10, 14]

fig = exe.visual_solve_2d(A, b, method='pseudo')
fig.show()

The solution found by the pseudoinverse minimizes the perpendicular distance from the point to all three lines simultaneously. It mathematically finds the “center of mass” of the intersection triangle!

Example: Inconsistent overdetermined system with linearly dependent rows¶

Let’s look at one more inconsistent system.

(1) $2 c + t = 15$
(2) $c + -t = 10$
(3) $2 c + 1 t = 10$

Now two equations are parallel (but not overlapping).

A = [[2, 1], [1,-1], [2, 1]]
b = [15, 10, 10]

fig = exe.visual_solve_2d(A, b, method='pseudo')
fig.show()

Example: Inconsistent underdetermined system¶

Consider if we only had one equation but two unknowns!

(1) $2 c + t = 15$

Let’s see what ‘solution’ the pseudoinverse finds!

A = [[2, 1]]
b = [15]

fig = exe.visual_solve_2d(A, b, method='pseudo')
fig.show()

Conditioning of a Rectangular Matrix¶

Just as we saw with square matrices, we can calculate the condition number of a rectangular matrix to understand how sensitive our best-fit solution is to noise in the data $\mathbf{b}$ .

We simply extend the definition using the pseudoinverse:

\text{cond}(\mathbf{A}) = \|\mathbf{A}\| \|\mathbf{A}^+\|

(8)

Modern approach to pseduoinverse calculation¶

Baring a few special cases, the pseudoinverse is most efficiently solved through Singular Value Decomposition (the generalization of eigen-decomposition for rectangular matricies.)