Computing `AB⁻¹` with `np.linalg.solve()`

Question

I need to compute AB⁻¹ in Python / Numpy for two matrices A and B (B being square, of course).

I know that np.linalg.inv() would allow me to compute B⁻¹, which I can then multiply with A. I also know that B⁻¹A is actually better computed with np.linalg.solve().

Inspired by that, I decided to rewrite AB⁻¹ in terms of np.linalg.solve(). I got to a formula, based on the identity (AB)ᵀ = BᵀAᵀ, which uses np.linalg.solve() and .transpose():

np.linalg.solve(a.transpose(), b.transpose()).transpose()

that seems to be doing the job:

import numpy as np


n, m = 4, 2
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))

print(np.matmul(b, np.linalg.inv(a)))
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
#  [-1.08733434  1.00110176  0.79683577  0.67487591]]
print(np.linalg.solve(a.transpose(), b.transpose()).transpose())
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
#  [-1.08733434  1.00110176  0.79683577  0.67487591]]
print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True

and also comes up much faster for sufficiently large inputs:

n, m = 400, 200
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))

print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True

%timeit np.matmul(b, np.linalg.inv(a))
# 100 loops, best of 3: 13.3 ms per loop
%timeit np.linalg.solve(a.transpose(), b.transpose()).transpose()
# 100 loops, best of 3: 7.71 ms per loop

---

My question is: does this identity always stand correct or there are some corner cases I am overlooking?

Answer 1

In general, np.linalg.solve(B, A) is equivalent to B-1A. The rest is just math.

In all cases, (AB)T = BTAT: https://math.stackexchange.com/q/1440305/295281.

Not necessary for this case, but for invertible matrices, (AB)-1 = B-1A-1: https://math.stackexchange.com/q/688339/295281.

For an invertible matrix, it is also the case that (A-1)T = (AT)-1: https://math.stackexchange.com/q/340233/295281.

From that it follows that (AB-1)T = (B-1)TAT = (BT)-1AT. As long as B is invertible, you should have no issues with the transformation you propose in any case.