Type inference and exponentially large compile times

[fguevaravas]

Issue #97 seemed to fix problems with type inference taking an inordinate amount of time when stringing together several LinearOperators. As far as I could tell, issue #97 was fixed, but the problem seems to be back in Julia 1.8.1. Please see below for a minimal working example + output on a fresh session of Julia 1.8.1 (I am using v2.4.1). Of course after running this another time in the same Julia session, the type inference is done and the whole thing takes very little time (as it would be expected)

# put in file lops.jl
# LinearOperator type stability
using LinearOperators
Ns = 4:17
function test_type_stability(N,n=5)
  lops = [ LinearOperator(randn(n,n)) for i=1:N ]
  C = hcat(lops...)
  stats = @timed y = C*randn(n*length(lops))
  return stats.time
end
ts = test_type_stability.(Ns)
for i=1:length(Ns)
 println("t=$(ts[i]), N=$(Ns[i])")
end

Julia takes 46s to do type inference on a hcat of 17 operators:

julia> include("lops.jl")
t=0.160799892, N=4
t=0.063678277, N=5
t=0.053124796, N=6
t=0.052513037, N=7
t=0.061400427, N=8
t=0.061368384, N=9
t=0.072215448, N=10
t=0.11183163, N=11
t=0.242638285, N=12
t=0.639014602, N=13
t=1.658588125, N=14
t=5.229325223, N=15
t=14.906830724, N=16
t=46.837328005, N=17

[dpo]

Thank you for the report. Would an implementation of hcat and vcat for an arbitrary number of operators solve the issue? (As opposed to one that concatenates ones operator at a time.)

[fguevaravas]

Yes would solve the issue with hcat and cat in my particular usage case (and I ended up doing something similar as a workaround).

Unfortunately, the problem will still persist for any operation that combines several linear operators together. For example the following:

A = opEye(n)
for i=1:N
 A = A*LinearOperator(randn(n,n))
end
y = A*randn(n)

takes way too long inferring types for say N>=15. This is not a surprise considering that a simple product:

A = LinearOperator(randn(10,10))
B = LinearOperator(randn(10,10))
C = A*B
typeof(C)

gives a long type that is not included here because it is not human-readable. BTW, when developing code using LinearOperators, this causes error logs that are sometimes so long that they overrun my terminal's buffer (and not informative... but that is more of a Julia issue)

As far as I understand, this is a consequence of using @closure?

[geoffroyleconte]

Yes, this is probably because the parametric type LinearOperator has 3 of its parameters F, Ft, Fct that are the types of prod!, tprod!, ctprod!, and when using closures these types can be really long:

LinearOperators.jl/src/abstract.jl

Lines 47 to 63 in a204372

	mutable struct LinearOperator{T, I <: Integer, F, Ft, Fct, S} <: AbstractLinearOperator{T}
	nrow::I
	ncol::I
	symmetric::Bool
	hermitian::Bool
	prod!::F
	tprod!::Ft
	ctprod!::Fct
	nprod::I
	ntprod::I
	nctprod::I
	args5::Bool
	use_prod5!::Bool # true for 5-args mul! and for composite operators created with operators that use the 3-args mul!
	Mv5::S
	Mtu5::S
	allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
	end

After looking at #97, I see that the fix was to remove the parameters F, Ft, Fct (that have been reintroduced with v2.0), but I am not sure that removing these parameters is the best solution for performance.

[fguevaravas]

I really don't see the benefit of using closures. I do see however that using closures affects performance the first time a LinearOperator is applied (and very significantly for complicated combinations of operators).

Maybe it is good to add warning somewhere in the documentation about this limitation, i.e. operations involving more than 10 operators can be long to compile. In Julia, compilation is also a performance problem...

FYI, I just realized that LinearMaps.jl does not have the same problem. The types generated by composite operations are also more humanly readable. Anyway, I have some testing to do to see if it is worth migrating.