# Avoiding computing inverse of matrix, don't think it's possible in my case

I have a system of (6+3N) equations given by:

`A*inv(B)*C*dX = K`

which I then solve using

`dX = (A*inv(B)*C) \ K`

A*inv(B)*C form a square (6+3N) matrix but the submatrices have sizes:

A: (6+3N) x (6N)
B: (6N) x (6N)
C: (6N) x (6+3N)

The bulk of time taken to solve this sytem of equations is spent calculating the inverse of B. I know generally it is advised to avoid calculating inverses of matrices and transform the problem into a solution of a linear system. However, in this case, I don’t think it’s possible to do this? I just wanted a quick opinion to see if there’s anything clever I can do here. Cheers

What’s wrong with `dX = (B \ A * C) \ K`

2 Likes

Because A, B and C are different sizes, B \ A * C cannot be done sadly. If they were all square then I think this would be much simpler.

``````julia> size.((A,B,C,K))
((66, 120), (120, 120), (120, 66), (66,))

julia> B \ A*C
ERROR: DimensionMismatch: arguments must have the same number of rows
Stacktrace:
 \(F::LU{Float64, Matrix{Float64}, Vector{Int64}}, B::Matrix{Float64})
@ LinearAlgebra C:\Users\hj20005\AppData\Local\Programs\Julia-1.8.3\share\julia\stdlib\v1.8\LinearAlgebra\src\LinearAlgebra.jl:472
 \(A::Matrix{Float64}, B::Matrix{Float64})
@ LinearAlgebra C:\Users\hj20005\AppData\Local\Programs\Julia-1.8.3\share\julia\stdlib\v1.8\LinearAlgebra\src\generic.jl:1110
 top-level scope
@ REPL:1
``````

Oops, sorry. It’s supposed to be `(A/B * C) \ K` (or `A * (B \ C)) \ K` depending on how you want to compute it).

4 Likes

Amazing, thanks. I clearly don’t have a good understanding of the relevant linear algebra here…will read up a bit…

Straightforward >2x speedup for free! I expect for larger system sizes the speedup will be even greater!

``````@btime dX1 = (A*inv(B)*C) \ K
@btime dX2 = (A/B * C) \ K
@btime dX3 = (A * (B \ C)) \ K

40.912 ms (13 allocations: 12.22 MiB)
17.451 ms (13 allocations: 12.21 MiB)
15.032 ms (11 allocations: 11.63 MiB)
``````