I ran across a question that was posted on the Internet many years ago, but not truly answered. The question is, “Can you use QR decomposition to find the inverse of a matrix?” By truly answered, I mean a complete explanation and an actual demo of the technique.

The answer is yes, but I could not find a concrete example, with code, so I set out to put together such an example. The key ideas are as follows. If you have a matrix A and perform QR decomposition, you will get a matrix Q and a matrix R. If you can find the inverse of Q and R and multiply them together, you’ll get the inverse of A.

Because R is upper triangular, there is a very clever algorithm to find its inverse. See my post at https://jamesmccaffreyblog.com/2023/08/01/how-to-invert-an-upper-triangular-matrix-from-scratch-using-python/. Because of the special way Q is constructed, the inverse of Q is the transpose of Q.

A = QR
inv(A) = inv(QR)
       = inv(R) * inv(Q)
       = inv(R) * trans(Q)

Here’s my demo function to compute the inverse of a matrix using QR decomposition. My demo uses a from-scratch implementation of QR decomposition but it’s easier to use the scipy.linalg.qr() function.

def qr_inverse(A):
  q, r = my_qr(A, standard=True)  # could use scipy.linalg.qr()
  r_inv = upper_tri_inverse(r)  # invert R
  q_tr = np.transpose(q)  # invert Q (transpose)
  return np.matmul(r_inv, q_tr)  # the inverse of A

The “standard” parameter in the qr_inverse() function means to return Q and R in a standard form where the diagonal elements of R are all positive. This isn’t necessary. Somewhat annoyingly, the scipy.linalg.qr() function doesn’t have an option to return Q and R in standard form.

My demo code uses a clever but complicated, special algorithm to find the inverse of upper triangular matrix R. In a non-demo scenario you’d want to use a non-clever, standard algorithm:

def mat_inverse_upper(U):
  # inverse of upper triangular directly
  n = len(U)
  result = np.eye(n)  # Identity matrix
  for k in range(n):
    for j in range(n):
      for i in range(k):
        result[j][k] -= result[j][i] * U[i][k]
      result[j][k] /= U[k][k]
  return result

Two images of artsy mathematical matrices I generated using the DALL-E text-to-image tool. I’m clearly not an artist, even with the help of AI. I’ll stick to mathematics and computer science.

Demo program. Replace “lt” with less-than Boolean operator symbols.

# qr_matrix_inverse.py

# find the inverse of a matrix using QR decomposition

import numpy as np

# -----------------------------------------------------------

def qr_inverse(A):
  q, r = my_qr(A, standard=True)  # could scipy.linalg.qr()
  r_inv = upper_tri_inverse(r)  # invert R
  q_tr = np.transpose(q)  # invert Q (transpose)
  return np.matmul(r_inv, q_tr)  # the inverse of A

# -----------------------------------------------------------

def my_qr(A, standard=False):
  # QR decomposition. standard: make R diag positive
  # Householder algorithm

  m, n = A.shape  # check m == n
  Q = np.eye(n)
  R = np.copy(A)
  for i in range(n-1):
    H = np.eye(m)

    # -------------------
    a = R[i:, i]  # extract partial column
    norm_a = np.linalg.norm(a)  # 
    if a[0] "lt" 0.0: norm_a = -norm_a
    v = a / (a[0] + norm_a)
    v[0] = 1.0
    h = np.eye(len(a))  # H reflection
    h -= (2 / np.dot(v, v)) * \
      np.matmul(v[:,None], v[None,:])  # WTF. thank you np
    # -------------------

    H[i:, i:] = h  # copy h to low right corner H
    Q = np.matmul(Q, H)
    R = np.matmul(H, R)

  if standard == True:
    S = np.zeros((n,n))  # signs of diagonal
    for i in range(n):
      if R[i][i] "lt" 0.0: S[i][i] = -1.0
      else: S[i][i] = +1.0
    Q = np.matmul(Q, S)
    R = np.matmul(S, R)

  return Q, R

# -----------------------------------------------------------

def upper_tri_inverse(A):
  # A is upper triangular -- values below diagonal all zero
  # clever special case algorithm
  n = len(A) 

  Di = np.zeros((n,n))  # diag element inverses
  for i in range(n):
    Di[i][i] = 1 / A[i][i]

  Tu = np.zeros((n,n))  # strictly upper triangular
  for i in range(0,n):
    for j in range(0,n):
      if i "lt" j:
        Tu[i][j] = A[i][j]
  
  DiTu0 = np.eye(n)  # math magic
  DiTu1 = np.matmul(-Di, Tu)  # more magic

  sum = DiTu0 + DiTu1
  prev = DiTu1
  for i in range(2,n):
    curr = np.matmul(DiTu1, prev)
    sum += curr
    prev = curr
  
  inv = np.matmul(sum, Di)
  return inv

# -----------------------------------------------------------

def main():
  print("\nBegin use QR decomp to find matrix inverse ") 
  np.set_printoptions(precision=8, suppress=True)

  A = np.array([[1, 2, 3, 4],
                [6, 5, 6, 7],
                [8, 9, 8, 9],
                [2, 4, 3, 7]], dtype=np.float64) 

  print("\nSource matrix A = ")
  print(A)

  A_inv = qr_inverse(A)
  print("\nInverse A_inv = ")
  print(A_inv)

  # verify inverse gives I matrix
  AiA = np.matmul(A_inv, A)
  print("\nA_inv * A = ")
  print(AiA)  # should be I matrix

  print("\nEnd demo ")

if __name__ == "__main__":
  main()