There are tons of Internet resources that go something like, “blah, blah, blah and so matrix R is upper triangular and so it’s easy to find its inverse.” But then there is either no further information about how to invert the upper triangular matrix, or there is not particularly helpful information, such as, “Use back-substitution”.

Note: It’s perfectly possible to invert an upper triangular matrix using a standard approach but I came across a clever algorithm that intrigued me. Here’s code for the standard, non-clever approach:

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

Here’s Python/NumPy code to invert the special case of an upper triangular matrix, using the clever algorithm:

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

  Di = np.zeros((n,n))  # diag scalar 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

The function can be called like:

A = np.array([[1, 2, 3],
              [0, 4, 5],
              [0, 0, 6]])

Ainv = upper_tri_inverse(A)

This code is an implementation of a very clever algorithm posted by a Robert Lewis in a Math Stack Exchange answer from 2014. See https://math.stackexchange.com/questions/1003801/inverse-of-an-invertible-upper-triangular-matrix-of-order-3. Expressed as a math equation the algorithm is:

A is the source upper triangular matrix to invert. Lambda-1 is a diagonal matrix of the scalar inverses of the diagonal elements of A. Tu is the strictly upper triangular elements of A. The * is matrix multiplication. I is the identity matrix.

I use explicit looping (which is inefficient in Python) so that I can easily refactor my code to a language like C#. I omit all normal error checking for clarity.

An alternative implementation of inverse for an upper triangular matrix saves all the intermediate matrices so you can debug easier:

def upper_tri_inverse(A):
  # A is upper triangular -- values below diagonal all zero
  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

  lst = []  # save intermediate terms for debugging
  lst.append( DiTu0 )  # [0]
  lst.append( DiTu1 )  # [1]

  sum = lst[0] + lst[1]
  for i in range(2,n):
    lst.append( np.matmul (DiTu1, lst[i-1] )) 
    sum += lst[i]
  
  inv = np.matmul(sum, Di)
  return inv

As it turns out, if you know how to invert an upper triangular matrix, you can use such a function in conjunction with QR matrix decomposition to find the inverse of an arbitrary matrix. Briefly, if you want to find the inverse of a matrix A:

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

The idea is to decompose A into Q and R. R is upper triangular, so you can find the inverse of R using the code in this post. Because of the special way Q is constructed, the inverse of Q is just the transpose of Q. If you multiply the two inverse together, you get the inverse of A. Put another way, inverting an arbitrary matrix is extremely difficult, but QR decomposition is only moderately difficult, and inverting an upper triangular matrix is only moderately difficult.

Very nice idea. I’ll post an example of inverting a matrix using QR decomposition plus upper triangular inverse when I get a chance.

An upper triangular matrix is a special, simple case of a matrix. I’ve always thought that black and white photography is sort of a special, simple case of an image of reality. Here are three such special cases with a theme of a face matrix. I fully understand math matrices; I do not fully understand face matrices.

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

# invert_triangular.py

# clever algorithm to invert an upper triangular matrix
# https://math.stackexchange.com/questions/
#  1003801/inverse-of-an-invertible-upper-
#  triangular-matrix-of-order-3 

import numpy as np

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

def upper_tri_inverse(A):
  # A is upper triangular -- values below diagonal all zero
  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("\nInverse upper triangular matrix from scratch ")
  np.set_printoptions(precision=8, suppress=True)

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


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

  Ai = upper_tri_inverse(A)
  print("\nInverse Ai: ")
  print(Ai)

  AiA = np.matmul(Ai, A)
  print("\nAi * A = ")
  print(AiA)  # should be Identity

  AAi = np.matmul(A, Ai)
  print("\nA * Ai = ")
  print(AAi)   # should be Identity

  print("\nEnd demo ")

if __name__ == "__main__":
  main()