Mathematics

Pseudoinverses

Rectangular matrices do not have inverses or determinants. At least one of the equations AX = I and XA = I does not have a solution. A partial replacement for the inverse is provided by the Moore-Penrose pseudoinverse, which is computed by the pinv function.

• C = fix(10*rand(3,2));
  X = pinv(C)
  
  X =
      0.0401   -0.1492    0.1050
      0.0110    0.1657   -0.0055
  

The matrix

• Q = X*C
  
  Q =
        1.0000    0.0000
        0.0000    1.0000
  

is the 2-by-2 identity, but the matrix

• P = C*X
  
  P =
      0.2044    0.0663    0.3978
      0.0663    0.9945   -0.0331
      0.3978   -0.0331    0.8011
  

is not the 3-by-3 identity. However, P acts like an identity on a portion of the space in the sense that P is symmetric, P*C is equal to C and X*P is equal to X.

Solving a Rank-Deficient System

If A is m-by-n with m > n and full rank n, then each of the three statements

• x = A\b
  x = pinv(A)*b
  x = inv(A'*A)*A'*b
  

theoretically computes the same least squares solution x, although the backslash operator does it faster.

However, if A does not have full rank, the solution to the least squares problem is not unique. There are many vectors x that minimize

• norm(A*x -b)
  

The solution computed by x = A\b is a basic solution; it has at most r nonzero components, where r is the rank of A. The solution computed by x = pinv(A)*b is the minimal norm solution because it minimizes norm(x). An attempt to compute a solution with x = inv(A'*A)*A'*b fails because A'*A is singular.

Here is an example that illustrates the various solutions.

• A = [ 1  2  3
        4  5  6
        7  8  9
       10 11 12 ]
  

does not have full rank. Its second column is the average of the first and third columns. If

• b = A(:,2)
  

is the second column, then an obvious solution to A*x = b is x = [0 1 0]'. But none of the approaches computes that x. The backslash operator gives

• x = A\b
  
  Warning: Rank deficient, rank = 2.  
  
  x =
        0.5000
        0
        0.5000
  

This solution has two nonzero components. The pseudoinverse approach gives

• y = pinv(A)*b
  
  y =
        0.3333
        0.3333
        0.3333
  

There is no warning about rank deficiency. But norm(y) = 0.5774 is less than norm(x) = 0.7071. Finally

• z = inv(A'*A)*A'*b
  

fails completely.

• Warning: Matrix is singular to working precision.
  
  z =
       Inf
       Inf
       Inf
  

Overview

Cholesky, LU, and QR Factorizations