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

General Usage Performance
,
1

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

2

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

2 Likes
3

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:
 [1] \(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
 [2] \(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
 [3] top-level scope
   @ REPL[46]:1
4

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
5

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)