5. Generalized inverse

Michael Friendly

2024年11月24日

Source: vignettes/a5-ginv.Rmd

a5-ginv.Rmd

In matrix algebra, the inverse of a matrix is defined only for square matrices, and if a matrix is singular, it does not have an inverse.

The generalized inverse (or pseudoinverse) is an extension of the idea of a matrix inverse, which has some but not all the properties of an ordinary inverse.

A common use of the pseudoinverse is to compute a ‘best fit’ (least squares) solution to a system of linear equations that lacks a unique solution.

library (matlib )

Construct a square, singular matrix [See: Timm, EX. 1.7.3]

A <-matrix (c (4, 4, -2,
 4, 4, -2,
 -2, -2, 10), nrow=3, ncol=3, byrow=TRUE)
det (A)

## [1] 0

The rank is 2, so inv(A) won’t work

R (A)

## [1] 2

In the echelon form, this rank deficiency appears as the final row of zeros

echelon (A)

## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 0 0 1
## [3,] 0 0 0

inv() will throw an error

try (inv (A))

## Error in Inverse(X, tol = sqrt(.Machine$double.eps), ...) : 
## X is numerically singular

A generalized inverse does exist for any matrix, but unlike the ordinary inverse, the generalized inverse is not unique, in the sense that there are various ways of defining a generalized inverse with various inverse-like properties. The function matlib::Ginv() calculates a Moore-Penrose generalized inverse.

(AI <- Ginv (A))

## [,1] [,2] [,3]
## [1,] 0.27778 0 0.05556
## [2,] 0.00000 0 0.00000
## [3,] 0.05556 0 0.11111

We can also view this as fractions:

Ginv (A, fractions=TRUE)

## [,1] [,2] [,3]
## [1,] 5/18 0 1/18
## [2,] 0 0 0
## [3,] 1/18 0 1/9

Properties of generalized inverse (Moore-Penrose inverse)

The generalized inverse is defined as the matrix $A^-$ such that

$A * A^- * A = A$ and
$A^- * A * A^- = A^-$

A %*%  AI %*%  A

## [,1] [,2] [,3]
## [1,] 4 4 -2
## [2,] 4 4 -2
## [3,] -2 -2 10

AI %*%  A %*%  AI

## [,1] [,2] [,3]
## [1,] 0.27778 0 0.05556
## [2,] 0.00000 0 0.00000
## [3,] 0.05556 0 0.11111

In addition, both $A * A^-$ and $A^- * A$ are symmetric, but neither product gives an identity matrix, A %*% AI != AI %*% A != I

zapsmall (A %*%  AI)

## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 1 0 0
## [3,] 0 0 1

zapsmall (AI %*%  A)

## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 0 0 0
## [3,] 0 0 1

Rectangular matrices

For a rectangular matrix, $A^- = (A^{T} A)^{-1} A^{T}$ is the generalized inverse of $A$ if $(A^{T} A)^-$ is the ginv of $(A^{T} A)$ [See: Timm: EX 1.6.11]

A <- cbind ( 1, matrix (c (1, 0, 1, 0, 0, 1, 0, 1), nrow=4, byrow=TRUE))
A

## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 1 1 0
## [3,] 1 0 1
## [4,] 1 0 1

This $4 \times 3$ matrix is not of full rank, because columns 2 and 3 sum to column 1.

R (A)

## [1] 2

(AA <- t (A) %*%  A)

## [,1] [,2] [,3]
## [1,] 4 2 2
## [2,] 2 2 0
## [3,] 2 0 2

(AAI <- Ginv (AA))

## [,1] [,2] [,3]
## [1,] 0.5 -0.5 0
## [2,] -0.5 1.0 0
## [3,] 0.0 0.0 0

The generalized inverse of $A$ is $(A^{T} A)^- A^{T}$ , AAI * t(A)

AI <- AAI %*%  t (A)

Show that it is a generalized inverse:

A %*%  AI %*%  A

## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 1 1 0
## [3,] 1 0 1
## [4,] 1 0 1

AI %*%  A %*%  AI

## [,1] [,2] [,3] [,4]
## [1,] 0.0 0.0 0.5 0.5
## [2,] 0.5 0.5 -0.5 -0.5
## [3,] 0.0 0.0 0.0 0.0