rmatrixmatrix-inverse

Recover a matrix after matrix multiplication


I am given two matrices (A and B), where the second has rows with unit length. For example in R:

k <- 1000
m <- 10
n <- 100

A <- matrix(rnorm(m*k), m, k) # m x k
B <- matrix(rnorm(k*n), k, n) # k x n
B <- B/sqrt(rowSums(B^2)) # make rows have unit length

What is the best way to find a third matrix (C) such that the multiplication of the matrices ABC is as close as possible to the original matrix A? That is:

C <- MASS::ginv(B) # what is the best way to get C?

# [m x k] = [m x k] * [k x n] * [n x k]
A2 <- A %*% B %*% C

plot(A, A2, pch=46)
abline(a=0, b=1)
sum((A - A2)^2) # goal is to minimize this value

When B is full row rank (k <= n in the example of a random matrix above) the pseudoinverse gives nearly zero error. But when B is not full row rank (e.g., k > n) then the pseudoinverse (ginv(B)) is nearly as bad as the transpose (t(B)), especially when k >> n. Is the pseudoinverse always the best possible solution, even when k >> n?


Solution

  • Yes, pseudo inverse is the best possible solution.

    Obviously, the best C is B^T(B B^T)^{-1}, which is known as the right inverse. However this is only applicable in full row rank matrix.

    For k>n, the best we can do is (B^T B)^{-1} B^T, which is known as the left inverse. This is because B (B^T B)^-1 B^T is the projection matrix that projects onto the column space of B. This is indeed the pseudoinverse of B for this problem.

    C <- MASS::ginv(B)
    left_inverse <- solve(t(B) %*%B) %*% t(B)
    all.equal(left_inverse,C)
    [1] TRUE
    

    In addition, pseudoinverse is the best approximation for least square problem. I think the results can be extended to this problem.